Class Specsyn

Inheritance Relationships

Derived Types

Class Documentation

class Specsyn

A base class defining a common interface for spectral synthesis.

This class defines a common interface to spectral synthesis classes, whose purpose is to take as input a set of stellar parameters produced by interpolation on the tracks, and return as output stellar spectra on a pre-defined wavelength grid.

Subclassed by specsyn::SpecsynBlackbody, specsyn::SpecsynLib< Policy >, specsyn::SpecsynLibChained

Public Types

using Segment = interp::Interpolator1D<static_cast<size_t>(tracks::FieldIdx::nTrackQty)>
using Isochrone = std::vector<std::unique_ptr<Segment>>
using StarData = std::array<double, static_cast<size_t>(tracks::FieldIdx::nTrackQty)>

Public Functions

inline explicit Specsyn(const io::SimControls &controls)

Construct a Specsyn.

Parameters:

controls – Simulation controls this synthesizer reads its integrator tolerances (intRelTol(), intAbsTol(), intMaxIter()) and redshift (see wlObs()) from, live, for the rest of its lifetime – see controls_’s own comment. Must outlive this Specsyn.

virtual ~Specsyn() = default
Specsyn(const Specsyn&) = default
Specsyn(Specsyn&&) = default
auto operator=(const Specsyn&) -> Specsyn& = delete
auto operator=(Specsyn&&) -> Specsyn& = delete
inline auto wl() const -> const std::vector<double>&

Return the rest-frame wavelength grid.

Returns:

A const reference to the wavelength grid, in Angstrom

inline auto controls() const -> const io::SimControls&

Get the SimControls this Specsyn was constructed against.

Returns:

A const reference to controls_ – see its own comment. Exposed so a caller replacing a SimControls’s own specsyn (see SimControls::setSpecsyn()) can verify it is installing a Specsyn actually built against that same SimControls, rather than one whose controls_ points somewhere else entirely.

auto wlObs() const -> std::vector<double>

Return the observed-frame wavelength grid.

Defined out-of-line, in Specsyn.cpp – see that file’s own comment for why (same reason as intRelTol()/intAbsTol()/ intMaxIter(): controls_.z() needs io::SimControls’s complete type, which this header can only forward-declare).

Returns:

The wavelength grid, in Angstrom, redshifted by (1 + z), with z read live from controls_

inline virtual auto fehMin() const -> double

Return the minimum [Fe/H] this synthesizer has real spectral data for.

Unlike getMin() on a pdfs::PDF, this describes a hard data boundary, not a probability domain: a synthesizer with no [Fe/H] axis at all (SpecsynBlackbody, SpecsynLibWD) has nothing to restrict, so the base class default of -infinity (paired with fehMax()’s own +infinity) leaves any feh unclamped. Exists for SpecsynLibChained’s own benefit, which scans every library it chains together via this (see its own fehMin_ member) to build a per-GridType [Fe/H] range, ultimately used by its own specForIntegration() override – see that function’s own comment for why a single library’s own fehMin()/fehMax() aren’t clamped against directly by continuousSpecIntegrand() (which doesn’t, in general, know which concrete Specsyn it’s talking to).

Returns:

-infinity by default (no restriction) – overridden by SpecsynLibNoWind/SpecsynLibWR (their own FeH_.front())

inline virtual auto fehMax() const -> double

Return the maximum [Fe/H] this synthesizer has real spectral data for.

Returns:

+infinity by default (no restriction) – see fehMin()’s own comment for the overriding classes and full rationale, which applies here symmetrically

inline auto requestedFehMin() const -> double

Return the minimum [Fe/H] this synthesizer was constructed to cover.

Unlike fehMin(), which reports the real data a library holds, this is the range the caller asked for, fixed at construction. SimControls::setSpecsyn() compares it against the currently installed synthesizer’s own, to reject a replacement built for a narrower [Fe/H] range.

Returns:

The fehMin its constructor was given, or -infinity for a synthesizer with no [Fe/H] axis (SpecsynBlackbody, SpecsynLibWD)

inline auto requestedFehMax() const -> double

Return the maximum [Fe/H] this synthesizer was constructed to cover.

Returns:

The fehMax its constructor was given, or +infinity for a synthesizer with no [Fe/H] axis – see requestedFehMin()’s own comment

inline virtual auto loggMin() const -> double

Return the minimum log(g) this synthesizer has real spectral data for.

Mirrors fehMin()’s own rationale exactly, for SpecsynLibChained’s per-library log(g) clamp rather than its per-library [Fe/H] one – see clampFehForLibrary()’s own comment for why a star whose own log(g) falls outside every chained library’s real coverage needs the same kind of per-library rescue a star outside every chained library’s [Fe/H] coverage already gets.

Returns:

-infinity by default (no restriction) – overridden by SpecsynLibNoWind/SpecsynLibWD (their own logg().front()); SpecsynLibWR keeps the base default, since Wolf-Rayet atmospheres aren’t parameterized by log(g) at all

inline virtual auto loggMax() const -> double

Return the maximum log(g) this synthesizer has real spectral data for.

Returns:

+infinity by default (no restriction) – see loggMin()’s own comment for the overriding classes and full rationale, which applies here symmetrically

virtual auto spec(const StarData &props, double feh) const -> std::vector<double> = 0

Compute the spectrum of a single star.

Parameters:

  • props – Stellar properties, as produced by evaluating the Interpolator1D returned by Tracks2D::getIsochrone at this star’s mass

  • feh – [Fe/H] value of the star, needed because it is not carried by props itself (e.g. by a SpecsynLib, to locate props in its spectral library’s [Fe/H] direction)

Returns:

The star’s spectrum, evaluated on the wavelength grid returned by wl(), in units of erg/s/Angstrom

inline virtual auto specForce(const StarData &props, double feh) const -> std::vector<double>

Compute a star’s spectrum, forcing a result even if it falls outside this Specsyn’s own domain.

Unlike spec(), which may return an empty vector (depending on OOBPolicy) for a star outside this Specsyn’s own domain, specForce() always either returns a genuine spectrum or throws. Used by SpecsynLibChained as a last resort on the last library in a star’s own type-specific chain, once every library in that chain (tried via ordinary spec()) has failed to produce one.

The default implementation here just calls spec() and throws if it comes back empty – appropriate for a Specsyn (like SpecsynLibWR) whose own spec() already forces a match internally (e.g. by moving a star to the nearest populated grid point along one axis) rather than ever returning empty for a star within its own domain. SpecsynLibNoWind and SpecsynLibWD override this with a real forcing mechanism of their own (moving log(g) to the nearest populated grid value), since their own spec() can genuinely return empty for a star whose log(g) falls in a gap ordinary interpolation can’t bridge.

Parameters:

  • props – Stellar properties, as produced by evaluating the Interpolator1D returned by Tracks2D::getIsochrone at this star’s mass

  • feh – [Fe/H] value of the star, needed because it is not carried by props itself

Throws:

std::runtime_error – if no spectrum can be produced for this star even by forcing

Returns:

The star’s spectrum, evaluated on the wavelength grid returned by wl(), in units of erg/s/Angstrom – unlike spec(), never a size-0 vector

inline virtual auto specForIntegration(const StarData &props, double feh) const -> std::vector<double>

Compute a star’s spectrum for Specsyn::continuousSpecIntegrand()’s own benefit.

The default implementation here just calls spec(props, feh) unchanged – correct for every Specsyn with a single, uniform notion of its own [Fe/H] coverage (SpecsynBlackbody, SpecsynLibWD, and a standalone SpecsynLibNoWind/SpecsynLibWR alike), which is exactly what spec()/specForce() already handle via their own existing OOBPolicy. Overridden only by SpecsynLibChained, whose own [Fe/H] coverage is not uniform (it varies by which chained library, of possibly several different real coverages, a given star’s own GridType dispatches to) – see its own override’s comment for why that needs a real clamp, and why this exists as a separate function from spec()/specForce() rather than changing their own behavior: [Fe/H] deliberately has no forcing/rescue mechanism there (unlike log(Teff)/log(g)), since a genuinely out-of-range [Fe/H] request reflects a real physical difference no nearest-neighbor substitution can paper over – a distinction this function exists specifically to preserve for every caller except Specsyn::continuousSpecIntegrand() (via specWlForIntegration() below), which evaluates at whatever [Fe/H] values Specsyn::specCtsHelper()’s own [Fe/H] quadrature chooses, and so still needs a usable spectrum there even where one chained library’s own real coverage stops short.

Parameters:

  • props – See spec()’s own props parameter

  • feh – [Fe/H] value of the star – unlike spec()’s/ specForce()’s own feh, not necessarily guaranteed to fall within this synthesizer’s own real data coverage; see this function’s own

Returns:

The star’s spectrum – see spec()’s own return value

inline auto spec(double m, const Segment &segment, double feh) const -> std::vector<double>

Compute the spectrum of a single star, given its mass and isochrone segment.

Evaluates segment at m to get m’s stellar properties, and returns spec() of those properties. This overload exists mainly so spec() can be handed to PDFIntegrator, which expects a callable taking the integration variable (here, mass) as its first argument. It takes a single segment, rather than a full Isochrone, so that specCts() can integrate each segment separately over its own valid domain – an Isochrone as a whole may have gaps where no segment is defined, and the quadrature routine has no way to know to avoid evaluating there.

Parameters:

  • m – Stellar mass, in Msun; must lie within segment’s valid domain (segment.xMin() <= m <= segment.xMax())

  • segment – A single isochrone segment (one element of the Isochrone returned by Tracks2D::getIsochrone) to evaluate at mass m

  • feh – [Fe/H] value of the segment’s isochrone, passed through to spec() unchanged

Returns:

The star’s spectrum, evaluated on the wavelength grid returned by wl(), in units of erg/s/Angstrom

inline auto specWl(double m, const Segment &segment, double feh) const -> std::vector<double>

Compute the wavelength-weighted spectrum of a single star, given its mass and isochrone segment.

Exists so specCts() can integrate lambda * dL/dlambda instead of dL/dlambda: near the peak of a stellar population’s SED, lambda * dL/dlambda is order-unity in reasonable physical units (comparable to the population’s total luminosity), whereas dL/dlambda itself spans many orders of magnitude across a wide wavelength grid, making a single absolute tolerance on it effectively meaningless.

Parameters:

  • m – Stellar mass, in Msun; must lie within segment’s valid domain (segment.xMin() <= m <= segment.xMax())

  • segment – A single isochrone segment (one element of the Isochrone returned by Tracks2D::getIsochrone) to evaluate at mass m

  • feh – [Fe/H] value of the segment’s isochrone, passed through to spec() unchanged

Returns:

spec(m, segment, feh), multiplied elementwise by wl(): lambda * dL/dlambda rather than dL/dlambda, in units of erg/s

inline auto specWlForIntegration(const StarData &props, double feh) const -> std::vector<double>

Compute the wavelength-weighted spectrum of a single star, for Specsyn::continuousSpecIntegrand()’s own benefit.

Identical to specWl(), except taking a star’s properties directly (rather than a mass and an isochrone segment to evaluate) and calling specForIntegration() in place of spec() – see specForIntegration()’s own comment for why Specsyn::continuousSpecIntegrand() needs this, rather than specWl() itself.

Parameters:

  • props – See spec()’s own props parameter

  • feh – See specForIntegration()’s own feh parameter – unlike specWl()’s own feh, not necessarily within this synthesizer’s own real data coverage

Returns:

specForIntegration(props, feh), multiplied elementwise by wl() – see specWl()’s own return value

virtual auto specCts(const Isochrone &isochrone, const pdfs::PDF &imf, double mTot, double mMin, double mMax, double feh) const -> std::vector<double>

Compute the spectrum of a continuously-sampled stellar population.

Combines the star-by-star spectral synthesis provided by spec() with the population’s IMF, continuously sampled over [mMin, mMax] rather than drawn as a discrete set of stars.

isochrone may consist of several disjoint segments (see Tracks2D::getIsochrone), each valid only over its own [segment.xMin(), segment.xMax()] range, with gaps possibly in between where no segment is defined. The quadrature routine has no way to know to avoid evaluating in those gaps, so this integrates each segment separately, over the intersection of its own domain with [mMin, mMax], skipping segments whose domain does not overlap [mMin, mMax] at all, and sums the per-segment results – this guarantees spec() is only ever evaluated at masses where some segment is actually defined. The summed result is then scaled by mTot.

Parameters:

  • isochrone – The isochrone for the population, as returned by Tracks2D::getIsochrone

  • imf – The initial mass function of the population

  • mTot – Total mass of the population, in Msun – i.e. of its stars in [mMin, mMax]; the integral of the per-star spectrum against imf over [mMin, mMax] is scaled by mTot divided by the integral of m imf(m) dm over that same range, since imf is normalized by number, not mass

  • mMin – Minimum stellar mass in the population, in Msun

  • mMax – Maximum stellar mass in the population, in Msun

  • feh – [Fe/H] value of the population, passed through to spec() unchanged

Returns:

The specific luminosity of the population, evaluated on the wavelength grid returned by wl(), in units of erg/s/Angstrom

auto specCts(const pdfs::PDF &sfr, const pdfs::PDF &imf, const pdfs::PDF &fehDist, double curTime, double fCluster, double mMin, double mMax) const -> std::vector<double>

Compute the integrated spectrum of a continuously-formed stellar population.

Computes the integrated spectrum from a stellar population that is continuously distributed in age from 0 to curTime, in mass over imf, and in metallicity over fehDist – unlike the other specCts() overload, which integrates a single fixed-age isochrone over mass alone, this one integrates over age (equivalently, formation time) and metallicity as well, so a different isochrone applies at different points in the integration domain, computed on demand – see continuousSpecIntegrand()’s own comment for how.

Called by Galaxy::computeSpec() when fCluster() < 1, to cover the star-forming mass Galaxy’s own individual Cluster objects don’t – see io::SimControls::fCluster()’s own comment.

Structured as a pair of nested 1D integrals, not a single joint 2D (age, mass) one: a utils::PDFIntegrator (GKOrder::GK15) over age alone, whose own integrand – continuousSpecIntegrand() – builds the isochrone for its own age on demand and immediately hands it to specCtsImpl() to perform a complete inner 1D integral over mass, reusing the same proven mass integrator the isochrone specCts() overload itself uses. An earlier design instead integrated (age, mass) jointly as a single 2D domain; that domain has an internal, curved live/dead-mass boundary (wherever the shortest-lived stars at a given age reach the end of their own lives) that 2D adaptive cubature – h- or p-adaptive alike – resolves very poorly once curTime exceeds those stars’ own lifetimes, producing a sharp, empirically observed cost cliff (some real, permanent-test output times converging in seconds, others not converging within minutes). Nesting the integrals instead removes that boundary from the integration domain entirely: each age’s own isochrone already encodes exactly which masses are alive, as ordinary segment edges specCtsImpl() naturally integrates up to and no further – see continuousSpecIntegrand()’s own comment.

If fehDist is degenerate (fehDist.getMin() == fehDist.getMax(), a single metallicity shared by every star), runs this nested integral once. Otherwise, integrates over [Fe/H] with a third, outermost 1D PDFIntegrator weighted by fehDist, whose integrand is the complete nested (age, mass) integral at each [Fe/H] it visits – see specCtsHelper()’s own comment.

The time dimension is integrated as age instead, via a pdfs::PDFReflect view of sfr pivoted at curTime / 2 (so reflect(x) = curTime - x, i.e. this coordinate already is age), and log-transformed whenever curTime exceeds min(1e4 yr, 1e-3 * curTime) – see specCtsHelper()’s own comment for why: the luminosity this integral is dominated by (young, hot, blue stars) is concentrated in a narrow sliver of age near 0, which a linear axis resolves increasingly poorly as curTime grows.

As with the other specCts() overload, the actual quantity integrated is lambda * dL/dlambda (via continuousSpecIntegrand(), which mirrors specWl()’s own multiplication), undone by dividing back out by wl_ elementwise once the integral itself is complete, before scaling by (1 - fCluster) * (1 - controls_.fracStochMass()): fCluster is the fraction of star-forming mass diverted to individual, stochastically-sampled Cluster objects (see fCluster’s own comment); of the (1 - fCluster) that remains, mMin/mMax already restrict this call’s own integral to whatever mass range the caller asked for (e.g. Galaxy::computeSpec() passes [imf.getMin(), controls_.minStochMass()], the purely continuous, sub-minStochMass share of the non-clustered population – see its own comment), but controls_.fracStochMass() – the IMF-integral fraction of mass at or above controls_.minStochMass() – is what the caller’s own mMax choice implicitly excludes, so it must be divided back out of the (1 - fCluster) share explicitly here, rather than folding it into mMax/mTot the way the isochrone specCts() overload’s own mTot already does. In the historical case where mMax >= imf.getMax() (mMin/mMax cover the whole IMF, i.e. every non-clustered star is treated continuously, with no field-star population above minStochMass at all), fracStochMass() is exactly 0, and this factor is a no-op. Finally, the result is divided by the integral of m imf(m) dm over [mMin, mMax]: imf is normalized by number, so the integral against it is per star, and this converts it into one per unit mass of stars in [mMin, mMax], matching the mass-based (1 - fCluster) * (1 - fracStochMass()) share it is scaled by – mirroring the isochrone specCts() overload’s own mTot scaling.

A thin wrapper over specCtsHelper() (computeLbol = false) – see specAndLbolCts() for the sibling that also returns the same population’s bolometric luminosity, computed alongside the spectrum in the same integral rather than a second one.

Parameters:

  • sfr – Star formation rate, in Msun/yr, as a function of time (see io::SimControls::sfr()) – unlike a strictly normalized probability density, integrating this over a time range gives the actual stellar mass (in Msun) formed over that range, matching how io::SimControls::sfr() is already used elsewhere (e.g. Galaxy::advance())

  • imf – The initial mass function of the population

  • fehDist – The [Fe/H] distribution of the population

  • curTime – The current time, in yr; star formation is integrated from 0 to curTime

  • fCluster – The fraction of star-forming mass treated stochastically, in individual Cluster objects (see io::SimControls::fCluster()); the result is scaled by (1 - fCluster) * (1 - controls_.fracStochMass()), divided by the integral of m imf(m) dm over [mMin, mMax] – see this function’s own comment for the second and third factors

  • mMin – Minimum stellar mass to integrate over, in Msun – see the isochrone specCts() overload’s own mMin parameter; passed straight through to the inner mass integral (via specCtsImpl()), exactly like that overload’s own mMin

  • mMax – Maximum stellar mass to integrate over, in Msun – see mMin’s own comment

Returns:

The specific luminosity of the continuously-formed population, evaluated on the wavelength grid returned by wl(), in units of erg/s/Angstrom

auto specAndLbolCts(const pdfs::PDF &sfr, const pdfs::PDF &imf, const pdfs::PDF &fehDist, double curTime, double fCluster, double mMin, double mMax) const -> std::pair<std::vector<double>, double>

Compute the integrated spectrum and bolometric luminosity of a continuously-formed stellar population, together.

Identical to specCts() in every respect but what it returns: calls specCtsHelper() with computeLbol = true, so the same PDFIntegrator evaluation – the expensive part of which is building an isochrone at every distinct (age, feh) the integral visits, not the handful of extra scalar quantities integrated alongside the spectrum at each of those isochrones – also integrates the population’s bolometric luminosity, rather than requiring a second, separate integral (sharing none of the same isochrones) to get it. See continuousSpecIntegrand()’s own comment for exactly how the extra quantity is folded into the integrand itself.

Used by Galaxy::computeSpec() in place of specCts() itself when io::SimControls::computeLbol() is true, so that a caller requesting both a spectrum and Lbol together (the common case) never pays for isochrone construction twice.

Parameters:

  • sfr – See specCts()’s own sfr parameter

  • imf – See specCts()’s own imf parameter

  • fehDist – See specCts()’s own fehDist parameter

  • curTime – See specCts()’s own curTime parameter

  • fCluster – See specCts()’s own fCluster parameter

  • mMin – See specCts()’s own mMin parameter

  • mMax – See specCts()’s own mMax parameter

Returns:

A pair of (1) the same specific luminosity specCts() itself returns, and (2) the population’s own bolometric luminosity, in Lsun (matching Cluster::lbol()’s own units)

auto intRelTol() const -> double

Return the relative tolerance for PDF integration.

Returns:

Relative tolerance passed to PDFIntegrator, read live from controls_

auto intAbsTol() const -> double

Return the absolute tolerance for PDF integration.

Returns:

Absolute tolerance passed to PDFIntegrator, read live from controls_

auto intMaxIter() const -> std::size_t

Return the maximum number of evaluations for PDF integration.

Returns:

Max bisection iterations passed to PDFIntegrator (0 = unlimited), read live from controls_

Protected Functions

inline void setRequestedFehRange(const double fehMin, const double fehMax)

Record the [Fe/H] range this synthesizer was constructed to cover.

Called by every derived constructor that takes an [Fe/H] range; see requestedFehMin()’s own comment.

Parameters:

  • fehMin – The fehMin its constructor was given

  • fehMax – The fehMax its constructor was given

Protected Attributes

const io::SimControls &controls_

Simulation controls this synthesizer reads its integrator tolerances and redshift from.

Read live, not snapshotted, every time specCts() integrates or wlObs() is called (see intRelTol()/intAbsTol()/ intMaxIter()/wlObs(), which just forward to the same-named methods on this reference) – changing controls_’s own tolerances or redshift after this Specsyn is built takes effect immediately, with no need to rebuild the synthesizer. Bound once, at construction, from whichever SimControls actually built this synthesizer (see SimControls::readSpectra()); never reseated afterward, so that SimControls must outlive this Specsyn.

std::vector<double> wl_

Wavelength grid for the spectral synthesizer, in Angstrom

Protected Static Functions

static auto describeStar(const StarData &props, double feh) -> std::string

Describe a star’s properties, for use in an error message.

Parameters:

  • props – Stellar properties, as produced by evaluating the Interpolator1D returned by Tracks2D::getIsochrone at this star’s mass

  • feh – [Fe/H] value of the star

Returns:

A one-line, human-readable listing of feh and every field of props, by name, at full precision

static inline auto getSAandLogg(const StarData &props) -> std::pair<double, double>

Compute a star’s surface area and log(g) from its stellar properties.

The stellar radius (and hence surface area) is derived from the star’s luminosity and effective temperature via the Stefan-Boltzmann law, L = 4 pi R^2 sigma T^4. log(g) is then computed as log10(G M / R^2), with the stellar mass (given in Msun) converted to grams via utils::Msun and G taken from utils::G, so that g comes out in cgs units.

Parameters:

props – Stellar properties, as produced by evaluating the Interpolator1D returned by Tracks2D::getIsochrone at this star’s mass

Returns:

A pair containing the star’s surface area, in cm^2, and log10(g), with g in cgs units (cm/s^2)