Autodiff

Project.toml

@@ -8,10 +8,11 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 Interpolations = "0.13"
 NLsolve = "4.5"
-SpecialFunctions = "1.3"
 StaticArrays = "1"
+SpecialFunctions = "1.4"
 julia = "1.0"

src/Korg.jl

@@ -1,5 +1,6 @@
 module Korg
-    export synthesize, constant_R_LSF, read_line_list, read_model_atmosphere
+    export synthesize, constant_R_LSF, rectify, air_to_vacuum, vacuum_to_air, read_line_list, 
+            read_model_atmosphere
 
     _data_dir = joinpath(@__DIR__, "../data") 
 
@@ -18,6 +19,6 @@ module Korg
 
     include("continuum_opacity/continuum_opacity.jl") #Define continuum opacity functions.
     include("synthesize.jl")                          #solve radiative transfer equation
-    include("LSF.jl")                                 #functions to apply LSF
+    include("utils.jl")                               #functions to apply LSF, vac<->air wls, etc
 
 end # module

src/continuum_opacity/continuum_opacity.jl

@@ -28,5 +28,5 @@ I adopted the formula described in section 5.12 of Kurucz (1970) and the equatio
 scattering subsection of Gray (2005); the actual coefficient value comes from the latter. It turns
 out that the coefficient in Kurucz (1970) has a typo (it's a factor of 10 too large).
 """
-electron_scattering(nₑ::F, ρ::F) where {F<:AbstractFloat} = 0.6648e-24*nₑ/ρ
+electron_scattering(nₑ::F, ρ::F) where {F<:Real} = 0.6648e-24*nₑ/ρ
 end

src/continuum_opacity/hydrogenic_bf_ff.jl

@@ -10,20 +10,19 @@ Calculates the LTE number density (in cm^-3) of a hydrogenic species at a given
 
 # Arguments
 - `n::Integer`: The quantum number of the state
-- `nsdens_div_partition::Flt`: The total (including all states) number density of the current 
+- `nsdens_div_partition`: The total (including all states) number density of the current 
 ionization species (in cm^-3) divided by the species's partition function.
-- `ion_energy::Flt`: The minimum energy (in eV) required to ionize the species (from the ground 
+- `ion_energy`: The minimum energy (in eV) required to ionize the species (from the ground 
 state). This can be estimated as Z²*ion_energy of hydrogen or Z²*Rydberg_H.
-- `T::Flt`: Temperature
+- `T`: Temperature
 
 # Note
 This assumes that all hydrogenic species have a statistical weight of gₙ = 2*n².
 
 This was taken from equation (5.4) of Kurucz (1970) (although this comes directly from the 
 boltzmann equation).
 """
-function ndens_state_hydrogenic(n::Integer, nsdens_div_partition::Flt, T::Flt,
-                                ion_energy::Flt) where {Flt<:AbstractFloat}
+function ndens_state_hydrogenic(n::Integer, nsdens_div_partition::Real, T::Real, ion_energy::Real)
     n2 = n*n
     g_n = 2.0*n2
     energy_level = ion_energy - ion_energy/n2
@@ -46,7 +45,7 @@ const _bf_σ_coef = SA[(0.9916,  2.719e13, -2.268e30),
 
 # this is helper function is the most likely part of the calculation to change.
 # This uses double precision to be safe about the polynomial coefficients.
-function _hydrogenic_bf_cross_section(Z::Integer, n::Integer, ν::Float64, ion_freq::Float64)
+function _hydrogenic_bf_cross_section(Z::Integer, n::Integer, ν::Real, ion_freq::Real)
     # this implements equation 5.5 from Kurucz (1970)
     # - Z is the atomic number
     # - n is the energy level (remember, they start from n=1)
@@ -88,20 +87,18 @@ where α(n=n') is the absorption coefficient for the bound-free atomic absorptio
 # Arguments
 - `Z::Integer`: Z is the atomic number of the ion (1 for HI)
 - `nmin::Integer`: The lowest energy level (principle quantum number) included in the calculation
-- `nsdens_div_partition::AbstractFloat` is the number density of the current species divided by the
+- `nsdens_div_partition` is the number density of the current species divided by the
    partition function.
-- `ν::Flt`: frequency in Hz
-- `ρ::Flt`: mass density in g/cm³
-- `T::Flt`: temperature in K
-- `ion_energy::AbstractFloat`: the ionization energy from the ground state (in eV).
+- `ν`: frequency in Hz
+- `ρ`: mass density in g/cm³
+- `T`: temperature in K
+- `ion_energy`: the ionization energy from the ground state (in eV).
 
 # Notes
 This implements equation (5.6) from Kurucz (1970). I think ρ was simply omitted from that equation.
 """
-function _hydrogenic_bf_high_n_opacity(Z::Integer, nmin::Integer,
-                                       nsdens_div_partition::AbstractFloat,
-                                       ν::AbstractFloat, ρ::AbstractFloat, T::AbstractFloat,
-                                       ion_energy::AbstractFloat)
+function _hydrogenic_bf_high_n_opacity(Z::Integer, nmin::Integer, nsdens_div_partition::Real, 
+                                       ν::Real, ρ::Real, T::Real, ion_energy::Real)
 
     # this function corresponds to the evaluation of a integral. We subdivide the solution into two
     # parts: (i) consts and (ii) integral. The solution is the product of both parts
@@ -149,12 +146,12 @@ integral.
 
 # Arguments
 - `Z::Integer`: Z is the atomic number of the species (e.g. 1 for H I or 2 for He II)
-- `nsdens_div_partition::Flt` is the total number density of the species divided by the species's
+- `nsdens_div_partition` is the total number density of the species divided by the species's
    partition function.
-- `ν::Flt`: frequency in Hz
-- `ρ::Flt`: mass density in g/cm³
-- `T::Flt`: temperature in K
-- `ion_energy::AbstractFloat`: the ionization energy from the ground state (in eV). This can be 
+- `ν`: frequency in Hz
+- `ρ`: mass density in g/cm³
+- `T`: temperature in K
+- `ion_energy`: the ionization energy from the ground state (in eV). This can be 
    estimated as Z²*Rydberg_H (Rydberg_H is the ionization energy of Hydrogen)
 - `nmax_explicit_sum::Integer`: The highest energy level whose opacity contribution is included in
    the explicit sum. The contributions from higher levels are included in the integral.
@@ -164,9 +161,9 @@ integral.
 # Notes
 This follows the approach described in section 5.1 of Kurucz (1970).
 """
-function hydrogenic_bf_opacity(Z::Integer, nsdens_div_partition::Flt, ν::Flt, ρ::Flt, T::Flt,
-                               ion_energy::Flt, nmax_explicit_sum::Integer,
-                               integrate_high_n::Bool = true) where {Flt<:AbstractFloat}
+function hydrogenic_bf_opacity(Z::Integer, nsdens_div_partition::Real, ν::Real, ρ::Real, T::Real, 
+                               ion_energy::Real, nmax_explicit_sum::Integer, 
+                               integrate_high_n::Bool=true)
     ionization_freq = _eVtoHz(ion_energy)
 
     # first, directly sum individual the opacity contributions from H I atoms at each of the lowest
@@ -221,8 +218,8 @@ computes the thermally averaged free-free gaunt factor by interpolating the tabl
 section 5.1 of Kurucz (1970). The table was derived from a figure in Karsas and Latter (1961).
 
 # Arguments
-- `log_u::F`: Equal to log₁₀(u) = log₁₀(h*ν/(k*T))
-- `log_γ2::F`: Equal to log₁₀(γ²) = log₁₀(RydbergH*Z²/(k*T))
+- `log_u`: Equal to log₁₀(u) = log₁₀(h*ν/(k*T))
+- `log_γ2`: Equal to log₁₀(γ²) = log₁₀(RydbergH*Z²/(k*T))
 
 # Note
 There is some ambiguity over whether we should replace RydbergH*Z² in the definition of γ² with the
@@ -250,11 +247,11 @@ means that `ni` should refer to:
 
 # Arguments
 - `Z::Integer`: the charge of the ion. For example, this is 1 for ionized H.
-- `ni::Flt`: the number density of the ion species in cm⁻³.
-- `ne::Flt`: the number density of free electrons.
-- `ν::Flt`: frequency in Hz
-- `ρ::Flt`: mass density in g/cm³
-- `T::Flt`: temperature in K
+- `ni`: the number density of the ion species in cm⁻³.
+- `ne`: the number density of free electrons.
+- `ν`: frequency in Hz
+- `ρ`: mass density in g/cm³
+- `T`: temperature in K
 
 # Note
 This approach was adopted from equation 5.8 from section 5.1 of Kurucz (1970). Comparison against
@@ -269,10 +266,9 @@ 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² * g_ff / (sqrt(T) * ν³).
 """
-function hydrogenic_ff_opacity(Z::Integer, ni::Flt, ne::Flt, ν::Flt, ρ::Flt,
-                               T::Flt) where {Flt<:AbstractFloat}
+function hydrogenic_ff_opacity(Z::Integer, ni::Real, ne::Real, ν::Real, ρ::Real, T::Real)
     β = 1.0/(kboltz_eV * T)
-    Z2 = convert(Flt, Z*Z)
+    Z2 = Z*Z
 
     hν_div_kT = hplanck_eV * ν * β
     log_u = log10(hν_div_kT)

src/continuum_opacity/opacity_H.jl

@@ -51,7 +51,7 @@ Compute the H⁻ bound-free cross-section, which has units of cm^2 per H⁻ part
 
 The cross-section does not include a correction for stimulated emission.
 """
-function _Hminus_bf_cross_section(ν::AbstractFloat)
+function _Hminus_bf_cross_section(ν::Real)
     λ = c_cgs*1e8/ν # in Angstroms
     # we need to somehow factor out this bounds checking
     if !(2250 <= λ <= 15000.0)
@@ -113,8 +113,8 @@ data from Wishart (1979). While Gray (2005) claims that the polynomial fits the
 precision for 2250 Å ≤ λ ≤ 15000 Å, in practice we find that it fits the data to better than 0.25%
 precision. Wishart (1979) expects the tabulated data to have better than 1% percent accuracy.
 """
-function Hminus_bf(nH_I_div_partition::Flt, ne::Flt, ν::Flt, ρ::Flt, T::Flt,
-                   ion_energy_H⁻::Flt = 0.7552) where {Flt<:AbstractFloat}
+function Hminus_bf(nH_I_div_partition::Real, ne::Real, ν::Real, ρ::Real, T::Real, 
+                   ion_energy_H⁻::Real = 0.7552)
     αbf_H⁻ = _Hminus_bf_cross_section(ν) # does not include contributions from stimulated emission
     stimulated_emission_correction = (1 - exp(-hplanck_cgs*ν/(kboltz_cgs*T)))
     n_H⁻ = _ndens_Hminus(nH_I_div_partition, ne, T, _H⁻_ion_energy)
@@ -162,8 +162,7 @@ n(H I, n = 1) over the temperature range where the polynomial is valid.
 We also considered the polynomial fit in Section 5.3 from Kurucz (1970). Unfortunately, it seems
 wrong (it gives lots of negative numbers).
 """
-function Hminus_ff(nH_I_div_partition::Flt, ne::Flt, ν::Flt, ρ::Flt,
-                   T::Flt) where {Flt<:AbstractFloat}
+function Hminus_ff(nH_I_div_partition::Real, ne::Real, ν::Real, ρ::Real, T::Real)
     λ = c_cgs*1e8/ν # in Angstroms
     # we need to somehow factor out this bounds checking
     if !(2604 <= λ <= 113918.0)
@@ -207,12 +206,15 @@ This uses polynomial fits from Gray (2005) that were derived from data tabulated
 [Bates (1952)](https://ui.adsabs.harvard.edu/abs/1952MNRAS.112...40B/abstract).
 
 # Arguments
-- `nH_I_div_partition::Float64`: the total number density of H I divided by its partition 
+- `nH_I_div_partition`: the total number density of H I divided by its partition 
    function.
-- `nH_II::Float64`: the number density of H II (not of H₂⁺).
-- `ν::Float64`: frequency in Hz
-- `ρ::Float64`: mass density in g/cm³
-- `T::Float64`: temperature in K
+- `nH_II`: the number density of H II (not of H₂⁺).
+- `ν`: frequency in Hz
+- `ρ`: mass density in g/cm³
+- `T`: temperature in K
+
+While the formal type signature requires that these be `Real`, `Float32`s (or `Float32`-derived 
+types) may introduce numerical instability.
 
 # Notes
 This follows equation 8.15 of Gray (2005), which involves 2 polynomials that were fit to data
@@ -247,8 +249,7 @@ times larger than the max λ the polynomials are fit against). He suggests that
 probably correct "to well within one part in ten even at the lower temperatures and [lower
 wavelengths]."
 """
-function H2plus_bf_and_ff(nH_I_div_partition::Float64, nH_II::Float64, ν::Float64, ρ::Float64,
-                          T::Float64)
+function H2plus_bf_and_ff(nH_I_div_partition::Real, nH_II::Real, ν::Real, ρ::Real, T::Real)
     λ = c_cgs*1e8/ν # in ångstroms
     if !(3846.15 <= λ <= 25000.0) # the lower limit is set to include 1.e5/26 Å
         throw(DomainError(λ, "The wavelength must lie in the interval [3847 Å, 25000 Å]"))

src/continuum_opacity/opacity_He.jl

@@ -19,7 +19,7 @@ He_II_ff(nHe_III, ne, ν, ρ, T) = hydrogenic_ff_opacity(2, nHe_III, ne, ν, ρ,
 
 # Compute the number density of atoms in different He I states
 # taken from section 5.5 of Kurucz (1970)
-function ndens_state_He_I(n::Integer, nsdens_div_partition::Flt, T::Flt) where {Flt<:AbstractFloat}
+function ndens_state_He_I(n::Integer, nsdens_div_partition::Real, T::Real)
     g_n, energy_level = if n == 1
         (1.0, 0.0)
     elseif n == 2
@@ -84,8 +84,7 @@ approximations for 5.06e3 Å ≤ λ ≤ 1e4 Å "are expected to be well below 10
 
 An alternative approach using a fit to older data is provided in section 5.7 of Kurucz (1970).
 """
-function Heminus_ff(nHe_I_div_partition::Flt, ne::Flt, ν::Flt, ρ::Flt,
-                    T::Flt) where {Flt<:AbstractFloat}
+function Heminus_ff(nHe_I_div_partition::Real, ne::Real, ν::Real, ρ::Real, T::Real)
 
     λ = c_cgs * 1.0e8 / ν # Å
     if !(5063.0 <= λ <= 151878.0)

src/fit.jl

@@ -0,0 +1,24 @@
+using Random
+using Optim
+
+function best_fit_params(atm, linelist, wls, f_obs; ivar=ones(length(f_obs)),
+                         free_parameters=["metallicity", "vmic"], fixed_values::Dict=Dict(),
+                         initial_values::Dict=Dict(), wl_deflation_factor=0.01)
+
+    function chi2(params, wl_mask)
+        f_synth = synthesize(atm, linelist, wls[wl_mask]; metallicity=params[1],
+                             abundances=Dict(["Ca"=>params[2]]), verbose=false).flux
+        f_synth = Korg.constant_R_LSF(f_synth, wls[wl_mask], 11500.0)
+        sum((f_obs[wl_mask] .- f_synth).^2 .* ivar[wl_mask])
+    end
+
+    mask = rand(length(wls)) .< wl_deflation_factor
+
+    result = optimize(p->chi2(p, mask), [0.0, 0.0], NewtonTrustRegion(),
+        Optim.Options(show_trace=true, extended_trace=true, x_tol=1e-1);
+        autodiff=:forward)
+
+    optimize(p->chi2(p, ones(Bool, length(wls))), result.minimizer, NewtonTrustRegion(),
+        Optim.Options(show_trace=true, extended_trace=true, x_tol=1e-3);
+        autodiff=:forward)
+end

src/line_opacity.jl

@@ -24,7 +24,16 @@ otherwise specified
 """
 function line_absorption(linelist, λs, temp, nₑ, n_densities::Dict, partition_fns::Dict, ξ 
                          ; α_cntm=nothing, cutoff_threshold=1e-3, window_size=20.0*1e-8)
-    α_lines = zeros(length(λs))
+    if length(linelist) == 0
+        return zeros(length(λs))
+    end
+
+    #type shenanigans to allow autodiff to do its thing
+    α_type = typeof(promote(linelist[1].wl, λs[1], temp, nₑ, n_densities["H_I"], ξ, cutoff_threshold
+                           )[1])
+    α_lines = zeros(α_type, length(λs))
+
+
     #lb and ub are the indices to the upper and lower wavelengths in the "window", i.e. the shortest
     #and longest wavelengths which feel the effect of each line 
     lb = 1
@@ -38,7 +47,7 @@ function line_absorption(linelist, λs, temp, nₑ, n_densities::Dict, partition
         #get all damping params from line list.  There may be better sources for this.
         Γ = line.gamma_rad 
         if !ismolecule(line.species) 
-            Γ += (nₑ*scaled_stark(line.gamma_stark, temp) + 
+            Γ += (nₑ*scaled_stark(line.gamma_stark, temp) +
                   (n_densities["H_I"] + 0.42n_densities["He_I"])*scaled_vdW(line.vdW, mass, temp))
         end
 
@@ -59,7 +68,8 @@ function line_absorption(linelist, λs, temp, nₑ, n_densities::Dict, partition
 
         if !isnothing(α_cntm)
             α_crit = α_cntm(line.wl) * cutoff_threshold
-            window_size = max(4Δλ_D, sqrt(line_amplitude * Δλ_L / α_crit))
+            Δλ_crit = sqrt(line_amplitude * Δλ_L / α_crit) #where α from lorentz component == α_crit
+            window_size = max(4Δλ_D, Δλ_crit)
         end
         lb, ub = move_bounds(λs, lb, ub, line.wl, window_size)
         if lb==ub
@@ -75,14 +85,13 @@ end
 
 "the stark broadening gamma scaled acording to its temperature dependence"
 scaled_stark(γstark, T; T₀=10_000) = γstark * (T/T₀)^(1/6)
-
 
 """
 the vdW broadening gamma scaled acording to its temperature dependence, using either simple scaling 
 or ABO
 """
-scaled_vdW(vdW::AbstractFloat, m, T, T₀=10_000) = vdW * (T/T₀)^0.3
-function scaled_vdW(vdW::Tuple{F, F}, m, T) where F <: AbstractFloat
+scaled_vdW(vdW::Real, m, T, T₀=10_000) = vdW * (T/T₀)^0.3
+function scaled_vdW(vdW::Tuple{F, F}, m, T) where F <: Real
     v₀ = 1e6 #σ is given at 10_000 m/s = 10^6 cm/s
     σ = vdW[1]
     α = vdW[2]
@@ -123,7 +132,7 @@ end
 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`.
 """
-function sigma_line(λ) where F <: AbstractFloat
+function sigma_line(λ::Real)
     #work in cgs
     e  = electron_charge_cgs
     mₑ = electron_mass_cgs
@@ -138,11 +147,11 @@ end
 A normalized voigt profile centered on λ₀ with doppler width `invΔλ_D` = 1/Δλ_D and lorentz width 
 `Δλ_L` evaluated at `λ` (cm).  Note that this returns values in units of cm^-1.
 """
-function line_profile(λ₀, invΔλ_D::F, Δλ_L, line_amplitude::F, λ::F) where F <: AbstractFloat
+function line_profile(λ₀, invΔλ_D::F, Δλ_L, line_amplitude::F, λ::F) where F <: Real
     _line_profile(λ₀, invΔλ_D, Δλ_L*invΔλ_D/(4π), line_amplitude*invΔλ_D/sqrt(π), λ)
 end
 function _line_profile(λ₀, invΔλ_D::F, Δλ_L_invΔλ_D_div_4π::F, amplitude_invΔλ_D_div_sqrt_π::F, λ::F
-                      ) where F <:  AbstractFloat
+                      ) where F <: Real
     voigt(Δλ_L_invΔλ_D_div_4π, abs(λ-λ₀) * invΔλ_D) * amplitude_invΔλ_D_div_sqrt_π
 end