Simulating Color-Magnitude Diagrams

generate_stars_mass(mini_vec::AbstractVector{<:Number},
                    mags, mag_names::AbstractVector{String},
                    limit::Number,
                    imf::Distributions.Sampleable{Distributions.Univariate, Distributions.Continuous};
                    dist_mod::Number=0,
                    rng::Random.AbstractRNG=Random.default_rng(),
                    mag_lim::Number = Inf,
                    mag_lim_name::String = "V",
                    binary_model::StarFormationHistories.AbstractBinaryModel =
                        StarFormationHistories.RandomBinaryPairs(0.3))

Generates a random sample of stars from an isochrone defined by the provided initial stellar masses mini_vec, absolute magnitudes mags, and filter names mag_names with total population birth stellar mass limit (e.g., 1e5 solar masses). Initial stellar masses are sampled from the provided imf.

Arguments

• mini_vec::AbstractVector{<:Number} contains the initial masses (in solar masses) for the stars in the isochrone. The isochrone points must be in evolutionary order (e.g., in order of equivalent evolutionary points), as given by isochrone libraries; see StarFormationHistories.ingest_isochrone.
• mags contains the absolute magnitudes from the isochrone in the desired filters corresponding to the same stars as provided in mini_vec. mags is internally interpreted and converted into a standard format by StarFormationHistories.ingest_mags. Valid inputs are:
  • mags::AbstractVector{AbstractVector{<:Number}}, in which case the length of the outer vector length(mags) can either be equal to length(mini_vec), in which case the length of the inner vectors must all be equal to the number of filters you are providing, or the length of the outer vector can be equal to the number of filters you are providing, and the length of the inner vectors must all be equal to length(mini_vec); this is the more common use-case.
  • mags::AbstractMatrix{<:Number}, in which case mags must be 2-dimensional. Valid shapes are size(mags) == (length(mini_vec), nfilters) or size(mags) == (nfilters, length(mini_vec)), with nfilters being the number of filters you are providing.
• mag_names::AbstractVector{String} contains strings describing the filters you are providing in mags; an example might be ["B","V"]. These are used when mag_lim is finite to determine what filter you want to use to limit the faintest stars you want returned.
• limit::Number gives the total birth stellar mass of the population you want to sample. See the "Notes" section on population masses for more information.
• imf::Distributions.Sampleable{Distributions.Univariate, Distributions.Continuous} is a sampleable continuous univariate distribution implementing a stellar initial mass function with a defined rand(rng::Random.AbstractRNG, imf) method to use for sampling masses. All instances of Distributions.ContinuousUnivariateDistribution are also valid. Implementations of commonly used IMFs are available in InitialMassFunctions.jl.

Keyword Arguments

• dist_mod::Number=0 is the distance modulus (see StarFormationHistories.distance_modulus) you wish to have added to the returned magnitudes to simulate a population at a particular distance.
• rng::Random.AbstractRNG=Random.default_rng() is the rng instance that will be used to sample the stellar initial masses from imf.
• mag_lim::Number=Inf gives the faintest apparent magnitude for stars you want to be returned in the output. Stars fainter than this magnitude will still be sampled and contribute properly to the total mass of the population, but they will not be returned.
• mag_lim_name::String="V" gives the filter name (as contained in mag_names) to use when considering if a star is fainter than mag_lim. This is unused if mag_lim is infinite.
• binary_model::StarFormationHistories.AbstractBinaryModel=StarFormationHistories.RandomBinaryPairs(0.3) is an instance of a model for treating binaries; currently provided options are NoBinaries, RandomBinaryPairs, and BinaryMassRatio.

Returns

(sampled_masses, sampled_mags) defined as

• sampled_masses::Vector{SVector{N,eltype(imf)}} is a vector containing the initial stellar masses of the stars sampled by sample_system; see that method's documentation for details on format. In short, each StaticArrays.SVector represents one stellar system and each entry in each StaticArrays.SVector is one star in that system. Entries will be 0 when companions could have been sampled but were not (i.e., when using a binary_model that supports multi-star systems) and for stars that have died (i.e., with initial masses above maximum(mini_vec)), so that summing the entries gives the birth mass of the returned stars that are still alive. A system is returned only when at least one of its stars has a mass within extrema(mini_vec) and the combined system is brighter than mag_lim. Stars with initial masses below minimum(mini_vec) keep their masses if they are returned at all (e.g., if they are a secondary of a living primary), but contribute no light. Stars keep their positions in each system, so the entry for a primary that has died is 0 while its companion's is not.
• sampled_mags::Vector{SVector{N,<:Number}} is a vector containing StaticArrays.SVectors with the multi-band magnitudes of the sampled stars. To get the magnitude of star i in band j, you index as sampled_mags[i][j]. This can be reinterpreted as a 2-dimensional Matrix with reduce(hcat,sampled_mags).

Notes

Population Masses

Given a particular isochrone with an initial mass vector mini_vec, it will never cover the full range of stellar birth masses because stars that die before present-day are not included in the isochrone. However, these stars were born, and so contribute to the total birth mass of the system. There are two ways to properly account for this lost mass when sampling:

1. Set the upper limit for masses that can be sampled from the imf distribution to a physical value for the maximum birth mass of stars (e.g., 100 solar masses). In this case, these stars will be sampled from imf, and will contribute their masses to the system, but they will not be returned if their birth mass is greater than maximum(mini_vec). This is typically easiest for the user and only results in ∼15% loss of efficiency for 10 Gyr isochrones. This approach is preferred when sampling with binaries.
2. Set the upper limit for masses that can be sampled from the imf distribution to maximum(mini_vec) and adjust limit to respect the amount of initial stellar mass lost by not sampling higher mass stars. This can be calculated as new_limit = limit * ( QuadGK.quadgk(x->x*pdf(imf,x), minimum(mini_vec), maximum(mini_vec))[1] / QuadGK.quadgk(x->x*pdf(imf,x), minimum(imf), maximum(imf))[1] ), with the multiplicative factor being the fraction of the population stellar mass contained in stars with initial masses between minimum(mini_vec) and maximum(mini_vec); this factor is the ratio

\[\frac{\int_a^b \ m \times \frac{dN(m)}{dm} \ dm}{\int_0^∞ \ m \times \frac{dN(m)}{dm} \ dm}.\]

Note that, if binaries are included, this approach only forms binary pairs between stars whose masses are less than maximum(mini_vec). This is probably not desired, so we recommend the previous approach when including binaries.

(sampled_masses, sampled_mags) =  generate_stars_mag(mini_vec::AbstractVector{<:Number}, mags, mag_names::AbstractVector{String}, absmag::Real, absmag_name::String, imf::Distributions.Sampleable{Distributions.Univariate,Distributions.Continuous}; dist_mod::Number=0, rng::AbstractRNG=default_rng(), mag_lim::Number=Inf, mag_lim_name::String="V", binary_model::StarFormationHistories.AbstractBinaryModel=RandomBinaryPairs(0.3))

Generates a mock stellar population from an isochrone defined by the provided initial stellar masses mini_vec, absolute magnitudes mags, and filter names mag_names. The population is sampled to a total absolute magnitude absmag::Real (e.g., -7 or -12) in the filter absmag_name::String (e.g., "V" or "F606W") which is contained in the provided mag_names::AbstractVector{String}. Other arguments are shared with generate_stars_mass, which contains the main documentation.

Notes

Population Magnitudes

Unlike when sampling a population to a fixed initial birth mass, as is implemented in generate_stars_mass, when generating a population up to a fixed absolute magnitude, only stars that survive to present-day contribute to the flux of the population. If you choose to limit the apparent magnitude of stars returned by passing the mag_lim and mag_lim_name keyword arguments, stars fainter than your chosen limit will still be sampled and will still contribute their luminosity to the total population, but they will not be contained in the returned output.

(sampled_masses, sampled_mags) = generate_stars_mass_composite(mini_vec::AbstractVector{<:AbstractVector{<:Number}}, mags::AbstractVector, mag_names::AbstractVector{String}, limit::Number, massfrac::AbstractVector{<:Number}, imf::Sampleable{Univariate,Continuous}; rng::Random.AbstractRNG=Random.default_rng(), kws...)

Generates a random sample of stars with a complex star formation history using multiple isochrones. Very similar to generate_stars_mass except the isochrone-related arguments mini_vec and mags should now be vectors of vectors containing the relevant data for the full set of isochrones to be considered. The total birth stellar mass of the sampled population is given by limit. The proportion of this mass allotted to each of the individual isochrones is given by the entries of the massfrac vector. This basically just proportions limit according to massfrac and calls generate_stars_mass for each of the individual stellar populations; as such it is set up to multi-thread across the multiple stellar populations.

Arguments

• mini_vec::AbstractVector{<:AbstractVector{<:Number}} contains the initial masses (in solar masses) for the stars in each isochrone. The isochrone points must be in evolutionary order (e.g., in order of equivalent evolutionary points), as given by isochrone libraries; see StarFormationHistories.ingest_isochrone. The length of mini_vec should be equal to the number of isochrones.
• mags contains the absolute magnitudes from the isochrones in the desired filters corresponding to the same stars as provided in mini_vec. The length of mags should be equal to the number of isochrones. The individual elements of mags are each internally interpreted and converted into a standard format by StarFormationHistories.ingest_mags. The valid formats for the individual elements of mags are:
  • AbstractVector{AbstractVector{<:Number}}, in which case the length of the vector length(mags[i]) can either be equal to length(mini_vec[i]), in which case the length of the inner vectors must all be equal to the number of filters you are providing, or the length of the outer vector can be equal to the number of filters you are providing, and the length of the inner vectors must all be equal to length(mini_vec[i]); this is the more common use-case.
  • AbstractMatrix{<:Number}, in which case mags[i] must be 2-dimensional. Valid shapes are size(mags[i]) == (length(mini_vec[i]), nfilters) or size(mags[i]) == (nfilters, length(mini_vec[i])), with nfilters being the number of filters you are providing.
• mag_names::AbstractVector{String} contains strings describing the filters you are providing in mags; an example might be ["B","V"]. These are used when mag_lim is finite to determine what filter you want to use to limit the faintest stars you want returned. These are assumed to be the same for all isochrones.
• limit::Number gives the total birth stellar mass of the population you want to sample.
• massfrac::AbstractVector{<:Number} is a vector giving the relative fraction of mass allotted to each individual stellar population; length must be equal to the length of mini_vec and mags.
• imf::Distributions.Sampleable{Distributions.Univariate, Distributions.Continuous} is a sampleable continuous univariate distribution implementing a stellar initial mass function with a defined rand(rng::Random.AbstractRNG, imf) method to use for sampling masses. All instances of Distributions.ContinuousUnivariateDistribution are also valid. Implementations of commonly used IMFs are available in InitialMassFunctions.jl.

Keyword Arguments

• rng::Random.AbstractRNG is the random number generator. Because the populations are sampled in parallel, rng is not used to sample stars directly; instead, one seed per population is drawn from rng, and each population is sampled with a copy of rng seeded with Random.seed!(copy(rng), seed). Results for a given rng state are therefore reproducible regardless of the number of threads. rng must support copy and seeding from an integer, which is true of Random.default_rng(), Random.Xoshiro, Random.MersenneTwister, and StableRNGs.StableRNG, among others.

All other keyword arguments kws... are passed to generate_stars_mass; you should refer to that method's documentation for more information.

Returns

• sampled_masses::Vector{Vector{SVector{N,eltype(imf)}}} is a vector of vectors containing the initial stellar masses of the sampled stars. The outer vectors are separated by the isochrone the stars were generated from; i.e., all of sampled_masses[1] were sampled from mini_vec[1] and so on. These can be concatenated into a single vector with reduce(vcat,sampled_masses). The format of the contained StaticArrays.SVectors are as output from sample_system, except that stars that have died are recorded with zero mass as in generate_stars_mass; see those methods' documentation for more details.
• sampled_mags::Vector{Vector{SVector{N,<:Number}}} is a vector of vectors containing StaticArrays.SVectors with the multi-band magnitudes of the sampled stars. The outer vectors are separated by the isochrone the stars were generated from; i.e. all of sampled_mags[1] were sampled from mags[1] and so on. To get the magnitude of star i in band j sampled from isochrone k, you would do sampled_mags[k][i][j]. This can be concatenated into a Vector{SVector} with reduce(vcat,sampled_mags) and a 2-D Matrix with reduce(hcat,reduce(vcat,sampled_mags)).

(sampled_masses, sampled_mags) = generate_stars_mag_composite(mini_vec::AbstractVector{<:AbstractVector{<:Number}}, mags::AbstractVector, mag_names::AbstractVector{String}, absmag::Number, absmag_name::String, fracs::AbstractVector{<:Number}, imf::Sampleable{Univariate,Continuous}; frac_type::Symbol=:lum, rng::Random.AbstractRNG=Random.default_rng(), kws...)

Generates a random sample of stars with a complex star formation history using multiple isochrones. Very similar to generate_stars_mag except the isochrone-related arguments mini_vec and mags should now be vectors of vectors containing the relevant data for the full set of isochrones to be considered. The total absolute magnitude of the sampled population is given by absmag. The proportion of the luminosity allotted to each of the individual isochrones is given by the entries of the frac vector. This basically just proportions the luminosity according to frac and calls generate_stars_mag for each of the individual stellar populations; as such it is set up to multi-thread across the multiple stellar populations.

Arguments

• mini_vec::AbstractVector{<:AbstractVector{<:Number}} contains the initial masses (in solar masses) for the stars in each isochrone. The isochrone points must be in evolutionary order (e.g., in order of equivalent evolutionary points), as given by isochrone libraries; see StarFormationHistories.ingest_isochrone. The length of mini_vec should be equal to the number of isochrones.
• mags contains the absolute magnitudes from the isochrones in the desired filters corresponding to the same stars as provided in mini_vec. The length of mags should be equal to the number of isochrones. The individual elements of mags are each internally interpreted and converted into a standard format by StarFormationHistories.ingest_mags. The valid formats for the individual elements of mags are:
  • AbstractVector{AbstractVector{<:Number}}, in which case the length of the vector length(mags[i]) can either be equal to length(mini_vec[i]), in which case the length of the inner vectors must all be equal to the number of filters you are providing, or the length of the outer vector can be equal to the number of filters you are providing, and the length of the inner vectors must all be equal to length(mini_vec[i]); this is the more common use-case.
  • AbstractMatrix{<:Number}, in which case mags[i] must be 2-dimensional. Valid shapes are size(mags[i]) == (length(mini_vec[i]), nfilters) or size(mags[i]) == (nfilters, length(mini_vec[i])), with nfilters being the number of filters you are providing.
• mag_names::AbstractVector{String} contains strings describing the filters you are providing in mags; an example might be ["B","V"]. These are used when mag_lim is finite to determine what filter you want to use to limit the faintest stars you want returned. These are assumed to be the same for all isochrones.
• absmag::Number gives the total absolute magnitude of the complex population to be sampled.
• absmag_name::String is the name of the filter for which the desired absolute magnitude is absmag; must be contained in mag_names.
• fracs::AbstractVector{<:Number} is a vector giving the relative fraction of luminosity or mass (determined by the frac_type keyword argument) allotted to each individual stellar population; length must be equal to the length of mini_vec and mags.
• imf::Distributions.Sampleable{Distributions.Univariate, Distributions.Continuous} is a sampleable continuous univariate distribution implementing a stellar initial mass function with a defined rand(rng::Random.AbstractRNG, imf) method to use for sampling masses. All instances of Distributions.ContinuousUnivariateDistribution are also valid. Implementations of commonly used IMFs are available in InitialMassFunctions.jl.

Keyword Arguments

• frac_type::Symbol either :lum, in which case fracs is assumed to contain the relative luminosity fractions for each individual isochrone, or :mass, in which case it is assumed that fracs contains birth stellar mass fractions. Mass fractions are converted to luminosity fractions in the absmag_name filter using StarFormationHistories.luminosity_per_mass.
• rng::Random.AbstractRNG is the random number generator. As in generate_stars_mass_composite, one seed per population is drawn from rng and each population is sampled with a seeded copy of rng, so results for a given rng state are reproducible regardless of the number of threads. rng must support copy and seeding from an integer.

All other keyword arguments kws... are passed to generate_stars_mag; you should refer to that method's documentation for more information.

Returns

• sampled_masses::Vector{Vector{SVector{N,eltype(imf)}}} is a vector of vectors containing the initial stellar masses of the sampled stars. The outer vectors are separated by the isochrone the stars were generated from; i.e., all of sampled_masses[1] were sampled from mini_vec[1] and so on. These can be concatenated into a single vector with reduce(vcat,sampled_masses). The format of the contained StaticArrays.SVectors are as output from sample_system, except that stars that have died are recorded with zero mass as in generate_stars_mass; see those methods' documentation for more details.
• sampled_mags::Vector{Vector{SVector{N,<:Number}}} is a vector of vectors containing StaticArrays.SVectors with the multi-band magnitudes of the sampled stars. The outer vectors are separated by the isochrone the stars were generated from; i.e. all of sampled_mags[1] were sampled from mags[1] and so on. To get the magnitude of star i in band j sampled from isochrone k, you would do sampled_mags[k][i][j]. This can be concatenated into a Vector{SVector} with reduce(vcat,sampled_mags) and a 2-D Matrix with reduce(hcat,reduce(vcat,sampled_mags)).

result::Vector{NamedTuple} = add_metadata(starcat, mag_symbols; kws...)

Processes output star catalogs from the composite generate_stars... methods that create complex stellar populations into Vector{NamedTuple} with added metadata which can then be used to construct tables (e.g., TypedTables.Table or DataFrames.DataFrame).

Arguments

• starcat: output from generate_stars_mass_composite or generate_stars_mag_composite.
• mag_symbols: the list of magnitude names to use when unpacking the magnitudes in starcat.

Keyword Arguments

All keyword arguments must have length equal to length(starcat[1]) and length(starcat[2]); that is, equal to the number of SSPs used when constructing the mock catalog. This metadata will then be added to the output result. For example, if the log10(age [yr]) of the SSPs used are in a vector ages and the metallicities in a vector MH, you can call

julia> starcat = generate_stars_mass_composite(...)

julia> filters = (:F475W, :F814W) # Names of filters to associate with the magnitudes in `starcat[2]`

julia> add_metadata(starcat, filters, logAge = ages, MH = MH)
5866911-element Vector{@NamedTuple{mini::Float64, F475W::Float64, F814W::Float64, logAge::Float64, MH::Float64}}:
 (mini = 0.25293134503128467, F475W = 33.53553000158083, F814W = 31.415817763596586, logAge = 6.6, MH = -2.4)
[...]

Initial Masses

The initial masses of each sampled system are split into scalar fields. mini is always the initial mass of the primary star. When starcat was sampled with a binary model (e.g., RandomBinaryPairs), the companion initial mass is given in mini2, which is zero for systems without a companion. Stars that have died are recorded with zero initial mass (see generate_stars_mass), so mini is zero for a primary that has died while its companion is alive, and mini + mini2 is the initial mass of the living stars in the system.

luminosity_per_mass(mini_vec::AbstractVector{<:Number}, mags, imf; mean_mass::Number=mean(imf))

Returns the expected present-day luminosity per unit birth stellar mass of the stellar population defined by the isochrone initial masses mini_vec and absolute magnitudes mags for the provided imf, as a StaticArrays.SVector with one entry per filter. mags accepts the same formats as generate_stars_mass. As there, the isochrone points must be in evolutionary order.

\[\frac{L_i}{M_*} = \frac{1}{\langle m \rangle} \int_{m_\text{min}}^{m_\text{max}} \phi(m) \, 10^{-0.4 \, M_i(m)} \, dm\]

where $\phi(m)$ is the IMF probability density, $M_i(m)$ is the absolute magnitude in filter $i$ along the isochrone at initial mass $m$, $\langle m \rangle$ is mean_mass, and $m_\text{min}$ and $m_\text{max}$ are extrema(mini_vec).

Luminosities are returned as 10^(-0.4 M) for absolute magnitude M, i.e., in units of the luminosity of a source with absolute magnitude zero in the same filter and photometric system (the absolute-magnitude analog of maggies). These are the units returned by StarFormationHistories.mag2flux with its default zeropoint, so a population with total birth stellar mass M_* has expected integrated absolute magnitude flux2mag(M_* * luminosity_per_mass(mini_vec, mags, imf)[i]) in filter i. To convert to solar luminosities, multiply by 10^(0.4 M_⊙), where M_⊙ is the Sun's absolute magnitude in that filter.

As in partial_cmd_smooth, imf can be a Distributions.ContinuousUnivariateDistribution or any callable returning the IMF probability density at a given mass (see StarFormationHistories.dispatch_imf). mean_mass is the mean initial stellar mass of imf; it is computed for you if imf is a Distributions.Distribution, but must be provided otherwise.

The integral is computed with the trapezoidal rule over the isochrone points, so only stars with initial masses in extrema(mini_vec) contribute; stars more massive than the end of the isochrone are assumed to be dead. For binary models that pair stars independently drawn from the IMF (e.g., RandomBinaryPairs), the expected luminosity per unit mass is unchanged by the binary model. For BinaryMassRatio this is an approximation.

mass_from_absmag(mini_vec::AbstractVector{<:AbstractVector{<:Number}}, mags::AbstractVector, mag_names::AbstractVector{String}, absmag::Number, absmag_name::String, massfrac::AbstractVector{<:Number}, imf; mean_mass::Number=mean(imf))

Returns the expected total birth stellar mass of a complex stellar population whose integrated absolute magnitude in the filter absmag_name is absmag. The population consists of the isochrones given by mini_vec and mags with relative birth stellar mass fractions massfrac, which need not be normalized. Arguments are as in generate_stars_mag_composite with frac_type=:mass, and imf and mean_mass are as in StarFormationHistories.luminosity_per_mass, which is used to compute the luminosity per unit birth stellar mass of each isochrone.

\[M_* = \frac{\langle m \rangle \, 10^{-0.4 \, M} \sum_j x_j}{\sum_j x_j \int_{m_{\text{min},j}}^{m_{\text{max},j}} \phi(m) \, 10^{-0.4 \, M_j(m)} \, dm}\]

where $M$ is absmag, $x_j$ is massfrac[j], $M_j(m)$ is the absolute magnitude in filter absmag_name along isochrone $j$ at initial mass $m$, $m_{\text{min},j}$ and $m_{\text{max},j}$ are extrema(mini_vec[j]), and $\phi(m)$ and $\langle m \rangle$ are as in StarFormationHistories.luminosity_per_mass. Equivalently, $M_*$ is the mass for which $M_* \sum_j (x_j / \sum_k x_k) (L/M_*)_j = 10^{-0.4 \, M}$, where $(L/M_*)_j$ is the luminosity per unit birth stellar mass of isochrone $j$.

This is the expected mass; a population sampled to absmag with generate_stars_mag_composite will scatter around it due to stochastic sampling of the IMF.

new_mags [, good_idxs] = model_cmd(mags::AbstractVector{<:AbstractVector{<:Number}}, 
                                   errfuncs, 
                                   completefuncs,
                                   biasfuncs = [zero for i in first(mags)]; 
                                   rng::Random.AbstractRNG = Random.default_rng(), 
                                   ret_idxs::Bool = false)

Simple method for modeling photometric error and incompleteness to "mock observe" a pure catalog of stellar photometry, such as those produced by generate_stars_mass and generate_stars_mag. This method assumes Gaussian photometric errors that are uncorrelated between filters. Each of errfuncs, completefuncs, and biasfuncs may be either an iterable of single-filter callables as described below, or a single joint callable taking one magnitude per filter, as described in the "Joint Photometric Models" section of partial_cmd_smooth; for example, completefuncs(mags[i]...) then gives the detection probability of star i, and errfuncs(mags[i]...) gives its photometric error in each filter.

Arguments

• mags::AbstractVector{<:AbstractVector{<:Number}}: a vector of vectors giving the magnitudes of each star to be modeled. The first index is the per-star index and the second index is the per-filter index (so mags[10][2] would give the magnitude of the tenth star in the second filter). This is the same format as the magnitudes returned by generate_stars_mass and generate_stars_mag; to use output from the composite versions, you must first reduce(vcat,mags) before passing to this function.
• errfuncs: an iterable (typically a vector or tuple) of callables (typically functions or interpolators) with length equal to the number of filters contained in the elements of mags. This iterable must contain callables that, when called with the associated magnitudes from mags, will return the expected 1-σ photometric error at that magnitude. The organization is such that the photometric error for star i in band j is σ_ij = errfuncs[j](mags[i][j]).
• completefuncs: an iterable (typically a vector or tuple) of callables (typically functions or interpolators) with length equal to the number of filters contained in the elements of mags. This iterable must contain callables that, when called with the associated magnitudes from mags, will return the probability that a star with that magnitude in that band will be found in your color-magnitude diagram (this should include the original detection probability and any post-detection quality, morphology, or other cuts). The organization is such that the detection probability for star i in band j is c_ij = completefuncs[j](mags[i][j]).
• biasfuncs = [zero for i in first(mags)]: an iterable (typically a vector or tuple) of callables (typically functions or interpolators) with length equal to the number of filters contained in the elements of mags. This iterable must contain callables that, when called with the associated magnitudes from mags, will return the expectation value of the photometric bias at that magnitude. We define the photometric bias as the typical difference between measured and input magnitudes; i.e., ⟨measured - input⟩. This definition matches that of the bias functions returned by process_ASTs). By default this method assumes no photometric bias.

Keyword Arguments

• rng::AbstractRNG = Random.default_rng(): The object to use for random number generation.
• ret_idxs::Bool: whether to return the indices of the input mags for the stars that were successfully "observed" and are represented in the output new_mags.

Returns

• new_mags: an object similar to mags (i.e., a Vector{Vector{<:Number}}, Vector{SVector{N,<:Number}}, etc.) containing the magnitudes of the mock-observed stars. This will be shorter than the provided mags vector as we are modeling photometric incompleteness, and the magnitudes will also have random photometric errors added to them. This can be reinterpreted as a 2-dimensional Matrix with reduce(hcat,new_mags).
• good_idxs: if ret_idxs is true, the vector of indices into the input mags for the stars that were successfully "observed" and are represented in the output new_mags.

Notes

• This is a simple implementation that seeks to show a simple example of how one can post-process catalogs of "pure" stars from methods like generate_stars_mass and generate_stars_mag to include observational effects. This method assumes Gaussian photometric errors, which may not, in general, be accurate. When completefuncs contains single-filter callables, it also assumes that the total detection probability can be modeled as the product of the single-filter detection probabilities (i.e., that the completeness functions are separable across filters). This can be a reasonable assumption when you have separate photometric catalogs derived for each filter and you only collate them afterwards, but it is generally not a good assumption for detection algorithms that operate simultaneously on multi-band photometry or for catalogs with selection criteria applied in multiple filters – the completeness functions for these cases are generally not separable and are better represented by a joint completefuncs.

surviving_fraction(mmax::Number, imf, binary_model::StarFormationHistories.AbstractBinaryModel=NoBinaries())

Returns the fraction of the stars born in a simple stellar population that are still alive (i.e., that have not yet ended their lives as white dwarfs, neutron stars, or black holes), given the initial mass mmax of the most massive living star (e.g., maximum of the initial masses of an isochrone). Stars with initial masses below mmax are alive. Each star in a binary system is counted separately, so this is not the fraction of stellar systems that are still alive.

imf is a Distributions.ContinuousUnivariateDistribution with a defined cdf, and binary_model is as in generate_stars_mass. With NoBinaries or RandomBinaryPairs, each star is drawn independently from imf, so the result is cdf(imf, mmax). With BinaryMassRatio, imf is the distribution of system masses, and a binary system can lose its primary while its secondary is still alive. A system of mass $M$ and mass ratio $q$ has a primary that is alive for $M \leq m_\text{max} (1 + q)$ and a secondary that is alive for $M \leq m_\text{max} (1 + q) / q$, so for a binary system fraction $b$ (StarFormationHistories.binary_system_fraction) and mass ratio distribution $g(q)$ the result is

\[\frac{1}{1 + b} \left[ (1 - b) \, F(m_\text{max}) + b \int g(q) \left( F \left( m_\text{max} (1 + q) \right) + F \left( \frac{m_\text{max} (1 + q)}{q} \right) \right) dq \right]\]

where $F$ is the cdf of imf and $1 + b$ is the mean number of stars per system.

surviving_fraction(massvec::AbstractVector{<:AbstractVector}, coeffs::AbstractVector, imf, binary_model::StarFormationHistories.AbstractBinaryModel=NoBinaries())

Returns the fraction of the stars born in a complex stellar population that are still alive. massvec contains the initial masses of the living stars in each simple stellar population (e.g., their isochrones, which must not include remnants), and coeffs gives the birth stellar mass of each population (up to a common normalization), which is proportional to the number of stars born as all populations share imf and binary_model. The result is the average of surviving_fraction(maximum(massvec[i]), imf, binary_model) over the populations weighted by coeffs, so only the most massive star of each population is used.

recycling_fraction(args...)

Returns the fraction of the stars that were born that have ended their lives (i.e., have become white dwarfs, neutron stars, or black holes). Arguments are as in either method of surviving_fraction. By definition this is equal to 1 - surviving_fraction(args...).

surviving_mass_fraction(mmax::Number, imf, binary_model::StarFormationHistories.AbstractBinaryModel=NoBinaries(); mean_mass::Number=mean(imf))

Returns the fraction of the birth stellar mass of a simple stellar population that is in stars that are still alive (i.e., that have not yet ended their lives as white dwarfs, neutron stars, or black holes), given the initial mass mmax of the most massive living star (e.g., maximum of the initial masses of an isochrone). Stars with initial masses below mmax are alive. The mass counted is the initial (birth) mass of the living stars, without mass lost during their evolution or the mass of remnants.

imf is a Distributions.ContinuousUnivariateDistribution, whose partial mean is computed with InitialMassFunctions.partial_mean (in closed form for the IMFs of InitialMassFunctions.jl and by numerical integration otherwise), and mean_mass is its mean, as in partial_cmd_smooth. binary_model is as in generate_stars_mass. With NoBinaries or RandomBinaryPairs, each star is drawn independently from imf, so the result is

\[\frac{1}{\langle m \rangle} \int_{m_\text{min}}^{m_\text{max}} m \, \phi(m) \, dm\]

where $\phi(m)$ is the IMF probability density and $\langle m \rangle$ is mean_mass. With BinaryMassRatio, imf is the distribution of system masses, and a binary system can lose its primary while its secondary is still alive, so the fraction of the stellar mass in binary systems (StarFormationHistories.binary_mass_fraction) is integrated over the distribution of mass ratios with the primary and secondary counted separately. This is the expected fraction of the mass in the stars sampled by generate_stars_mass that are still alive, including those not returned because they are below the isochrone.

surviving_mass_fraction(massvec::AbstractVector{<:AbstractVector}, coeffs::AbstractVector, imf, binary_model::StarFormationHistories.AbstractBinaryModel=NoBinaries(); mean_mass::Number=mean(imf))

Returns the fraction of the birth stellar mass of a complex stellar population that is in stars that are still alive. massvec contains the initial masses of the living stars in each simple stellar population (e.g., their isochrones, which must not include remnants), and coeffs gives the birth stellar mass of each population (up to a common normalization). The result is the average of surviving_mass_fraction(maximum(massvec[i]), imf, binary_model; mean_mass) over the populations weighted by coeffs, so only the most massive star of each population is used.

recycling_mass_fraction(args...; kws...)

Returns the fraction of the birth stellar mass that was in stars that have ended their lives (i.e., have become white dwarfs, neutron stars, or black holes). Arguments are as in either method of surviving_mass_fraction. By definition this is equal to 1 - surviving_mass_fraction(args...; kws...).

new_mags = ingest_mags(mini_vec::AbstractVector, mags::AbstractVector{T}) where {S <: Number, T <: AbstractVector{S}}
new_mags = ingest_mags(mini_vec::AbstractVector, mags::AbstractMatrix{S}) where S <: Number

Reinterprets provided mags to be in the correct format for input to Interpolations.interpolate.

Returns

• new_mags::Base.ReinterpretArray{StaticArrays.SVector}: a length(mini_vec) vector of StaticArrays.SVectors containing the same data as mags but formatted for input to Interpolations.interpolate.

(new_mini_vec, new_mags) = ingest_isochrone(mini_vec::AbstractVector, mags::AbstractVector)

Prepares the initial masses mini_vec and magnitudes mags (one element per isochrone point, either a number or a vector of magnitudes in several filters) of an isochrone for interpolation in initial mass. The arguments are not mutated.

Isochrone points should be provided in evolutionary order (e.g., in order of equivalent evolutionary points), as given by isochrone libraries. Initial mass then increases along the isochrone, but interpolated isochrones can have small decreases in initial mass between evolutionary phases (e.g., from the tip of the red giant branch to the horizontal branch). Sorting by initial mass would interleave the points of the two phases, so instead each decrease is removed by shifting the initial masses of all later points up by the size of the decrease. This preserves the differences in initial mass within each phase, which set the relative numbers of stars in each phase. Repeated initial masses are then separated with Interpolations.deduplicate_knots!.

Points in evolutionary order trace the isochrone, while points out of order jump across it. If the path through the magnitudes of consecutive points is more than twice as long in the given order as in order of initial mass, the points are assumed not to be in evolutionary order; a warning is issued and they are sorted by initial mass instead.

(mmin, mmax) = mass_limits(mini_vec::AbstractVector{<:Number}, mags::AbstractVector{T},
                 mag_names::AbstractVector{String}, mag_lim::Number,
                 mag_lim_name::String) where T <: AbstractVector{<:Number}

Calculates initial mass limits that reflect a given faint-end magnitude limit.

Arguments

• mini_vec::AbstractVector{<:Number}: a length nstars vector containing initial stellar masses.
• mags::AbstractVector{<:AbstractVector{<:Number}}: a length nstars vector, with each element being a length nfilters vector giving the magnitudes of each star in the filters mag_names.
• mag_names::AbstractVector{String}: a vector giving the names of each filter as strings.
• mag_lim::Number: the faint-end magnitude limit you wish to use.
• mag_lim_name::String: the name of the filter in which mag_lim is to be applied. Must be contained in mag_names.

Returns

• mmin::eltype(mini_vec): the initial mass corresponding to your requested mag_lim in the filter mag_lim_name. If all stars provided are brighter than your requested mag_lim, then this will be equal to minimum(mini_vec).
• mmax::eltype(mini_vec): the maximum valid mass in mini_vec; simply maximum(mini_vec).

Examples

julia> mass_limits([0.05,0.1,0.2,0.3], [[4.0],[3.0],[2.0],[1.0]], ["F090W"], 2.5, "F090W")
(0.15, 0.3)

julia> mass_limits([0.05,0.1,0.2,0.3], [[4.0,3.0],[3.0,2.0],[2.0,1.0],[1.0,0.0]], ["F090W","F150W"], 2.5, "F090W")
(0.15, 0.3)