Developer documentation
This page contains info for people interested in contributing code to Korg. If you have questions, do not hesitate to ask.
Julia development
If you've never developed a package in Julia, here are some tips.
- Once you have a local copy of Korg, you can make it importable locally by
dev
ing it from the Julia environment in which you would like to run it:Pkg.dev("/path/to/Korg")
. If you runPkg.dev("Korg")
instead, Julia will automatically clone the repo to~/.julia/dev/Korg
. (See the Pkg documentation for details.) - When working on your local copy of Korg, Revise is easiest way to make and test changes without constantly restarting your Julia session.
- To run the test suite locally, start a Julia session in the Korg root directory and run
]activate .
thentest
. This will automatically runtest/runtests.jl
in the test environment. - Documentation is generated with Documentor. Continuous integration on github will ensure that the documentation is generated without errors, but it won't catch all formatting problems. If you wish to generate documentation locally in order to check that everything is as expected, start a Julia session in
Korg/docs
, activate the test environment (]activate .
), and runinstantiate
. This downloads and installs the docs dependences, and will only have to be run once. To generate documentation, runjulia --project make.jl
on the command line (the--project
flag activates the local environment). The generated docs can be served fromdocs/build
.
Code guidelines
- Try to be explicit about units throughout the code, particularly when not using CGS.
- Whenever possible, calculations should be precise up to a factor of $10^{-3}$. When it's easy and inexpensive, they should be precise to $10^{-5}$ or better. There are exceptions to this (e.g. the Voigt function), but ideally they will eventually be vanquished.
- Ensure types are generic enough to support dual numbers and autodifferentiation.
- Limit lines to 100 characters.
- Use ASCII characters in the names of functions that are part of the public API.
- Unless they will never be called elsewhere, provide docstrings describing the inputs, assumptions and outputs of any functions you write.
Continuum absorption
Continuum Absorption Kwargs
All bound-free and free-free absorption functions (other than absorption_coef_bf_TOPBase
) support the following keyword arguments:
error_oobounds::Bool
: specifies the function's behavior when it encounters invalidν
orT
values. Whenfalse
(the default), the function assumes that the absorption coefficient for these invalidν
orT
values is0.0
. Otherwise, aDomainError
is thrown when invalidν
orT
values are encountered.out_α::Union{Nothing,AbstractVector}
: When this isnothing
(the default case), the function will simply allocate a new vector to store the output continuum absorption coefficients. Alternatively, this can be a vector (with the same length as the function'sν
argument). In this case, the function will directly add the computed continuum absorption to the elements of this vector, in-place (the vector is also returned by the function).
Steps for implementing new continuum sources of absorption:
- Define a helper function that computes a single absorption coefficient (in units of cm⁻²). The function should accept
ν
(in Hz) andT
(in K) as the first and second arguments, respectively. The convention is for it should share a name with the corresponding public function, but have an underscore pre-appended (e.g. we define_H_I_bf
to help implementH_I_bf
). - The public function is the function constructed and returned by
Korg.ContinuumAbsorption.bounds_checked_absorption
that wraps the above helper function. This wrapper function implements bounds-checking forν
andT
and supports the keyword arguments described in Continuum Absorption Kwargs. - Add a docstring describing the new function. At the very least, please describe any non-standard arguments and include a reference in the docstring to the source where the function was taken from.
- Add a line to
doc/src/API.md
under theContinuum absorption
heading to render the docstring of your new function. - Add a line to
total_continuum_absorption
that calls the new public function for absorption.
The first two steps may not apply for sources that don't directly depend on ν
and T
(e.g. absorption from scattering).
Manifest.toml
and Project.toml
Read more about these files here. These files are used by the julia package manager. Project.toml
records dependencies, and you'll notice that the test and docs directories have their own Project.toml
s for test and documentation-specific dependencies. Manifest.toml
records exact package versions used when a package was run. It enables someone else to reproduce the exact environment and results later. There are a few directories containing scripts that generate Korg's data files which have there own Projects.toml
s and Manifest.toml
s, for example data/bf_cross-sections/
.
Where to put data
If you are adding data to Korg, data/README
provides an overview of the options and how to decide between them.
Complete API
Korg._alpha_5000_default_linelist
— ConstantThe default linelist for calculating the absorption coefficient at 5000 Å. This for internal use when the provided linelist doesn't cover the region and a radiative transfer scheme using τ_5000 is used.
See also _load_alpha_5000_linelist
.
Korg.default_alpha_elements
— ConstantThe "alpha elements" are defined as O, Ne, Mg, Si, S, Ar, Ca, Ti, i.e. all elements with even atomic numbers from 8-22. This definition is used by default in format_A_X
, get_metals_H
, and get_alpha_H
, but can be overridden with a keyword argument.
Korg.Formula
— TypeRepresents an atom or molecule, irrespective of its charge.
Korg.Line
— MethodLine(wl::F, log_gf::F, species::Species, E_lower::F,
gamma_rad::Union{F, Missing}=missing, gamma_stark::Union{F, Missing}=missing,
vdw::Union{F, Tuple{F, F}, Missing}, missing) where F <: Real
Arguments:
wl
: wavelength (Assumed to be in cm if < 1, otherwise in Å)log_gf
: (log base 10) oscillator strength (unitless)species
: theSpecies
associated with the lineE_lower
: The energy (excitiation potential) of the lower energy level (eV)
Optional Arguments (these override default recipes):
gamma_rad
: Fundemental widthgamma_stark
: per-perturber Stark broadening width at 10,000 K (s⁻¹).vdW
: If this is present, it may may belog10(γ_vdW)
, assumed if negative- 0, corresponding to no vdW broadening
- A fudge factor for the Unsoeld approximation, assumed if between 0 and 20
- The ABO parameters as packed float (assumed if >= 20) or a
Tuple
,(σ, α)
.
This behavior is intended to mirror that of Turbospectrum as closely as possible.
See approximate_gammas
for more information on the default recipes for gamma_stark
and vdW
.
Note the the "gamma" values here are FWHM, not HWHM, of the Lorenztian component of the line profile, and are in units of s⁻¹.
Korg.MolecularCrossSection
— MethodMolecularCrossSection(linelist, wl_params...; cutoff_alpha=1e-30, log_temp_vals=3:0.025:5, verbose=true)
Precompute the molecular absorption cross section for a given linelist and set of wavelengths. The MolecularCrossSection
object can be passed to synthesize
and potentially speed up the calculation significantly. At present, Korg only supports precomputed cross-sections created by this function, though they can be saved and loaded using save_molecular_cross_section
and read_molecular_cross_section
.
Arguments
linelist
: A vector ofLine
objects representing the molecular linelist. These must be of the same species.wl_params...
: Parameters specifying the wavelengths in the same format thatsynthesize
expects.
Keyword Arguments
cutoff_alpha
(default: 1e-30): The value of the single-line absorption coefficient (in cm^-1) at which to truncate the profile.log_temp_vals
(default: 3:0.025:5): The log10 of the temperatures at which to precompute the cross-section.verbose
(default: true): Whether to print progress information.
The default values of vmic_vals
, log_temp_vals
, and cutoffalpha
were chosen to ensure that lines in the APOGEE linelist ([
getAPOGEEDR17linelist`](@ref)) could be accurately reproduced (better than 10^-3 everywhere). You should verify that they yield acceptable accuracy for other applications by comparing spectra synthesize with and without precomputing the molecular cross-section.
Korg.PlanarAtmosphere
— MethodPlanarAtmosphere(atm::ShellAtmosphere)
Construct a planar atmosphere with the data from a shell atmosphere. Mostly useful for testing.
Korg.ShellAtmosphere
— MethodShellAtmosphere(atm::PlanarAtmosphere, R)
Construct a shell atmosphere with the data from a planar atmosphere and an outer radius. Mostly useful for testing.
Korg.Species
— TypeRepresents an atom or molecule (a Formula
) with a particular number of electrons (regardless of their configuration).
Korg.Species
— MethodSpecies(code)
Parse the "species code" in many of the forms in which it is often specified and return an object representing the species. code
can be either a string or a float.
Examples
- "H I" -> H I
- "H 1" -> H I
- "H 1" -> H I
- "H_1" -> H I
- "H.I" -> H I
- "H 2" -> H II
- "H2" -> H₂
- "H" -> H I
- "01.00" → H I
- "02.01" → He II
- "02.1000" → He II
- "0608" → CO I
To parse at compile time, use the species
string macro, i.e. species"H I"
. This is important in hot inner loops.
MOOG codes which include isotopic information will not be parsed correctly by this function, though read_linelist
handles them correctly. # leading 0s are safe to remove
Korg.SynthesisResult
— TypeSynthesisResult
The result of a synthesis. Returned by synthesize
.
Fields
flux
: the output spectrumcntm
: the continuum at each wavelengthintensity
: the intensity at each wavelength and mu value, and possibly each layer in the model atmosphere, depending on the radiative transfer scheme.alpha
: the linear absorption coefficient at each wavelength and atmospheric layer a Matrix of size (layers x wavelengths)mu_grid
: a vector of tuples containing the μ values and weights used in the radiative transfer calculation. Can be controlled with themu_values
keyword argument.number_densities
: A dictionary mappingSpecies
to vectors of number densities at each atmospheric layerelectron_number_density
: the electron number density at each atmospheric layerwavelengths
: The vector of vacuum wavelengths (in Å) over which the synthesis was performed. Ifair_wavelengths=true
this will not be the same as the input wavelengths.subspectra
: A vector of ranges which can be used to index intoflux
to extract the spectrum for each range provided inwavelength_ranges
. If you use the standardλ_start
,λ_stop
,λ_step
arguments, this will be a vector containing only one range.
Korg.Wavelengths
— TypeKorg.Wavelengths(wl_params...; air_wavelengths=false, wavelength_conversion_warn_threshold=1e-4)
Construct a Wavelengths
object which represents the (possibly non-contiguous) wavelengths for which to compute a spectrum. The wavelengths can be specified with an upper and lower bound, or a vector of upper and lower bounds. For example,
Korg.Wavelengths(5000, 5500)
Korg.Wavelengths([(5000, 5500), (6000, 6500)])
Keyword Arguments
air_wavelengths
(default:false
): Whether or not the input wavelengths are air wavelengths to be converted to vacuum wavelengths by Korg. The conversion will not be exact, so that the wavelength range can internally be represented by an evenly-spaced range. If the approximation error is greater thanwavelength_conversion_warn_threshold
, an error will be thrown. (To do wavelength conversions yourself, seeair_to_vacuum
andvacuum_to_air
.)wavelength_conversion_warn_threshold
(default: 1e-4): seeair_wavelengths
. (In Å.)
Korg.Qfactor
— MethodQfactor(synth_flux, synth_wl, obs_wl, LSF_mat; obs_mask=nothing)
Compute the Q factor from the high-resolution theoretical spectrum, the high resolution wavelength grid, the low resolution (observed) wavelength grid, and the LSF matrix. Based on work from Bouchy et al. 2001, A&A, 374, 733. Note that the Q factor is an approximation when the flux uncertainty is not photon-dominated.
Arguments
synth_flux
: High resolution theoretical spectrumsynth_wl
: High resolution wavelength gridobs_wl
: Low resolution wavelength gridLSF_mat
: LSF matrix (seecompute_LSF_matrix
)
Keyword Arguments
obs_mask::Vector{Bool}=nothing
: Mask for the low resolution spectrum to account for pixels masked from the observation
See also: RV_prec_from_Q
and RV_prec_from_noise
Korg.RV_prec_from_Q
— MethodRV_prec_from_Q(Q, RMS_SNR, Npix)
Compute the RV precision (in m/s) from the Q factor, the RMS SNR per pixel, and the number of pixels in the spectrum.
Arguments
Q::Real
: Q factor of the spectrumRMS_SNR::Real
: The root mean squared per pixel SNR of the spectrum. For spectra which are not line-blanketed, this is approximately equal to the average SNR. approximately equal to the SNR.Npix::Real
: Number of pixels in the spectrum
See also: Qfactor
and RV_prec_from_noise
Korg.RV_prec_from_noise
— MethodRV_prec_from_noise(synth_flux, synth_wl, obs_wl, LSF_mat, obs_err; obs_mask=nothing)
Compute the best achievable RV precision given the a spectrum with uncertainties.
Arguments
synth_flux
: High-resolution theoretical spectrumsynth_wl
: High-resolution wavelength gridobs_wl
: Low-resolution wavelength gridLSF_mat
: LSF matrix (seecompute_LSF_matrix
)obs_err
: Noise in the continuum normalized spectrum
Keyword Arguments
obs_mask=nothing
: Mask for the low resolution spectrum to account for pixels masked from the observation
See also: RV_prec_from_Q
and Qfactor
Korg._get_multi_X_H
— MethodGiven a vector of abundances, A_X
, get [I+J+K/H], where Zs = [I,J,K]
is a vector of atomic numbers. This is used to calculate, for example, [α/H] and [metals/H].
Korg._load_alpha_5000_linelist
— Function_load_alpha_5000_linelist([path])
Load the default linelist for calculating the absorption coefficient at 5000 Å. This for internal use when the provided linelist doesn't cover the region and a radiative transfer scheme using τ_5000 is used.
This linelist loaded into Korg._alpha_5000_default_linelist
when Korg is imported.
Korg.air_to_vacuum
— Methodair_to_vacuum(λ; cgs=λ<1)
Convert λ from an air to vacuum. λ is assumed to be in Å if it is ⩾ 1, in cm otherwise. Formula from Birch and Downs (1994) via the VALD website.
See also: vacuum_to_air
.
Korg.all_atomic_species
— Methodall_atomic_species()
Returns an iterator that runs over all atomic species supported by Korg.
Korg.apply_LSF
— Methodapply_LSF(flux, wls, R; window_size=4)
Applies a gaussian line spread function the the spectrum with flux vector flux
and wavelengths wls
in any format accepted by synthesize, e.g. as a pair (λstart, λstop)
) with constant spectral resolution (, $R = \lambda/\Delta\lambda$, where $\Delta\lambda$ is the LSF FWHM. The window_size
argument specifies how far out to extend the convolution kernel in standard deviations.
For the best match to data, your wavelength range should extend a couple $\Delta\lambda$ outside the region you are going to compare.
If you are convolving many spectra defined on the same wavelenths to observational resolution, you will get much better performance using compute_LSF_matrix
.
Keyword Arguments
window_size
(default: 4): how far out to extend the convolution kernel in units of the LSF width (σ, not HWHM)renormalize_edge
(default:true
): whether or not to renormalize the LSF at the edge of the wl range. This doen't matter as long assynth_wls
extends to large and small enough wavelengths.
This is a naive, slow implementation. Do not use it when performance matters.
apply_LSF
will have weird behavior if your wavelength grid is not locally linearly-spaced. It is intended to be run on a fine wavelength grid, then downsampled to the observational (or otherwise desired) grid.
Korg.apply_rotation
— Functionapply_rotation(flux, wls, vsini, ε=0.6)
Given a spectrum flux
sampled at wavelengths wls
for a non-rotating star, compute the spectrum that would emerge given projected rotational velocity vsini
(in km/s) and linear limb-darkening coefficient ε
: $I(\mu) = I(1) (1 - \varepsilon + arepsilon \mu))$. See, for example, Gray equation 18.14.
Korg.approximate_gammas
— Methodapproximate_gammas(wl, species, E_lower; ionization_energies=Korg.ionization_energies)
A simplified form of the Unsoeld (1955) approximation for van der Waals broadening and the Cowley 1971 approximation for Stark broadening, evaluated at 10,000 K. Used for atomic lines with no vdW and stark broadening info in the linelist.
Returns (γ_stark
, log10(γ_vdW))
in Hz, where these are the per-perturber quantities. For autoionizing lines (those for which Eupper > χ), Returns 0.0 for γvdW. Note the the "gamma" values here are FWHM, not HWHM, of the Lorenztian component of the line profile.
In the calculation of n*²
, uses the approximation that $\overbar{r^2} = 5/2 {n^*}^4 / Z^2$ which neglects the dependence on the angular momentum quantum number, l, in the the form given by Warner 1967 (the earliest english work reporting the Unsoeld result).
Korg.approximate_line_strength
— Methodapproximate_line_strength(line::Line, T)
Approximate the line strength (log10(gfλ) - θχ
, in arbitrary units) of a line at temperature T
(K). This can be used to very quickly filter a linelist, and is used to filter very large molecular linelists from ExoMol (see load_ExoMol_linelist
).
Korg.approximate_radiative_gamma
— Methodapproximate_radiative_gamma(wl, log_gf)
Approximate radiate broadening parameter. When using this, make sure that log_gf
is the true value (not adjusted for isotopic abundance).
Korg.autodiffable_conv
— Methodautodiffable_conv(f, g)
Compute the convolution of two vectors, f
and g
. This is a wrapped version of DSP.conv
that handles ForwardDiff.Dual
types via the chain rule.
Korg.blackbody
— Methodblackbody(T, λ)
The value of the Planck blackbody function for temperature T
at wavelength λ
[cm].
Korg.bracket_line_interpolator
— Functionbracket_line_interpolator(m, λ₀, T, nₑ, ξ, λmin, λmax; kwargs...)
This routine numerically convolves the two components of the Brackett line stark profile (quasistatic/Holtsmark and impact) and the doppler profile, if necessary. It returns a pair containing the interpolator and the distance from the line center at which it is defined.
Arguments
m
: the principle quantum number of the upper levelλ₀
: the line center in ÅT
: the temperature [K]nₑ
: the electron number density [cm^-3]ξ
: the microturbulence [cm/s]λmin
: the minimum wavelength at which the profile should be computed (used to avoid calculations outside the required region)λmax
: the mxinimum wavelength at which the profile should be computed
Keyword Arguments
n_wavelength_points
(default=201): the number of wavelengths at which to sample the profiles quasistatic profiles to be convolved.window_size
(default=5): the size of the wavelength range over which the profiles should be calculated, in units of the characteristic profile width
Korg.brackett_line_stark_profiles
— Methodbrackett_line_stark_profiles(m, λs, λ₀, T, nₑ)
Stark-broadened line profile (specialized to Brackett series). Translated and heavily adapted from HLINOP.f by Barklem, who adapted it from Peterson and Kurucz. Mostly follows Griem 1960, and Griem 1967. Ions and distant electrons have E fields which can be treated quasi-statically, leading to a Holtsmark broadening profile.
Returns a pair of vectors containing the impact and quasistatic profiles, respectively.
Arguents:
m
: the upper level of the transitionλs
: the wavelengths at which to calculate the profile [cm]λ₀
: the line center [cm]T
: temperature [K]nₑ
: electron number density [cm^-3]
Korg.brackett_oscillator_strength
— Methodbrackett_oscillator_strength(n, m)
The oscillator strength of the transition from the 4th to the mth energy level of hydrogen. Adapted from HLINOP.f by Peterson and Kurucz.
Comparison to the values in Goldwire 1968 indicates that this is accurate to 10^-4 for the Brackett series.
Korg.chemical_equilibrium
— Methodchemical_equilibrium(T, nₜ, nₑ, absolute_abundances, ionization_energies,
partition_fns, log_equilibrium_constants; x0=nothing)
Iteratively solve for the number density of each species. Returns a pair containing the electron number density and Dict
mapping species to number densities.
arguments:
- the temperature,
T
, in K - the total number density
nₜ
- the electron number density
nₑ
- A Dict of
absolute_abundances
, NX/Ntotal - a Dict of ionization energies,
ionization_energies
. The keys of act as a list of all atoms. - a Dict of partition functions,
partition_fns
- a Dict of log molecular equilibrium constants,
log_equilibrium_constants
, in partial pressure form. The keys ofequilibrium_constants
act as a list of all molecules.
keyword arguments:
x0
(default:nothing
) is an initial guess for the solution (in the format internal tochemical_equilibrium
). If not supplied, a good guess is computed by neglecting molecules.electron_number_density_warn_threshold
(default:0.1
) is the fractional difference between the calculated electron number density and the model atmosphere electron number density at which a warning is issued.electron_number_density_warn_min_value
(default:1e-4
) is the minimum value of the electron number density at which a warning is issued. This is to avoid warnings when the electron number density is very small.
The system of equations is specified with the number densities of the neutral atoms as free parameters. Each equation specifies the conservation of a particular species, e.g. (simplified)
n(O) = n(CO) + n(OH) + n(O I) + n(O II) + n(O III).
In this equation:
n(O)
, the number density of oxygen atoms in any form comesabsolute_abundances
and the total number density (supplied later)n(O I)
is a free parameter. The numerical solver is varying this to satisfy the system of equations.n(O II)
, andn(O III)
come from the Saha (ionization) equation givenn(O I)
n(CO)
andn(OH)
come from the molecular equilibrium constants K, which are precomputed over a range of temperatures.
Equilibrium constants are defined in terms of partial pressures, so e.g.
K(OH) == (p(O) p(H)) / p(OH) == (n(O) n(H)) / n(OH)) kT
Korg.compute_LSF_matrix
— Methodcompute_LSF_matrix(synth_wls, obs_wls, R; kwargs...)
Construct a sparse matrix, which when multiplied with a flux vector defined over wavelenths synth_wls
, applies a gaussian line spead function (LSF) and resamples to the wavelenths obswls
.
Arguments
synth_wls
: the synthesis wavelengths in any form recognized bysynthesize
, e.g. a pair containing a lower and upper bound in Å, a vector of pairs, or a vector of Julia range objects.obs_wls
: the wavelengths of the observed spectrumR
: the resolving power, $R = \lambda/\Delta\lambda$
Keyword Arguments
window_size
(default: 4): how far out to extend the convolution kernel in units of the LSF width (σ, not HWHM)verbose
(default:true
): whether or not to emit warnings and information to stdout/stderr.step_tolerance
: the maximum difference between adjacent wavelengths insynth_wls
for them to be considered linearly spaced. This is only used ifsynth_wls
is a vector of wavelengths rather than a range or vector or ranges.
For the best match to data, your wavelength range should extend a couple $\Delta\lambda$ outside the region you are going to compare.
Korg.apply_LSF
can apply an LSF to a single flux vector efficiently. This function is relatively slow, but one the LSF matrix is constructed, convolving spectra to observational resolution via matrix multiplication is fast.
Korg.contained
— Methodcontained(value, interval)
Returns whether value
is contained by interval
.
Examples
julia> contained(0.5, Interval(1.0, 10.0))
false
julia> contained(5.0, Interval(1.0, 10.0))
true
Korg.contained_slice
— Methodcontained_slice(vals, interval)
Returns a range of indices denoting the elements of vals
(which are assumed to be sorted in increasing order) that are contained by interval
. When no entries are contained by interval, this returns (1,0)
(which is a valid empty slice).
Korg.doppler_width
— Methoddoppler_width(λ₀ T, m, ξ)
The standard deviation of of the doppler-broadening profile. In standard spectroscopy texts, the Doppler width often refers to σ√2, but this is σ
Korg.eachfreq
— Methodeachfreq(wls::Wavelengths)
Returns an array of the frequencies corresponding to the wavelengths in wls
(in Hz). They are sorted, i.e. in reverse order of the wavelengths.
Korg.eachwindow
— Methodeachwindows(wls::Wavelengths)
Returns an iterator over the wavelength ranges (λ_low, λ_hi)
in wls
(in cm, not Å).
Korg.exponential_integral_1
— Methodexponential_integral_1(x)
Compute the first exponential integral, E1(x). This is a rough approximation lifted from Kurucz's VCSE1F. Used in brackett_line_profile
.
Korg.filter_linelist
— Methodfilter_linelist(linelist, wls, line_buffer)
Return a new linelist containing only lines within the provided wavelength ranges.
Korg.format_A_X
— Methodformat_A_X(default_metals_H, default_alpha_H, abundances; kwargs... )
Returns a 92 element vector containing abundances in $A(X)$ ($\log_{10}(X/H) + 12$) format for elements from hydrogen to uranium.
Arguments
You can specify abundance with these positional arguments. All are optional, but if default_alpha_H
is provided, default_metals_H
must be as well.
default_metals_H
(default: 0), i.e. [metals/H] is the $\log_{10}$ solar-relative abundance of elements heavier than He. It is overridden bydefault_alpha
andabundances
on a per-element basis.default_alpha_H
(default: same asdefault_metals_H
), i.e. [alpha/H] is the $\log_{10}$ solar-relative abundance of the alpha elements (Seealpha_elements
, below). It is overridden byabundances
on a per-element basis.abundances
is aDict
mapping atomic numbers or symbols to [$X$/H] abundances. (Setsolar_relative=false
to use $A(X)$ abundances instead.) These overridedefault_metals_H
. This is the only way to specify an abundance of He that is non-solar.
Keyword arguments
solar_relative
(default: true): When true, interpret abundances as being in [$X$/H] ($\log_{10}$ solar-relative) format. When false, interpret them as $A(X)$ abundances, i.e. $A(x) = \log_{10}(n_X/n_\mathrm{H}) + 12$, where $n_X$ is the number density of $X$. Note that abundances not specified default to the solar value still depend on the solar value, as they are set according todefault_metals_H
anddefault_alpha_H
.solar_abundances
(default:Korg.asplund_2020_solar_abundances
) is the set of solar abundances to use, as a vector indexed by atomic number.Korg.asplund_2009_solar_abundances
andKorg.grevesse_2007_solar_abundances
are also provided for convenience.alpha_elements
(default:Korg.default_alpha_elements
): vector of atomic numbers of the alpha elements. (Useful since conventions vary.)
Korg.get_APOGEE_DR17_linelist
— Methodget_APOGEE_DR17_linelist(; include_water=true)
The APOGEE DR 17 linelist. It ranges from roughly 15,000 Å to 17,000 Å. It is nearly the same at the DR 16 linelist described in Smith+ 2021.
Korg.get_GALAH_DR3_linelist
— Methodget_GALAH_DR3_linelist()
The GALAH DR 3 linelist (also used for DR 4). It ranges from roughly 4,675 Å to 7,930 Å. This linelist is based on, but distinct from Heiter 2021. See Buder et al. 2021 for details.
Korg.get_GES_linelist
— Methodget_GES_linelist()
The Gaia-ESO survey linelist from Heiter et al. 2021. This linelist contains > 15 million lines, which means that it can take a while to synthesize spectra. If you don't need molecular lines, you can set include_molecules=false
to speed things up.
Keyword Arguments
include_molecules
(default:true
): whether to include molecular lines.
Korg.get_VALD_solar_linelist
— Methodget_VALD_solar_linelist()
Get a VALD "extract stellar" linelist produced at solar parameters, with the "threshold" value set to 0.01. It was downloaded on 2021-05-20. It is intended to be used for quick tests only.
Korg.get_alpha_5000_linelist
— Methodget_alpha_5000_linelist(linelist)
Arguments:
linelist
: the synthesis linelist, which will be used if it covers the 5000 Å region.
Return a linelist which can be used to calculate the absorption at 5000 Å, which is required for the standard radiative transfer scheme. If the provided linelist doesn't contain lines near 5000 Å, use the built-in one. (see _load_alpha_5000_linelist
)
Korg.get_alpha_H
— Methodget_alpha_H(A_X; kwargs...)
Calculate [α/H] given a vector, A_X
of absolute abundances, $A(X) = \log_{10}(n_α/n_\mathrm{H})$. Here, the alpha elements are defined to be O, Ne, Mg, Si, S, Ar, Ca, Ti. See also get_metals_H
.
Keyword Arguments
solar_abundances
(default:Korg.asplund_2020_solar_abundances
) is the set of solar abundances to use, as a vector indexed by atomic number.Korg.asplund_2009_solar_abundances
,Korg.grevesse_2007_solar_abundances
, andKorg.magg_2022_solar_abundances
are also provided for convenience.alpha_elements
(default:Korg.default_alpha_elements
): vector of atomic numbers of the alpha elements. (Useful since conventions vary.)
Korg.get_atoms
— Methodget_atoms(x)
Returns an array view containing the atomic number of each atom that makes up the formula or species x. E.g. get_atoms(Korg.species"H2O")
yields [1, 1, 8].
Korg.get_electron_number_densities
— Methodgetelectronnumberdensities(atm::ModelAtmosphere) = [l.electronnumber_density for l in atm.layers]
This is a convienince functions for making plots, etc. Note that it doesn't access quantities in a memory-efficient order.
Korg.get_gas_pressures
— Methodget_gas_pressures(atm::ModelAtmosphere)
This is a convienince functions for making plots, etc. Note that it doesn't access quantities in a memory-efficient order.
Korg.get_log_nK
— Methodget_log_nK(mol, T, log_equilibrium_constants)
Given a molecule, mol
, a temperature, T
, and a dictionary of log equilibrium constants in partial pressure form, return the base-10 log equilibrium constant in number density form, i.e. log10(nK)
where nK = n(A)n(B)/n(AB)
.
Korg.get_mass
— Methodget_mass(f::Formula)
Returns the mass [g] of f
.
Korg.get_metals_H
— Methodget_metals_H(A_X; kwargs...)
Calculate [metals/H] given a vector, A_X
of absolute abundances, $A(X) = \log_{10}(n_M/n_\mathrm{H})$. See also get_alpha_H
.
Keyword Arguments
solar_abundances
(default:Korg.asplund_2020_solar_abundances
) is the set of solar abundances to use, as a vector indexed by atomic number.Korg.asplund_2009_solar_abundances
,Korg.grevesse_2007_solar_abundances
, andKorg.magg_2022_solar_abundances
are also provided for convenience.ignore_alpha
(default:true
): Whether or not to ignore the alpha elements when calculating [metals/H]. Iftrue
, [metals/H] is calculated using all elements heavier than He. Iffalse
, then both carbon and the alpha elements are ignored.alpha_elements
(default:Korg.default_alpha_elements
): vector of atomic numbers of the alpha elements. (Useful since conventions vary.)
Korg.get_number_densities
— Methodget_number_densities(atm::ModelAtmosphere) = [l.number_density for l in atm.layers]
This is a convienince functions for making plots, etc. Note that it doesn't access quantities in a memory-efficient order.
Korg.get_tau_5000s
— Methodget_tau_5000s(atm::ModelAtmosphere) = [l.tau_5000 for l in atm.layers]
This is a convienince functions for making plots, etc. Note that it doesn't access quantities in a memory-efficient order.
Korg.get_temps
— Methodget_temps(atm::ModelAtmosphere) = [l.temp for l in atm.layers]
This is a convienince functions for making plots, etc. Note that it doesn't access quantities in a memory-efficient order.
Korg.get_zs
— Methodget_zs(atm::ModelAtmosphere) = [l.z for l in atm.layers]
This is a convienince functions for making plots, etc. Note that it doesn't access quantities in a memory-efficient order.
Korg.greim_1960_Knm
— Methodgreim_1960_Knm(n, m)
Knm constants as defined by Griem 1960 for the long range Holtsmark profile. This function includes only the values for Brackett lines.
$K_{nm} = C_{nm} 2 π c / λ^2$ where $C_{nm} F = Δω$ and $F$ is the ion field. See Griem 1960 EQs 7 and 12. This works out to $K_{nm} = λ/F$.
Korg.holtsmark_profile
— Methodholtsmark_profile(β, P)
Calculates the Holtsmark profile for broadening of hydrogen lines by quasistatic charged particles. Adapted from SOFBET in HLINOP by Peterson and Kurucz. Draws heavily from Griem 1960.
Korg.hummer_mihalas_U_H
— Methodhummer_mihalas_U_H(T, nH, nHe, ne)
!!!note This is experimental, and not used by Korg for spectral synthesis.
Calculate the partition function of neutral hydrogen using the occupation probability formalism from Hummer and Mihalas 1988. See hummer_mihalas_w
for details.
Korg.hummer_mihalas_w
— Methodhummer_mihalas_w(T, n_eff, nH, nHe, ne; use_hubeny_generalization=false)
Calculate the correction, w, to the occupation fraction of a hydrogen energy level using the occupation probability formalism from Hummer and Mihalas 1988, optionally with the generalization by Hubeny+ 1994. (Sometimes Daeppen+ 1987 is cited instead, but H&M seems to be where the theory originated. Presumably it was delayed in publication.)
The expression for w is in equation 4.71 of H&M. K, the QM correction used in defined in equation 4.24. Note that H&M's "N"s are numbers (not number densities), and their "V" is volume. These quantities apear only in the form N/V, so we use the number densities instead.
This is based partially on Paul Barklem and Kjell Eriksson's WCALC fortran routine (part of HBOP.f), which is used by (at least) Turbospectrum and SME. As in that routine, we do consider hydrogen and helium as the relevant neutral species, and assume them to be in the ground state. All ions are assumed to have charge 1. Unlike that routine, the generalization to the formalism from Hubeny+ 1994 is turned off by default because I haven't closely checked it. The difference effects the chargedterm only, and temperature is only used when `usehubeny_generalizationis set to
true`.
Korg.hydrogen_line_absorption!
— Methodhydrogen_line_absorption!(αs, λs, T, nₑ, nH_I, UH_I, ξ, window_size; kwargs...)
Calculate contribution to the the absorption coefficient, αs, from hydrogen lines in units of cm^-1, at wavelengths λs
(a vector of ranges).
Uses profiles from Stehlé & Hutcheon (1999), which include Stark and Doppler broadening. For Halpha, Hbeta, and Hgamma, the p-d approximated profiles from Barklem, Piskunovet, and O'Mara 2000 are added to the absortion coefficient. This "convolution by summation" is inexact, but true convolution is expensive.
Arguments:
αs
: the absorption coefficient [cm^-1] vector into which to add the line absorptionλs
: the wavelengths at which to calculate the absorption [cm]T
: temperature [K]nₑ
: electron number density [cm^-3]nH_I
: neutral hydrogen number density [cm^-3]UH_I
: the value of the neutral hydrogen partition functionξ
: microturbulent velocity [cm/s]. This is only applied to Hα-Hγ. Other hydrogen lines profiles are dominated by stark broadening, and the stark broadened profiles are pre-convolved with a doppler profile.window_size
: the max distance from each line center [cm] at which to calculate the Stark and self-broadening profiles absorption for Hα-Hγ (those dominated by self-broadening).
Keyword arguments:
stark_profiles
(default:Korg._hline_stark_profiles
): tables from which to interpolate Stark profilesuse_MHD
: whether or not to use the Mihalas-Daeppen-Hummer formalism to adjust the occupation probabilities of each hydrogen orbital for plasma effects. Default:true
.
Korg.interpolate_marcs
— Methodinterpolate_marcs(Teff, logg, A_X; kwargs...)
interpolate_marcs(Teff, logg, m_H=0, alpha_m=0, C_m=0; kwargs...)
Returns a model atmosphere computed by interpolating models from MARCS ((Gustafsson+ 2008)[https://ui.adsabs.harvard.edu/abs/2008A&A...486..951G/abstract]). Along with Teff
and logg
, the atmosphere is specified by m_H
, alpha_m
, and C_m
, which can be automatically determined from an A_X
abundance vector (the recommended method, see format_A_X
). Note that the MARCS atmosphere models were constructed with the Grevesse+ 2007 solar abundances (Korg.grevesse_2007_solar_abundances
). This is handled automatically when A_X
is provided.
interpolate_marcs
uses three different interpolation schemes for different stellar parameter regimes. In the standard case the model atmosphere grid is the one generated for SDSS, transformed and linearly interpolated. For cool dwarfs (Teff
≤ 4000 K, logg
≥ 3.5), the grid is resampled onto unchanging tau_5000
values and interpolated with a cubic spline. For low-metallicity models (-5 ≤ m_H
< -2.5), a grid of standard composition (i.e. fixed alpha and C) atmospheres is used. (The microturbulence is 1km/s for dwarfs and 2km/s for giants and the mass for spherical models is 1 solar mass.) The interpolation method is the same as in the standard case. See Wheeler+ 2024 for more details and a discussion of errors introduced by model atmosphere interpolation. (Note that the cubic scheme for cool dwarfs is referred to as not-yet-implemented in the paper but is now available.)
keyword arguments
spherical
: whether or not to return a ShellAtmosphere (as opposed to a PlanarAtmosphere). By default true whenlogg
< 3.5.solar_abundances
: (default:grevesse_2007_solar_abundances
) The solar abundances to use whenA_X
is provided instead ofM_H
,alpha_M
, andC_M
. The default is chosen to match that of the atmosphere grid, and if you change it you are likely trying to do something else.clamp_abundances
: (default:false
) allowed only when specifyingA_X
. Whether or not to clamp the abundance parameters to be within range to avoid throwing an out of bounds error.perturb_at_grid_values
(default:true
): whether or not to add or a subtract a very small number to each parameter which is exactly at a grid value. This prevents null derivatives, which can cause problems for minimizers.resampled_cubic_for_cool_dwarfs
(default:true
): whether or not to used specialized method for cool dwarfs.archives
: A tuple containing the atmosphere grids to use. For testing purposes.
Korg.interpolate_molecular_cross_sections!
— Methodinterpolate_molecular_cross_sections!(α, molecular_cross_sections, λs, Ts, vmic, number_densities)
Interpolate the molecular cross-sections and add them to the total absorption coefficient α
. See MolecularCrossSection
for more information.
Korg.inverse_gaussian_density
— Methodinverse_gaussian_density(ρ, σ)
Calculate the inverse of a (0-centered) Gaussian PDF with standard deviation σ
, i.e. the value of x
for which ρ = exp(-0.5 x^2/σ^2}) / √[2π]
, which is given by σ √[-2 log (√[2π]σρ)]
. Returns 0 when ρ is larger than any value taken on by the PDF.
See also: inverse_lorentz_density
.
Korg.inverse_lorentz_density
— Methodinverse_lorentz_density(ρ, γ)
Calculate the inverse of a (0-centered) Lorentz PDF with width γ
, i.e. the value of x
for which ρ = 1 / (π γ (1 + x^2/γ^2))
, which is given by √[γ/(πρ) - γ^2]
. Returns 0 when ρ is larger than any value taken on by the PDF.
See also: inverse_gaussian_density
.
Korg.ismolecule
— Methodismolecule(f::Formula)
true
when f
is composed of more than one atom
Korg.lazy_multilinear_interpolation
— Methodlazy_multilinear_interpolation(params, nodes, grid; kwargs...)
This function is does multidimensional linear interpolation on a grid where the first two dimensions represent the values being interpolated. In other words, it interpolates matrices. It is written to minimize the number of reads from the grid, and access things in an efficient order. It is much faster than Interpolations.linear_interpolation when applied to memory-mapped grids. At present it is used only for model atmosphere and departure coefficient interpolation.
See also: interpolate_marcs
, Korg.CubicSplines.CubicSpline
Korg.line_absorption!
— Methodline_absorption!(α, linelist, λs, temp, nₑ, n_densities, partition_fns, ξ
; α_cntm=nothing, cutoff_threshold=1e-3, window_size=20.0*1e-8)
Calculate the opacity coefficient, α, in units of cm^-1 from all lines in linelist
, at wavelengths λs
[cm^-1].
other arguments:
temp
the temerature in K (as a vector, for multiple layers, if you like)n_densities
, a Dict mapping species to absolute number density in cm^-3 (as a vector, if temp is a vector).partition_fns
, a Dict containing the partition function of each speciesξ
is the microturbulent velocity in cm/s (n.b. NOT km/s)α_cntm
is as a callable returning the continuum opacity as a function of wavelength. The window within which a line is calculated will extend to the wavelength at which the Lorentz wings or Doppler core of the line are atcutoff_threshold * α_cntm[line.wl]
, whichever is greater.
Keyword Arguments
cuttoff_threshold
(default: 3e-4): seeα_cntm
verbose
(default: false): if true, show a progress bar.
Korg.line_profile
— Methodline_profile(λ₀, σ, γ, amplitude, λ)
A voigt profile centered on λ₀ with Doppler width σ (NOT √[2] σ, as the "Doppler width" is often defined) and Lorentz HWHM γ evaluated at λ
(cm). Returns values in units of cm^-1.
Korg.load_ExoMol_linelist
— Methodload_ExoMol_linelist(specs, states_file, transitions_file, upper_wavelength, lower_wavelength)
Load a linelist from ExoMol. Returns a vector of Line
s, the same as read_linelist
.
Arguments
spec
: the species, i.e. the molecule that the linelist is forstates_file
: the path to the ExoMol states filetransitions_file
: the path to the ExoMol, transitions fileupper_wavelength
: the upper limit of the wavelength range to load (Å)lower_wavelength
: the lower limit of the wavelength range to load (Å)
Keyword Arguments
line_strength_cutoff
: the cutoff for the line strength (default: -15) used to filter the linelist. Seeapproximate_line_strength
for more information.T_line_strength
: the temperature (K) at which to evaluate the line strength (default: 3500.0)
This functionality is in beta.
Korg.load_atomic_partition_functions
— Functionload_atomic_partition_functions()
Loads saved tabulated values for atomic partition functions from disk. Returns a dictionary mapping species to interpolators over log(T).
Korg.load_exomol_partition_functions
— Methodload_exomol_partition_functions()
Loads the exomol partition functions for polyatomic molecules from the HDF5 archive. Returns a dictionary mapping species to interpolators over log(T).
Korg.merge_bounds
— Functionmerge_bounds(bounds, merge_distance=0.0)
Sort a vector of lower-bound, upper-bound pairs and merge overlapping ranges.
Returns a pair containing:
- a vector of merged bounds
- a vector of vectors of indices of the original bounds which were merged into each merged bound
Korg.merge_close_lines
— Methodmerge_close_lines(linelist; merge_distance=0.2)
Produce a list of the species and wavelengths of the lines in linelist
, merging lines of the same species that are within merge_distance
(default: 0.2 Å). This is useful for labeling lines in a plot after running prune_linelist
.
Arguments
linelist
: the linelist (a vector ofLine
objects)
Keyword Arguments
merge_distance=0.2
: The maximum distance in Å between lines of the same species to be merged into a single entry
Returns
A vector of tuples (wl, wl_low, wl_high, species)
where wl
is gf-weighted wavelength of each set of merged lines (Å), wl_low
and wl_high
are their highest and lowest wavelength, and species
is a string (not a Korg.Species
) identifying the species of the line. These will be in wavelength order.
Korg.n_atoms
— Methodn_atoms(x)
The number of atoms in the Korg.Species or Korg.Formula x.
Korg.prune_linelist
— Methodprune_linelist(atm, linelist, A_X, wls...; threshold=1.0, sort=true, synthesis_kwargs...)
Return the vector containing the strongest lines in linelist
, (optionally) sorted by approximate equivalent width.
Arguments
atm
: the atmosphere modellinelist
: the linelist (a vector ofLine
objects)A_X
: the abundance of each element (seeformat_A_X
)wls...
: the wavelength ranges to synthesize over. These are specified the same way as thewls
forsynthesize
.
Keyword Arguments
threshold=0.1
: The threshold ratio in the line center absorption to the continuum absorption computed at the photosphere for a line to be included in the returned list.0.1
is a reasonable default for getting a sense of what might be measurable in a high-res, high-quality spectrum, but it should not be used to create a linelist for synthesis.sort_by_EW=true
: Iftrue
, the returned linelist will be sorted by approximate reduced equivalent width. Iffalse
, the linelist will be in wavelength order. Leaving the list in wavelength order is much faster, but sorting by strength is useful for visualizing the strongest lines.verbose=true
: Iftrue
, a progress bar will be displayed while measuring the EWs. All other kwargs are passed to internal calls tosynthesize
.max_distance=0.0
, how far fromwls
lines can be (in Å) before they are excluded from the returned list.
While this function can be used to prune a linelist for synthesis, the default behavior too aggressive for this purpose. Set a much lower threshold (e.g. threshold=1e-4
) and use sort_by_EW=false
if you are pruning the linelist to speedup synthesis. Note that Korg will dynamically choose which lines to include even if you use a large linelist (see the line_cutoff_threshold
keyword argument to synthesize
).
See also merge_close_lines
if you are using this for plotting.
Korg.read_Barklem_Collet_logKs
— Methodread_Barklem_Collet_logKs(fname)
Reads the equilibrium constants from the HDF5 file produced by the Barklem and Collet 2016 paper. Returns a Dict from Korg.Species to Korg.CubicSplines from ln(T) to log10(K).
As recommended by Aquilina+ 2024, we modify the C2 equilibrium constant reflect the dissociation energy reported by Visser+ 2019.
Korg.read_Barklem_Collet_table
— Methodfunction read_Barklem_Collet_table(fname; transform=identity)
Constructs a Dict holding tables containing partition function or equilibrium constant values across ln(temperature). Applies transform (which you can use to, e.g. change units) to each example.
Korg.read_linelist
— Methodread_linelist(filename; format="vald", isotopic_abundances=Korg.isotopic_abundances)
Parse a linelist file, returning a vector of Line
s.
The format
keyword argument can be used to specify one of these linelist formats (default: "vald"
):
"vald"
for a VALD linelist. These can be either "short" or "long" format, "extract all" or "extract stellar". Air wavelengths will automatically be converted into vacuum wavelengths, and energy levels will be automatically converted from cm$^{-1}$ to eV."kurucz"
for an atomic or molecular Kurucz linelist (format=kurucz_vac if it uses vacuum wavelengths; be warned that Korg will not assume that wavelengths are vacuum below 2000 Å),"moog"
for a MOOG linelist (doesn't support broadening parameters or dissociation energies, assumed to be in vacuum wavelengths)."moog_air"
for a MOOG linelist in air wavelengths."turbospectrum"
for a Turbospectrum linelist in air wavelengths. Note that Korg doesn't make use of the (optional) orbital angular momentum quantum number, l, for the upper or lower levels, so it won't fall back on generic ABO recipes when the ABO parameters are not available. Korg's interpretation of thefdamp
parameter is also slightly different from Turbospectrum's. See the documentation of thevdW
parameter ofLine
for details. Korg will error if encounters an Unsoeld fudge factor, which it does not support.- "turbospectrum_vac" for a Turbospectrum linelist in vacuum wavelengths.
- "korg" for a Korg linelist (saved with hdf5). If the filename ends in
.h5
, this will be used by default.
For VALD and Turbospectrum linelists with isotope information available, Korg will scale log gf values by isotopic abundance (unless VALD has already pre-scaled them), using isotopic abundances from NIST ([Korg.isotopicabundances]). To use custom isotopic abundances, just pass `isotopicabundances` with the same structure: a dict mapping atomic number to a dict mapping from atomic weight to abundance.
Be warned that for linelists which are pre-scaled for isotopic abundance, the estimation of radiative broadening from log(gf) is not accurate.
See also: load_ExoMol_linelist
, save_linelist
.
Korg.read_model_atmosphere
— Methodread_model_atmosphere(filename)
Parse the provided model atmosphere file in MARCS ".mod" format. Returns either a PlanarAtmosphere
or a ShellAtmosphere
.
Korg.read_molecular_cross_section
— Methodread_molecular_cross_section(filename)
Read a precomputed molecular cross-section from a file created by save_molecular_cross_section
. See also MolecularCrossSection
.
Korg.saha_ion_weights
— Methodsaha_ion_weights(T, nₑ, atom, ionization_energies, partition_functions)
Returns (wII, wIII)
, where wII
is the ratio of singly ionized to neutral atoms of a given element, and wIII
is the ration of doubly ionized to neutral atoms.
arguments:
- temperature
T
[K] - electron number density
nₑ
[cm^-3] - atom, the atomic number of the element
ionization_energies
is a collection indexed by integers (e.g. aVector
) mapping elements' atomic numbers to their first three ionization energiespartition_funcs
is aDict
mapping species to their partition functions
Korg.save_linelist
— Methodsave_linelist(path, linelist)
Save a Korg linelist (a Vector{Line}
) to an HDF5 file which can be read by read_linelist
.
Korg.save_molecular_cross_section
— Methodsave_molecular_cross_section(filename, cross_section)
Save a precomputed molecular cross-section to a file. See also MolecularCrossSection
, read_molecular_cross_section
.
Korg.scaled_stark
— Methodthe stark broadening gamma scaled acording to its temperature dependence
Korg.scaled_vdW
— Methodscaled_vdW(vdW, m, T)
The vdW broadening gamma scaled acording to its temperature dependence, using either simple scaling or ABO. See Anstee & O'Mara (1995) or Paul Barklem's notes for the definition of the ABO γ.
vdW
should be either γ_vdW
evaluated at 10,000 K, or tuple containing the ABO params (σ, α)
. The species mass, m
, is ignored in the former case.
Korg.setup_ionization_energies
— Functionsetup_ionization_energies([filename])
Parses the table of ionization energies and returns it as a dictionary mapping elements to their ionization energies, [χ₁, χ₂, χ₃]
in eV.
Korg.setup_partition_funcs_and_equilibrium_constants
— Methodsetup_partition_funcs_and_equilibrium_constants()
Returns two dictionaries. One holding the default partition functions, and one holding the default log10 equilibrium constants.
Default partition functions
The partition functions are custom (calculated from NIST levels) for atoms, from Barklem & Collet 2016 for diatomic molecules, and from exomol for polyatomic molecules. For each molecule, we include only the most abundant isotopologue.
Note than none of these partition functions include plasma effects, e.g. via the Mihalas Hummer Daeppen occupation probability formalism. They are for isolated species. This can lead to a couple percent error for neutral alkalis and to greater errors for hydrogen in some atmospheres, particularly those of hot stars.
Default equilibrium constants
Molecules have equilibrium constants in addition to partition functions. For the diatomics, these are provided by Barklem and Collet, which extensively discusses the dissociation energies. For polyatomics, we calculate these ourselves, using atomization energies calculated from the enthalpies of formation at 0K from NIST's CCCDB.
Korg's equilibrium constants are in terms of partial pressures, since that's what Barklem and Collet provide.
Korg.sigma_line
— Methodsigma_line(wl)
The cross-section (divided by gf) at wavelength wl
in Ångstroms of a transition for which the product of the degeneracy and oscillator strength is 10^log_gf
.
Korg.subspectrum_indices
— Methodsubspectrum_indices(wls::Wavelengths)
Returns a vector of Julia ranges, which can be used to index into the full spectrum to get the sub-spectrum corresponding to each wavelength range in wls
.
Korg.synth
— Methodsynth(kwargs...)
This function creates a synthetic spectrum. It's easier to use than synthesize
, but it gives you less control. Unlike synthesize
, it returns a tuple of (wavelengths, rectified_flux, cntm)
(Wavelength in Å, rectified flux as a unitless number between 0 and 1, and continuum in erg/s/cm^5). Korg.synth
also provides shortcuts for some ways you might want to post-process the spectrum (applying a LSF, rotation, etc).
Keyword arguments
Teff
: effective temperature in K (default: 5000)logg
: surface gravity in cgs units (default: 4.5)m_H
: metallicity, [metals/H], (default: 0.0) (Seeformat_A_X
for precisely how this is interpreted.)alpha_H
: alpha enhancement, [α/H], (default:m_H
) (Seeformat_A_X
for precisely how this is interpreted.)- Any atomic symbol (e.g.
Fe
orC
) can be used to to specify a (solar relative, [X/H]) abundance. These overridem_H
andalpha_H
. Specifying an individual abundance means that the true metallicity and alpha will not correspond precisely to the values ofm_H
andalpha_H
. Seeformat_A_X
for details. linelist
: a linelist, (default:get_VALD_solar_linelist()
). See alsoread_linelist
.wavelengths
: a tuple of the start and end wavelengths (default: (5000, 6000)), or a vector of(λstart, λstop)
pairs. SeeWavelengths
for all the ways the wavelengths can be specified.rectify
: whether to rectify (continuum normalize) the spectrum (default: true)R
: resolution (default:Inf
, no LSF applied).R
can be a scalar, or a function from wavelength (in Å) to resolving power. Seeapply_LSF
for details on how to do this manually.vsini
: projected rotational velocity in km/s (default: 0). Seeapply_rotation
for details on how to do this manually.vmic
: microturbulent velocity in km/s (default: 1.0).synthesize_kwargs
: additional keyword arguments pass tosynthesize
.format_A_X_kwargs
: additional keyword arguments pass toformat_A_X
.
Korg.synthesize
— Methodsynthesize(atm, linelist, A_X, (λ_start, λ_stop); kwargs... )
synthesize(atm, linelist, A_X, wavelength_ranges; kwargs... )
Compute a synthetic spectrum. Returns a SynthesisResult
.
Arguments
atm
: the model atmosphere (seeread_model_atmosphere
)linelist
: A vector ofLine
s (seeread_linelist
,get_APOGEE_DR17_linelist
, andget_VALD_solar_linelist
).A_X
: a vector containing the A(X) abundances (log(X/H) + 12) for elements from hydrogen to uranium. (seeformat_A_X
)- The wavelengths at which to synthesize the spectrum. They can be specified either as a pair
(λstart, λstop)
, or as a list of pairs[(λstart1, λstop1), (λstart2, λstop2), ...]
. λ_start
: the lower bound (in Å) of the region you wish to synthesize.λ_stop
: the upper bound (in Å) of the region you wish to synthesize.λ_step
(default: 0.01): the (approximate) step size to take (in Å).
Example
to synthesize a spectrum between 5000 Å and 5100 Å, with all metal abundances set to 0.5 dex less than the solar value except carbon, except carbon, which we set to [C/H]=-0.25:
atm = read_model_atmosphere("path/to/atmosphere.mod")
linelist = read_linelist("path/to/linelist.vald")
A_X = format_A_X(-0.5, Dict("C" => -0.25))
result = synthesize(atm, linelist, A_X, 5000, 5100)
Optional arguments:
vmic
(default: 0) is the microturbulent velocity, $\xi$, in km/s.line_buffer
(default: 10): the farthest (in Å) any line can be from the provided wavelength range before it is discarded. If the edge of your window is near a strong line, you may have to turn this up.cntm_step
(default 1): the distance (in Å) between point at which the continuum opacity is calculated.hydrogen_lines
(default:true
): whether or not to include H lines in the synthesis.use_MHD_for_hydrogen_lines
(default:true
): whether or not to use the MHD occupation probability formalism for hydrogen lines. (MHD is always used for hydrogen bound-free absorption.)hydrogen_line_window_size
(default: 150): the maximum distance (in Å) from each hydrogen line center at which to calculate its contribution to the total absorption coefficient.mu_values
(default: 20): the number of μ values at which to calculate the surface flux, or a vector of the specific values to use when doing transfer in spherical geometry. Ifmu_points
is an integer, the values are chosen per Gauss-Legendre integration. If they are specified directly, the trapezoid rule is used for the astrophysical flux. The default values is sufficient for accuracy at the 10^-3 level. Note that if you are using the default radiative transfer scheme, with a plane-parallel model atmosphere, the integral over μ is exact, so this parameter has no effect. The points and weights are returned in themu_grid
field of the output.line_cutoff_threshold
(default:3e-4
): the fraction of the continuum absorption coefficient at which line profiles are truncated. This has major performance impacts, since line absorption calculations dominate more syntheses. Turn it down for more precision at the expense of runtime. The default value should effect final spectra below the 10^-3 level.electron_number_density_warn_threshold
(default:Inf
): if the relative difference between the calculated electron number density and the input electron number density is greater than this value, a warning is printed. By default, this warning is suppress (threshold isInf
) because it is very easily raised in cases where it is of no observable consequence. See alsoelectron_number_density_warn_min_value
, below.electron_number_density_warn_min_value
(default:1e-4
): The minimum value of the ratio of the electron number density to the total number density at which a warning is printed.return_cntm
(default:true
): whether or not to return the continuum at each wavelength. If this is false,solution.cntm
will benothing
.ionization_energies
, aDict
mappingSpecies
to their first three ionization energies, defaults toKorg.ionization_energies
.partition_funcs
, aDict
mappingSpecies
to partition functions (in terms of ln(T)). Defaults to data from Barklem & Collet 2016,Korg.default_partition_funcs
.equilibrium_constants
, aDict
mappingSpecies
representing diatomic molecules to the base-10 log of their molecular equilibrium constants in partial pressure form. Defaults to data from Barklem and Collet 2016,Korg.default_log_equilibrium_constants
.use_chemical_equilibrium_from
(default:nothing
): Takes another solution returned bysynthesize
. When provided, the chemical equilibrium solution will be taken from this object, rather than being recomputed. This is physically self-consistent only when the abundances,A_X
, and model atmosphere,atm
, are unchanged.molecular_cross_sections
(default:[]
): A vector of precomputed molecular cross-sections. SeeMolecularCrossSection
for how to generate these. If you are using the default radiative transfer scheme, your molecular cross-sections should cover 5000 Å only if your linelist does.tau_scheme
(default: "linear"): how to compute the optical depth. Options are "linear" and "bezier" (testing only–not recommended).I_scheme
(default:"linear_flux_only"
): how to compute the intensity. Options are"linear"
,"linear_flux_only"
, and"bezier"
."linear_flux_only"
is the fastest, but does not return the intensity values anywhere except at the top of the atmosphere. "linear" performs an equivalent calculation, but stores the intensity at every layer."bezier"
is for testing and not recommended.verbose
(default:false
): Whether or not to print information about progress, etc.
Korg.translational_U
— Methodtranslational_U(m, T)
The (possibly inverse) contribution to the partition function from the free movement of a particle. Used in the Saha equation.
arguments
m
is the particle massT
is the temperature in K
Korg.vacuum_to_air
— Methodvacuum_to_air(λ; cgs=λ<1)
convert λ from a vacuum to air. λ is assumed to be in Å if it is ⩾ 1, in cm otherwise. Formula from Birch and Downs (1994) via the VALD website.
See also: air_to_vacuum
.
Korg.voigt_hjerting
— Methodvoigt_hjerting(α, v)
The Hjerting function, $H$, somtimes called the Voigt-Hjerting function. $H$ is defined as H(α, v) = ∫^∞_∞ exp(-y^2) / ((u-y)^2 + α^2) dy
(see e.g. the unnumbered equation after Gray equation 11.47). It is equal to the ratio of the absorption coefficient to the value of the absorption coefficient obtained at the line center with only Doppler broadening.
If x = λ-λ₀
, Δλ_D = σ√2
is the Doppler width, and Δλ_L = 4πγ
is the Lorentz width,
voigt(x|Δλ_D, Δλ_L) = H(Δλ_L/(4πΔλ_D), x/Δλ_D) / (Δλ_D√π)
= H(γ/(σ√2), x/(σ√2)) / (σ√(2π))
Approximation from Hunger 1965.
Korg.λ_to_ν_bound
— Methodλ_to_ν_bound(λ_bound)
Converts a λ Inverval
(in cm) to an equivalent ν Interval
(in Hz), correctly accounting for tricky floating point details at the bounds.
Korg.Fit
— ModuleFunctions for fitting to data.
This submodule is in beta. It's API may change.
Korg.Fit.ews_to_abundances
— Methodews_to_abundances(atm, linelist, A_X, measured_EWs; kwargs... )
Functions using equivalent widths have been dissabled while we fix some bugs.
Compute per-line abundances on the linear part of the curve of growth given a model atmosphere and a list of lines with equivalent widths.
Arguments:
atm
: the model atmosphere (seeKorg.read_model_atmosphere
andKorg.interpolate_marcs
).linelist
: A vector ofKorg.Line
s (seeKorg.read_linelist
). The lines must be sorted by wavelength.A_X
: a vector containing the A(X) abundances (log(nX/nH) + 12) for elements from hydrogen to uranium (seeKorg.format_A_X
). All syntheses are done with these abundances, so if the resulting abundances deviate significantly from these, you may wish to iterate.measured_EWs
: a vector of equivalent widths (in mÅ)
Returns
A vector of abundances (A(X) = log10(n_X/n_H) + 12
format) for each line in linelist
.
Optional arguments:
wl_step
(default: 0.01) is the resolution in Å at which to synthesize the spectrum around each line.ew_window_size
(default: 2): the farthest (in Å) to consider equivalent width contributions for each line. It's very important that this is large enough to include each line entirely.blend_warn_threshold
(default: 0.01) is the minimum absorption between two lines allowed before triggering a warning that they may be blended. All other keyword arguments are passed toKorg.synthesize
when synthesizing each line.
Korg.Fit.ews_to_stellar_parameters
— Functionews_to_stellar_parameters(linelist, measured_EWs, [measured_EW_err]; kwargs...)
Functions using equivalent widths have been dissabled while we fix some bugs.
Find stellar parameters from equivalent widths the "old fashioned" way. This function finds the values of $T_\mathrm{eff}$, $\log g$, $v_{mic}$, and [m/H] which satisfy the following conditions (using a Newton-Raphson solver):
- The slope of the abundances of neutral lines with respect to lower excitation potential is zero.
- The difference between the mean abundances of neutral and ionized lines is zero.
- The slope of the abundances of neutral lines with respect to reduced equivalent width is zero.
- The difference between the mean abundances of all lines and the model-atmosphere input [m/H] is zero. Here the "slope" refers to the slope of a linear fit to the abundances of the lines in question.
Arguments:
linelist
: A vector ofKorg.Line
objects (seeKorg.read_linelist
). The lines must be sorted by wavelength.measured_EWs
: a vector of equivalent widths (in mÅ).measured_EW_err
(optional): the uncertainty inmeasured_EWs
. If not specified, all lines are assumed to have the same uncertainty. These uncertainties are used when evaluating the equations above, and are propagated to provide uncertainties in the resulting parameters.
Returns:
A tuple containing:
- the best-fit parameters:
[Teff, logg, vmic, [m/H]]
as a vector - the statistical uncertainties in the parameters, propagated from the uncertainties in the equivalent widths. This is zero when EW uncertainties are not specified.
- the systematic uncertainties in the parameters, estimated from the scatter the abundances computed from each line not accounted for by the EW uncertainties.
Keyword arguments:
Teff0
(default: 5000.0) is the starting guess for Tefflogg0
(default: 3.5) is the starting guess for loggvmic0
(default: 1.0) is the starting guess for vmic. Note that this must be nonzero in order to avoid null derivatives. Very small values are fine.m_H0
(default: 0.0) is the starting guess for [m/H]tolerances
(default:[1e-3, 1e-3, 1e-3, 1e-3]
) is the tolerance for the residuals each equation listed above. The solver stops when all residuals are less than the corresponding tolerance.max_step_sizes
(default:[1000.0, 1.0, 0.3, 0.5]
) is the maximum step size to take in each parameter direction. This is used to prevent the solver from taking too large of a step and missing the solution. Be particularly cautious with the vmic (third) parameter, as the unadjusted step size is often too large.parameter_ranges
(default:[(2800.0, 8000.0), (-0.5, 5.5), (1e-3, 10.0), (-2.5, 1.0)]
) is the allowed range for each parameter. This is used to prevent the solver from wandering into unphysical parameter space, or outside the range of the MARCS grid supported byKorg.interpolate_marcs
. The default ranges $T_\mathrm{eff}$, $\log g$, and [m/H] are the widest supported by the MARCS grid. Note that vmic must be nonzero in order to avoid null derivatives.callback
: is a function which is called at each step of the optimization. It is passed three arguments:- the current values of the parameters
- the residuals of each equation being solved
- the abundances of each line computed with the current parameters. You can pass a callback function, to e.g. make a plot of the residuals at each step.
max_iterations
(default: 30) is the maximum number of iterations to allow before stopping the optimization.
Korg.Fit.fit_spectrum
— Functionfit_spectrum(obs_wls, obs_flux, obs_err, linelist, initial_guesses, fixed_params; kwargs...)
Find the parameters and abundances that best match a rectified observed spectrum.
Arguments:
obs_wls
: the wavelengths of the observed spectrum in Å. These must be vacuum wavelengths.obs_flux
: the rectified flux of the observed spectrumobs_err
: uncertainty influx
linelist
: a linelist to use for the synthesisinitial_guesses
: a NamedTuple specifying initial guesses for the parameters to be fit. See "Specifying parameters" below.fixed_params
: a NamedTuple specifying parameters to be held fixed. See "Specifying parameters" below.
initial_guesses
and fixed_params
can also be specified as Dicts instead of NamedTuples, which is more convenient when calling Korg from python.
Specifying parameters
Parameters are specified as named tuples or dictionaries. Named tuples look like this: (Teff=5000, logg=4.5, m_H=0.0)
. Single-element named tuples require a semicolon: (; Teff=5000)
.
Required parameters
Teff
and logg
must be specified in either initial_guesses
or fixed_params
.
Optional Parameters
These can be specified in either initial_guesses
or fixed_params
, but if they are not default values are used.
m_H
: the metallicity of the star, in dex. Default:0.0
alpha_H
: the alpha enhancement of the star, in dex. Default:m_H
. Note that, because of the parameter range supported byKorg.interpolate_marcs
, only values within ±1 ofm_H
are supported.vmic
: the microturbulence velocity, in km/s. Default:1.0
vsini
: the projected rotational velocity of the star, in km/s. Default:0.0
. SeeKorg.apply_rotation
for details.epsilon
: the linear limb-darkening coefficient. Default:0.6
. Used for applying rotational broadening only. SeeKorg.apply_rotation
for details.- Individual elements, e.g.
Na
, specify the solar-relative ([X/H]) abundance of that element.
Keyword arguments
R
, the resolution of the observed spectrum. This is required. It can be specified as a function of wavelength, in which case it will be evaluated at the observed wavelengths.windows
is a vector of wavelength pairs, each of which specifies a wavelength "window" to synthesize and contribute to the total χ². If not specified, the entire spectrum is used. Overlapping windows are automatically merged.adjust_continuum
(default:false
) if true, adjust the continuum with the best-fit linear correction within each window, minimizing the chi-squared between data and model at every step of the optimization.wl_buffer
is the number of Å to add to each side of the synthesis range for each window.time_limit
is the maximum number of seconds to spend in the optimizer. (default:10_000
). The optimizer will only checks against the time limit after each step, so the actual wall time may exceed this limit.precision
specifies the tolerance for the solver to accept a solution. The solver operates on transformed parameters, soprecision
doesn't translate straightforwardly to Teff, logg, etc, but the default value,1e-4
, provides a theoretical worst-case tolerance of about 0.15 K inTeff
, 0.0002 inlogg
, 0.0001 inm_H
, and 0.0004 in detailed abundances. In practice the precision achieved by the optimizer is about 10x bigger than this.postprocess
can be used to arbitrarilly transform the synthesized (and LSF-convolved) spectrum before calculating the chi2. It should take the formpostprocess(flux, data, err)
and write its changes in-place to the flux array.LSF_matrix
: this can be provedided along withsynthesis_wls
in place of specifyingR
if you have a precomputed custom LSF matrix.synthesis_wls
: seeLSF_matrix
above. This can be a Korg.Wavelengths object or any arguments that can be passed to its constructor, e.g. a range or vector of ranges. Wavelengths are in Å.- Any additional keyword arguments will be passed to
Korg.synthesize
when synthesizing the spectra for the fit.
Returns
A NamedTuple with the following fields:
best_fit_params
: the best-fit parametersbest_fit_flux
: the best-fit flux, with LSF applied, resampled, and rectified.obs_wl_mask
: a bitmask forobs_wls
which selects the wavelengths used in the fit (i.e. those in thewindows
)solver_result
: the result object fromOptim.jl
trace
: a vector of NamedTuples, each of which contains the parameters at each step of the optimization. This is empty for single parameter fits, because the underlying solver doesn't supply it.covariance
: a pair(params, Σ)
whereparams
is vector of parameter name (providing an order), andΣ
is an estimate of the covariance matrix of the parameters. It is the approximate inverse hessian of the log likelihood at the best-fit parameter calculated by the BGFS algorithm, and should be interpreted with caution.
This function takes a long time to compile the first time it is called. Compilation performance is significantly better on Julia 1.10+ than previous versions, so if you are using an older version of Julia, you may want to upgrade.
Korg.Fit.linear_continuum_adjustment!
— Methodlinear_continuum_adjustment!(obs_wls, windows, model_flux, obs_flux, obs_err)
Adjust the model flux to match the observed flux by fitting a line (as a function of wavelength) to the residuals, and dividing it out. This can compensate for poorly done continuum normalization.
Note, obs_wls must be masked.
Korg.Fit.postprocessed_synthetic_spectrum
— MethodSynthesize a spectrum, apply the LSF, and postprocess it, catching and potentially rethrowing errors. This is used by fit_spectrum
.
Korg.Fit.scale
— MethodRescale each parameter so that it lives on (-∞, ∞).
Korg.Fit.synthetic_spectrum
— MethodSynthesize a spectrum, returning the flux, with LSF applied, resampled, and rectified. This is an internal function oned by fitting routines. See Korg.synthesize
to synthesize spectra as a Korg user.
Korg.Fit.unscale
— MethodUnscale each parameter so that it lives on the appropriate range instead of (-∞, ∞).
Korg.Fit.validate_params
— MethodValidate fitting parameters, and insert default values when needed. Used by fit_spectrum
.
these can be specified in either initialguesses or fixedparams, but if they are not, these values are inserted into fixed_params
Korg.CubicSplines.CubicSpline
— MethodCubicSpline(xs, ys; extrapolate=false)
Construct a interpolant using xs
and ys
as the knot coordinates. Assumes xs
is sorted. Apply this object as a function to interpolate at any x value in the domain. If extrapolate
is false, x values outside [xs[1]
, xs[end]
] throw errors, if extrapolate
is true, the interpolant uses flat extrapolation, i.e. it returns the extreme value.
Korg.CubicSplines.cumulative_integral!
— Methodcumulative_integral!(out, A, t1, t2)
Given a curve described by the spine, A, Calculates the integral from t1 to t for all t = t1, t2, and all spline knots in between. So if t1
is A.t[1]
and t2
is A.t[end]
, out
should have the same length as A.t
.
Korg.ContinuumAbsorption.H2plus_bf_and_ff
— FunctionH2plus_bf_and_ff(ν, T, nH_I_div_partition, n_HII; kwargs...)
Compute the combined H₂⁺ bound-free and free-free linear absorption coefficient α using the tables from Stancil 1994. Note that this refers to interactions between free-protons and neutral hydrogen, not electrons and doubly ionized H2, i.e. the b-f interaction is photodissociation, not photoionization.
Arguments
ν::AbstractVector{<:Real}
: sorted frequency vector in HzT
: temperature in KnH_I
: the total number density of H I divided by its partition function.nH_II
: the number density of H II (not of H₂⁺).
For a description of the kwargs, see Continuum Absorption Kwargs.
Notes
This computes the H₂⁺ number density from those of H I and H II, since ionized molecules are not included in the molecular equlibrium calculation. Stancil provides approximate equilibrium constants, K, for the molecule, but the Barklem and Collet values used elsewhere by Korg may be more reliable. Once ionized molecules are fully supported, those values should be used instead. While cross sections are tabulated down to only 3150 K, the cross sections could be linearly interpolated 1000 K or so lower, if reliable K values are availble.
Because n(H₂⁺) is computed on the fly, the "cross-sections" used (internally) by this function have units of cm^-5, since they must be multiplied by n(H I) and n(H II) to obtain absorption coefficients.
Stancil, Bates (1952), and Gray (2005) all use n(H I), but Kurucz (1970) notes in section 5.2 that the photodisociation of H₂⁺ produces a ground state H atom and that we should therefore use n(H I, n=1) instead. This difference causes a descrepancy of a fraction of a percent (at most) for T ≤ 1.2e4 K. Here we use n(H I).
Korg.ContinuumAbsorption.H_I_bf
— MethodH_I_bf(νs, T, nH, nHe, ne, invU_H; n_max_MHD=6, use_hubeny_generalization=false,
taper=false, use_MHD_for_Lyman=false)
The bound-free linear absorption coefficient contributed by all energy states of a neutral Hydrogen atom. Even though the Mihalas-Hummer-Daeppen (MHD) occupation probability formalism is not used in Korg when computing the hydrogen partition function, it is used here. That means that series limit jumps (e.g. the Balmer jump), are "rounded off", as photons with less than the "classical" ionization energy can ionize if the upper level is dissolved into the continuum.
Required Arguments
νs
: sorted frequency vector in HzT
: temperature in KnH_I
: the total number density of neutral Hydrogen (in cm⁻³)nHe_I
: the total number density of neutral Helium (in cm⁻³)ne
: the number density of electrons (in cm⁻³)invU_H
: The inverse of the neutral hydrogen partition function (neglecting contributions from the MHD formalism)
For n=1 through n=n_max_MHD
(default: 6), the cross-sections are computed using Nahar 2021. These are modified using the MHD formalism to account for level dissolution. For larger n, the cross-sections are calculated with a simple analytic formula (see simple_hydrogen_bf_cross_section
).
Because MHD level dissolution applied to the the Lyman series limit leads to inflated cross-sections in the visible, we don't use MHD for bf absorption from n=1. This can be overridden by setting use_MHD_for_Lyman=true
, in which case you will also want to set taper=true
, which the same tapering of the cross-section as HBOP to fix the problem.
The use_hubeny_generalization
keyword argument enables the generalization of the MHD from Hubeny 1994. It is experimental and switched off by default.
Korg.ContinuumAbsorption.Heminus_ff
— FunctionHeminus_ff(ν, T, nHe_I_div_partition, ne; kwargs...)
Compute the He⁻ free-free opacity κ.
The naming scheme for free-free absorption is counter-inutitive. This actually refers to the reaction: photon + e⁻ + He I -> e⁻ + He I.
Arguments
ν::AbstractVector{<:Real}
: sorted frequency vector in HzT
: temperature in KnHe_I_div_partition
: the total number density of H I divided by its partition function.ne
: the number density of free electrons.
For a description of the kwargs, see Continuum Absorption Kwargs.
Notes
This uses the tabulated values from John (1994). The quantity K is the same used by Bell and Berrington (1987). See Hminus_ff
for an explanation.
According to John (1994), improved calculations are unlikely to alter the tabulated data for λ > 1e4Å, "by more than about 2%." The errors introduced by the approximations for 5.06e3 Å ≤ λ ≤ 1e4 Å "are expected to be well below 10%."
Korg.ContinuumAbsorption.Hminus_bf
— FunctionHminus_bf(ν, T, nH_I_div_partition, ne; kwargs...)
Compute the H⁻ bound-free linear absorption coefficient α, $\alpha_\nu = \sigma_{bf}(H^-) n(H⁻) (1 - \exp \left( \frac{-h\nu}{k T}\right))$
Arguments
ν::AbstractVector{<:Real}
: sorted frequency vector in HzT
: temperature in KnH_I_div_partition
: the total number density of H I divided by its partition function.ne
: the electron number density
This uses cross-sections from McLaughlin 2017. For a description of the kwargs, see Continuum Absorption Kwargs.
Notes
This function assumes that n(H⁻) ≪ n(H I) + n(H II). The number density of n(H⁻) isn't precomputed as part of Korg's molecular equlibrium, it's computed here instead.
The McLaughlin tabulated values (downloaded as a text file) contained a stray line out of monotonic order. A version of the data file with the line deleted is saved at data/McLaughlin2017Hminusbf.dat
for archival purposes. (Korg doesn't read this file, it reads data/McLaughlin2017Hminusbf.h5
.)
Korg.ContinuumAbsorption.Hminus_ff
— FunctionHminus_ff(ν, T, nH_I_div_partition, ne; kwargs...)
Compute the H⁻ free-free linear absorption coefficient α
The naming scheme for free-free absorption is counter-inutitive. This actually refers to the reaction: photon + e⁻ + H I -> e⁻ + H I.
Arguments
ν::AbstractVector{<:Real}
: sorted frequency vector in HzT
: temperature in KnH_I_div_partition::Flt
: the total number density of H I divided by its partition function.ne
: the number density of free electrons.
For a description of the kwargs, see Continuum Absorption Kwargs.
Notes
This is based on Table 1, in Bell & Berrington (1987), which tabulates values for "the H⁻ absorption coefficient", K (including the correction for stimulated emission). This quantity is in units of cm^4/dyn, and must be multiplied by the electron partial pressure and the ground-state neutral hydrogen number density to obtain a linear absorption coefficent, α.
The stipulation that the hydrogen should be ground-state only is based on the beginning of Section 2 in Bell and Berrington (1987), or alternately, Section 5.3 from Kurucz (1970). When Gray (2005) refers to this, it implicitly assumes that n(H I, n = 1) ≈ n(H I)
. Note that
n(H I, n = 1) = n(H I)*gₙ₌₁/U(T)*exp(-Eₙ₌₁/(k*T)) = n(H I) * 2/U(T)*exp(0) = n(H I) * 2/U(T).
Since U(T) ≈ 2 up to fairly large temperatures, this is not unreasonable.
Korg.ContinuumAbsorption._load_gauntff_table
— Function_load_gauntff_table([fname])
Returns a table of thermally-averaged free-free Gaunt factors, and the values of log₁₀(γ²) and log₁₀(u) associated with each point.
This loads the non-relativistic free-free data published by van Hoof et al. (2014).
Note: This function code could trivially be adapted to load the relativistic free-free gaunt factors published by van Hoof et al (2015).
Korg.ContinuumAbsorption._ndens_Hminus
— Function_ndens_Hminus(nH_I_div_partition, ne, T, ion_energy = _H⁻_ion_energy)
Compute the number density of H⁻ (implements eqn 5.10 of Kurucz 1970). This is an application of the saha equation where the "ground state" is H⁻ and the "first ionization state" is H I. The partition function of H⁻ is 1 at all temperatures.
Korg.ContinuumAbsorption.bounds_checked_absorption
— Methodbounds_checked_absorption(func; ν_bound, temp_bound)
Constructs a wrapped function that implements bounds checking and extrapolation
Parameters
func
: a function that has a signaturef(ν::Real, T::Real, args...)::Real
, whereν
is frequency (in Hz) andT
is temperature (in K)ν_bound::Interval
: Interval of frequencies (in Hz) over whichfunc
is valid.temp_bound::Interval
: Interval of temperatures (in K) over whichfunc
is valid.
The resulting function will have this signature:
wrapped_func(ν::AbstractVector{<:Real}, T::Real, args...; kwargs...)
Wrapped Function Parameters
ν::AbstractVector{<:Real}
: sorted vector of frequencies (in Hz)T::Real
: temperature (in K)args...
: function-specific arguments
For a description for the kwargs...
, see Continuum Absorption Kwargs. # if generating the inner function in this way involves too much overhead, we could:
Korg.ContinuumAbsorption.electron_scattering
— Methodelectron_scattering(nₑ)
Compute the linear absorption coefficient, α, from scattering off of free electrons. This has no wavelength dependence. It assumes isotropic scattering. (See, e.g. Gray p 160.)
Arguments
nₑ::F
: number density of free electrons (in cgs)
Korg.ContinuumAbsorption.gaunt_ff_vanHoof
— Methodgaunt_ff_vanHoof(log_u, log_γ2)
computes the thermally averaged, non-relativistic free-free gaunt factor by interpolating the table provided by van Hoof et al. (2014).
Arguments
log_u
: Equal to log₁₀(u) = log₁₀(hν/(kTₑ))log_γ2
: Equal to log₁₀(γ²) = log₁₀(RydbergZ²/(kTₑ))
Rydberg
is the "infinite mass unit of energy" and Tₑ
is the temperature of free electrons (for our purposes, we assume that free electrons are in thermal equilibrium with ions and neutral species).
Notes
van Hoof et al. (2014) computed the associated data table with a non-relativistic approach, which is invalid at very high temperatures. They conclude (from comparisons with a different paper) that their "results should be accurate up to electron temperatures of roughly 100 MK". This is more than adequate for stellar atmospheres. In van Hoof et al. (2015), they find that relativistic effects introduce a ∼0.75% at 100MK, for Z = 1 (when Z > 1, the change is smaller).
This function currently uses linear interpolation. However, van Hoof et al. (2014) provides an implementation of a third-order Lagrange scheme, which "reaches a relative precision better than 1.5e-4 everywhere." The C and Fortran implementations of this scheme can be found here, and are copyrighted by a BSD-style license.
Earlier variants of this function used less-accurate data from section 5.1 of Kurucz (1970) that extended over a smaller interval of data. That table was originally derived from a figure in Karsas and Latter (1961) and it's now used for testing purposes.
Korg.ContinuumAbsorption.hydrogenic_ff_absorption
— Methodhydrogenic_ff_absorption(ν, T, Z, ni, ne)
computes the free-free linear absorption coefficient for a hydrogenic species
The naming convention for free-free absorption is counter-intuitive. A free-free interaction is named as though the species interacting with the free electron had one more bound electron (in other words it's named as though the free-electron and ion were bound together). In practice, this means that ni
should refer to:
- the number density of H II if computing the H I free-free absorption
- the number density of He III if computing the He II free-free absorption
- the number density of Li IV if computing the Li III free-free absorption
Arguments
Z::Integer
: the charge of the ion. For example, this is 1 for ionized H.ni
: the number density of the ion species in cm⁻³.ne
: the number density of free electrons.ν
: frequency in HzT
: temperature in K
Note
This approach was adopted from equation 5.8 from section 5.1 of Kurucz (1970). Comparison against equation 5.18b of Rybicki & Lightman (2004), reveals that the equation in Kurucz (1970) omits the dependence on ρ. According to Rybicki & Lightman (2004) the free-free absorption coefficient (corrected for stimulated emission) is:
α = coef * Z² * ne * ni * (1 - exp(-hplanck*ν/(kboltz*T))) * g_ff / (sqrt(T) * ν³)
Note that the g_ff is the free-free gaunt factor and coef is ∼3.7e8 (a more exact coefficient can be computed from eqn 5.18a).
With this in mind, equation 5.8 of Kurucz (1970) should actually read
κ = ne * n(H II) * F_ν(T) * (1 - exp(-hplanck*ν/(kboltz*T))) / ρ
where Fν(T) = coef * Z² * gff / (sqrt(T) * ν³).
See gaunt_ff_vanHoof
for details about where our gaunt factor data comes from. For simplicity, we enforce temperature and ν bounds constraints that don't include all of the available gaunt factor data. In practice, this should never be a concern for stellar spectroscopy.
Korg.ContinuumAbsorption.metal_bf_absorption!
— Methodmetal_bf_absorption!(α, νs, T, number_densities)
Adds to α the contributions of bf metal opacities. Uses precomputed tables from TOPBase for Li, Be, B, C, N, O, F, Ne, Na, Mg, Al, Si, S, Ar, and Ca. Uses tables from NORAD for Fe. For these elements, tables have been precomputed for the neutral and singly ionized species assuming and LTE distribution of energy levels. See Korg/data/metal_bf_cross-sections/
for the scripts which generate the tables.
Cross sections were computed for 100 K < T < 100,000 K and frequencies corresponding to 500 Å < λ < 30,000 Å. Outside of either of those ranges, flat extrapolation is used (i.e. the extreme value is returned).
Korg.ContinuumAbsorption.positive_ion_ff_absorption!
— Methodpositive_ion_ff_absorption!(α_out::Vector{Real}, ν::Real, T::Real, number_densities::Dict,
ne::Real, departure_coefficients=Peach1970.departure_coefficients)
Computes the linear absorption coefficient (in cm⁻¹) for all free-free interactions involving positively charged atomic species. Uses the provided departure coefficients when they are available, and the uncorrected hydrogenic approximation when they are not.
Arguments
ν
: frequency in HzT
: temperature in Knumber_densities
is aDict
mapping eachSpecies
to its number density.ne
: the number density of free electrons.departure_coefficients
(optional, defaults toContinuumAbsorption.Peach1970.departure_coefficients
: a dictionary mapping species to the departure coefficients for the ff process it participates in (e.g.species"C II"
maps to the C III ff departure coefficients–see note below). Departure coefficients should be callables taking temperature and (photon energy / RydbergH / Zeff^2).
A free-free interaction is named as though the species interacting with the free electron had one more bound electron (in other words it's named as though the free-electron and ion were bound together). To make matters more confusing, some sources present the free-free linear absorption coefficient's formula in terms of the number density of the species in the interaction's name (by including the Saha equation in the formula).
To compute the absorption for a given free-free interaction, this function accesses elements from the dictionary arguments that are associated with the species that participates in the interaction. For example:
- Si I ff absorption: uses the number density of Si II (
number_densities[species"Si II"]
) and checksdeparture_coefficients[species"Si II"]
for departure coefficients. - Si II ff absorption: uses the number density of Si III (
number_densities[species"Si III"]
) and checksdeparture_coefficients[species"Si III"]
for departure coefficients.
Korg.ContinuumAbsorption.rayleigh
— Methodrayleigh(λs, nH_I, nHe_II, nH2)
Absorption coefficient from Rayleigh scattering by neutral H, He, and H2. Formulations for H and He are via Colgan+ 2016. Formulation for H2 from Dalgarno and Williams 1962.
The Dalgarno and Williams H2 is applicable redward of 1300 Å. Since Rayleigh scattering breaks down when the particle size to wavelength ratio gets large, we that all frequencies passed to this function be equivalent to 1300 Å or greater.
The formulations for H is adapted from Lee 2005, which states that it is applicable redward of Lyman alpha. See Colgan 2016 for details on He.
Korg.ContinuumAbsorption.simple_hydrogen_bf_cross_section
— Methodsimple_hydrogen_bf_cross_section(n::Integer, ν::Real)
Calculate the H I bf cross section in megabarns using a very simple approximation. See, for example, Kurucz 1970 equation 5.5 (though see note below). This implementation is used to extrapolate the cross-section past the ionization energy of an unperturbed hydrogen atom, as is required to take level dissolution into account with the MHD formalism.
Equation 5.5 of Kurucz had a typo in it. In the numerator of the fraction that is multiplied by the entire polynomial, Z² should be Z⁴. This was discovered during comparisons with data from the Opacity Project, and can be confirmed by looking at eqn 5.6 of Kurucz (it uses Z⁴ instead of Z²) or by comparison against equation 10.54 of Rybicki & Lightman.
Korg.ContinuumAbsorption.total_continuum_absorption
— Methodtotal_continuum_absorption(νs, T, nₑ, number_densities, partition_funcs; error_oobounds)
The total continuum linear absoprtion coefficient, α, at many frequencies, ν.
Arguments
νs
are frequencies in HzT
is temperature in Knₑ
is the electron number density in cm^-3number_densities
is aDict
mapping eachSpecies
to its number densitypartition_funcs
is aDict
mapping eachSpecies
to its partition function (e.g.Korg.partition_funcs
)error_oobounds::Bool
specifies the behavior of most continnum absorption sources when passed frequencies or temperature values that are out of bounds for their implementation. Whenfalse
(the default), those absorption sources are ignored at those values. Otherwise, an error is thrown.
For efficiency reasons, νs
must be sorted. While this function technically supports any sorted AbstractVector
, it is most effient when passed an AbstractRange
.
Korg.ContinuumAbsorption.Stancil1994
— ModuleThis module contains transcriptions of the tables from Stancil 1994, who calculated ff and bf absorption coefficients for H₂⁺ and He₂⁺.
It also contains pre-constructed interpolation objects.
Korg.ContinuumAbsorption.Peach1970.departure_coefficients
— MethodPeach1970.departure_coefficients()
This module contains interpolators of the tabulated ff departure coeffients from Peach+ 1970, which we use to correct the hydrogenic ff absorption coefficient for H I ff, C I ff, Si I ff, and Mg I ff. It contains a dictionary (returned by departure_coefficients()
), which maps Species
to interpolator objects. Crucially, the dictionary is indexed by the species which actually participates in the interaction, not the one after which the interaction is named.
Outside the regime in which Peach 1970 provides data, the interpolators return 0, falling back to the hydreogenic approximation.
The species for which we use corrections are the same species which get corrected in MARCS/Turbospectrum (see Table 1 of Gustafsson+ 2008). The choices seem are largely motivated by which species have departure terms at normal stellar atmosphere conditions and which species are most abundant in the sun. For C II ff, we include only the contribution from the ¹S parent term, even though (in contrast to other speices) information is available for the ³Pᵒ term as well.
The free-free absorption coefficient (including stimulated emission) is given by:
$\alpha_{ m ff} = \alpha_{ m hydrogenic, ff}(\nu, T, n_i, n_e; Z) (1 + D(T, \sigma))$,
where
- $\alpha_{ m hydrogenic, ff}(\nu, T, n_i, n_e; Z)$ should include the correction for stimulated emission.
- $n_i$ is the number density fo the ion species that participates in the interation, not the species the interaction is named after.
- $n_e$ is the number density of free electrons.
- $D(T, \sigma)$ is specified as the
departure
arg, and is expected to interpolate over the tabulated values specified in Table III of Peach (1970). - σ denotes the energy of the photon in units of RydbergH*Zeff²
It might not be immediately obvious how the above equation relates to the equations presented in Peach (1970). Peach describes the calculation for $k_ u^F$, the free-free absorption coefficient (uncorrected for stimulated emission) per particle of the species that the interaction is named after. In other words, he computes:
$k_ u^F = \alpha_{ m ff}/n_{i-1} \left(1 - e^\frac{-h\nu}{k T}\right)^{-1}$,
where $n_{i-1}$ is the number density of the species that the interaction is named after. $k_ u^F$ can directly be computed, under LTE, from just $u$, $T$, and $n_{i-1}$ (the Saha Equation relates $\alpha_{ m ff}$'s dependence on $n_e$ and $n_i$ to $n_{i-1}$ and $T$. Gray (2005) follows a similar convention when describing free-free absorption.
The tabulated data in this module was taken from Peach 1970 using OCR software, and may contain mis-read values, although they produce reasonable behavior and there are no obvious problems.
Korg.RadiativeTransfer.calculate_rays
— Methodcalculate_rays(μ_surface_grid, spatial_coord, spherical)
Arguments
μ_surface_grid
: the μ values at the surface of the star corresponding to the rays along which the optical depth and intensity will be calculated.spatial_coord
: a physical distance coordinate. This is radius for a spherical atmosphere, and height above the photosphere for a plane-parallel atmosphere.
Returns
A vector of pairs (s, ds/dz)
, where s
is the distance along the ray and z is the model atmosphere spatial coordinate.
Korg.RadiativeTransfer.compute_F_flux_only_expint
— Methodcompute_F_flux_only_expint(τ, S)
Compute the astrophysical flux, F, by linearly interpolating the source function, S
across optical depths τ
. Handle the integral over μ analytically using E₂.
Korg.RadiativeTransfer.compute_I_bezier!
— Methodcompute_I_bezier!(I, τ, S)
Given τ and S along a ray (at a particular wavelength), compute the intensity at the end of the ray (the surface of the star). This uses the method from de la Cruz Rodríguez and Piskunov 2013.
Korg.RadiativeTransfer.compute_I_linear!
— Methodcompute_I_linear!(I, τ, S)
Compute exactly the solution to the transfer integral obtained be linearly interpolating the source function, S
across optical depths τ
, without approximating the factor of exp(-τ).
This breaks the integral into the sum of integrals of the form $\int (m\tau + b) \exp(-\tau)$ d\tau$ , which is equal to $ -\exp(-\tau) (m*\tau + b + m)$.
Korg.RadiativeTransfer.compute_I_linear_flux_only
— Methodcompute_I_linear_flux_only(τ, S)
Returns the intensity at the end of the ray (the surface of the star) given τ and S along the ray. Uses the same numerical method as compute_I_linear!
, but doesn't retain the intensity at each layer.
Korg.RadiativeTransfer.compute_tau_bezier!
— Methodcompute_tau_bezier(τ, s, α)
Compute optical depth (write to τ) along a ray with coordinate s and absorption coefficient α. This is the method proposed in de la Cruz Rodríguez and Piskunov 2013, but the
Korg.RadiativeTransfer.expint_transfer_integral_core
— Methodexpint_transfer_integral_core(τ, m, b)
The exact solution to $\int (m\tau + b) E_2(\tau)$ d\tau$. The exponential integral function, expint, captures the integral over the disk of the star to get the emergent astrophysical flux. You can verify it by substituting the variable of integration in the exponential integal, t, with mu=1/t.
Korg.RadiativeTransfer.exponential_integral_2
— Methodexponential_integral_2(x)
Approximate second order exponential integral, E_2(x). This stitches together several series expansions to get an approximation which is accurate within 1% for all x
Korg.RadiativeTransfer.fritsch_butland_C
— Methodfritsch_butland_C(x, y)
Given a set of x and y values, compute the bezier control points using the method of Fritch & Butland 1984, as suggested in de la Cruz Rodríguez and Piskunov 2013.
Used in compute_I_bezier!
.
Korg.RadiativeTransfer.generate_mu_grid
— Methodgenerate_mu_grid(n_points)
generate_mu_grid(μ_values)
Used by both radiative transfer schemes to compute quadrature over μ. Returns (μ_grid, μ_weights)
. If an integer is passed, generate a grid of Gauss-Legendre quadrature points of corresponding size. Otherwise, use the provided μ values.
Korg.RadiativeTransfer.radiative_transfer
— Methodradiative_transfer(atm, α, S, n_μ_points; kwargs...)
radiative_transfer(α, S, spatial_coord, n_μ_points, spherical; kwargs...)
Arguments (note that some of these are only required for one method):
atm
: the model atmosphere.α
: a matrix (atmospheric layers × wavelengths) containing the absorption coefficientS
: the source function as a matrix of the same shape. rescale the total absorption to match the model atmosphere. This value should be calculated by Korg.μ_points
: the number of quadrature points to use when integrating over I_surface(μ) to obtain the astrophysical flux, or, alternatively, the specific μ values to use (in which case the integral is done with the trapezoid rule).spatial_coord
: a physical distance coordinate (radius for a spherical atmosphere, height above the photosphere for a plane-parallel atmosphere).spherical
: whether the atmosphere is spherical or plane-parallel.
Keyword Arguments:
include_inward_rays
(default:false
): if true, light propagating into the star (negative μs) is included. If false, only those which are needed to seed the intensity at the bottom of the atmosphere are included.τ_scheme
(default: "anchored"): how to compute the optical depth. Options are "linear" and "bezier" (not recommended).I_scheme
(default: "linearfluxonly"): how to compute the intensity. Options are "linear", "linearfluxonly", and "bezier". "linearfluxonly" is the fastest, but does not return the intensity values anywhere except at the top of the atmosphere. "linear" performs an equivalent calculation, but stores the intensity at every layer. "bezier" is not recommended.
Korg.RadiativeTransfer.radiative_transfer_core
— MethodCompute the intensity along a single ray.
n.b. this function has an additional Ischeme ("linearfluxonlyexpint") that radiative_transfer will automatically switch to when appropriate.