Merge pull request #313 from ajwheeler/correct_tau_5

include line contributions to `alpha_5000`

src/Korg.jl

@@ -7,11 +7,11 @@ module Korg
     include("atomic_data.jl")              # symbols and atomic weights
     include("isotopic_data.jl")            # abundances and nuclear spins
     include("species.jl")                  # types for chemical formulae and species
+    include("read_statmech_quantities.jl") # approximate Us, Ks, chis
     include("linelist.jl")                 # parse linelists, define Line type
     include("line_absorption.jl")          # opacity, line profile, voigt function
     include("hydrogen_line_absorption.jl") # hydrogen lines get special treatment
     include("autodiffable_conv.jl")        # wrap DSP.conv to be autodiffable
-    include("read_statmech_quantities.jl") # approximate Us, Ks, chis
     include("statmech.jl")                 # statistical mechanics, molecular equilibrium
     include("atmosphere.jl")               # parse model atmospheres
     include("RadiativeTransfer/RadiativeTransfer.jl") # radiative transfer formal solution

src/line_absorption.jl

@@ -2,7 +2,7 @@ using SpecialFunctions: gamma
 using ProgressMeter: @showprogress
 
 """
-    line_absorption(linelist, λs, temp, nₑ, n_densities, partition_fns, ξ
+    line_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

src/linelist.jl

@@ -46,7 +46,7 @@ struct Line{F1, F2, F3, F4, F5, F6}
         if wl >= 1
             wl *= 1e-8 #convert Å to cm
         end
-        if ismissing(gamma_stark) || isnan(gamma_stark) || ismissing(vdW) || isnan(vdW)
+        if ismissing(gamma_stark) || isnan(gamma_stark) || ismissing(vdW) || (!(vdW isa Tuple) && isnan(vdW))
             gamma_stark_approx, vdW_approx = approximate_gammas(wl, species, E_lower)
             if ismissing(gamma_stark) || isnan(gamma_stark)
                 gamma_stark = gamma_stark_approx
@@ -61,7 +61,7 @@ struct Line{F1, F2, F3, F4, F5, F6}
 
         # if vdW is a tuple, assume it's (σ, α) from ABO theory
         # if it's a float, there are four possibilities
-        if !ismissing(vdW) && !isnan(vdW) && !(vdW isa Tuple) #F6 will not be defined if vdW is missing
+        if !ismissing(vdW) && !(vdW isa Tuple) !isnan(vdW) #F6 will not be defined if vdW is missing
             if vdW < 0 
                 vdW = 10^vdW  # if vdW is negative, assume it's log(γ_vdW) 
             elseif vdW == 0
@@ -83,6 +83,9 @@ function Base.show(io::IO, ::MIME"text/plain", line::Line)
     print(io, " ", round(line.wl*1e8, digits=6), " Å (log gf = ", round(line.log_gf, digits=2) ,")")
 end
 
+# make it broadcast like a scalar
+Base.broadcastable(l::Line) = Ref(l)
+
 """
     approximate_radiative_gamma(wl, log_gf)
 
@@ -630,4 +633,35 @@ function get_GES_linelist()
             tentotheOrMissing.(Float64.(read(f["gamma_stark"]))),
             read(f["vdW"]))
     end
-end
+end
+
+"""
+    _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.
+"""
+function _load_alpha_5000_linelist(path=joinpath(_data_dir, "linelists", "alpha_5000", "alpha_5000_lines.csv"))
+    csv = CSV.File(path)
+    map(csv) do row
+        vdW = if ',' in row.vdW
+            σ, α = split(row.vdW[2:end-1], ',')
+            (parse(Float64, σ), parse(Float64, α))
+        else
+            parse(Float64, row.vdW)
+        end
+        Line(row.wl, row.log_gf, Species(row.species), row.E_lower, row.gamma_rad, row.gamma_stark, vdW)
+    end
+end
+
+"""
+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.
+
+See also [`_load_alpha_5000_linelist`](@ref).
+"""
+const _alpha_5000_default_linelist = _load_alpha_5000_linelist()

src/molecular_cross_sections.jl

@@ -1,9 +1,9 @@
-using Interpolations: GriddedInterpolation, interpolate!, Gridded, NoInterp, Linear
+using Interpolations: interpolate!, Gridded, Linear, extrapolate
 using HDF5
 
 struct MolecularCrossSection
-    wls # just for debugging
-    itp :: GriddedInterpolation
+    wls 
+    itp
     species :: Korg.Species
 end
 
@@ -13,13 +13,13 @@ end
 Precompute the molecular absorption cross section for a given linelist and set of wavelengths. The 
 `MolecularCrossSection` object can be passed to [`synthesize`](@ref) and potentially speed up the 
 calculation significantly.  At present, Korg only supports precomputed cross-sections created by 
-this function.
+this function, though they can be saved and loaded using [`save_molecular_cross_section`](@ref) and 
+[`read_molecular_cross_section`](@ref).
 
 # Arguments
 - `linelist`: A vector of `Line` objects representing the molecular linelist.  These must be of the 
    same species.
-- `wls`: A vector of wavelength ranges (in Å) at which to precompute the cross section.  *These must
-  match the wavelengths used for any subsequent synthesis exactly*.
+- `wls`: A vector of wavelength ranges (in Å) at which to precompute the cross section.
 
 # Keyword Arguments
 - `cutoff_alpha` (default: 1e-30): The value of the single-line absorption coefficient (in cm^-1) at 
@@ -63,21 +63,21 @@ function MolecularCrossSection(linelist, wls; cutoff_alpha=1e-32,
     end
 
     species = all_specs[1]
-    itp = interpolate!((vmic_vals, log_temp_vals, 1:size(α, 3)),
-                       α .* cutoff_alpha, 
-                       (Gridded(Linear()), Gridded(Linear()), NoInterp()))
+    itp = extrapolate(interpolate!((vmic_vals, log_temp_vals, vcat(wls...)*1e-8),
+                      α .* cutoff_alpha, 
+                      (Gridded(Linear()), Gridded(Linear()), Gridded(Linear()))), 0.0)
     MolecularCrossSection(wls, itp, species)
 end
 
 
 """
-    interpolate_molecular_cross_sections!(α, molecular_cross_sections, Ts, vmic, number_densities)
+    interpolate_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`](@ref) for more information.
 """
-function interpolate_molecular_cross_sections!(α::Matrix{R}, molecular_cross_sections, Ts, vmic,
-                                               number_densities) where R <: Real
+function interpolate_molecular_cross_sections!(α::AbstractArray{R}, molecular_cross_sections, λs,
+                                               Ts, vmic, number_densities) where R <: Real
     if length(molecular_cross_sections) == 0
         return
     end
@@ -86,7 +86,7 @@ function interpolate_molecular_cross_sections!(α::Matrix{R}, molecular_cross_se
         for i in 1:size(α, 1)
             # this is an inefficient order in which to write to α, but doing the interpolations this
             # way uses less memory.
-            view(α, i, :) .+= sigma.itp(vmic, log10(Ts[i]), 1:size(α, 2)) * number_densities[sigma.species][i]
+            view(α, i, :) .+= sigma.itp.(vmic, log10(Ts[i]), λs) * number_densities[sigma.species][i]
         end
     end
     ;
@@ -105,9 +105,9 @@ function save_molecular_cross_section(filename, cross_section)
 
     HDF5.h5open(filename, "w") do file
         HDF5.write(file, "wls", [(l[begin], step(l), l[end]) for l in wls])
-        HDF5.write(file, "vmic_vals", collect(itp.knots[1]))
-        HDF5.write(file, "T_vals", collect(itp.knots[2]))
-        HDF5.write(file, "vals", itp.coefs)
+        HDF5.write(file, "vmic_vals", collect(itp.itp.knots[1]))
+        HDF5.write(file, "T_vals", collect(itp.itp.knots[2]))
+        HDF5.write(file, "vals", itp.itp.coefs)
         HDF5.write(file, "species", string(species))
     end
 end
@@ -125,13 +125,15 @@ function read_molecular_cross_section(filename)
             start:step:stop
         end
         vmic_vals = HDF5.read(file, "vmic_vals")
-        T_vals = HDF5.read(file, "T_vals")
-        alpha_vals = HDF5.read(file, "vals")
+        logT_vals = HDF5.read(file, "T_vals")
+        # doesn't actually matter that it's mmaped because it's passed to interpolate!
+        alpha_vals = HDF5.readmmap(file["vals"]) 
         species = Species(HDF5.read(file, "species"))
 
-        itp = interpolate!((vmic_vals, T_vals, 1:sum(length.(wls))),
-                           alpha_vals,
-                           (Gridded(Linear()), Gridded(Linear()), NoInterp()))
+        itp = extrapolate(interpolate!((vmic_vals, logT_vals, vcat(wls...)*1e8),
+                                        alpha_vals,
+                                        (Gridded(Linear()), Gridded(Linear()), Gridded(Linear()))),
+                          0.0)
 
         MolecularCrossSection(wls, itp, species)
     end

src/prune_linelist.jl

@@ -21,18 +21,21 @@ equivalent width.
    order is much faster, but sorting by strength is useful for visualizing the strongest lines.
 - `verbose=true`: If `true`, a progress bar will be displayed while measuring the EWs.
 All other kwargs are passed to internal calls to [`synthesize`](@ref).
+- `max_distance=0.0`, how far from `wls` lines can be before they are excluded from the returned 
+   list. 
 
 !!! caution
     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=false` if you are pruning the linelist to speedup synthesis.  Note that Korg will 
+    `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`](@ref)).
 
 See also [`merge_close_lines`](@ref) if you are using this for plotting.
 """
 function prune_linelist(atm, linelist, A_X, wls...; 
-                        threshold=0.1, sort_by_EW=true, verbose=true, synthesis_kwargs...)
+                        threshold=0.1, sort_by_EW=true, verbose=true, max_distance=0.0,
+                        synthesis_kwargs...)
     # linelist will be sorted after call to synthesize
     sol = synthesize(atm, linelist, A_X, wls...; synthesis_kwargs...)
     cntm_sol = synthesize(atm, [], A_X, wls...; synthesis_kwargs...) 
@@ -59,7 +62,7 @@ function prune_linelist(atm, linelist, A_X, wls...;
     strong_lines = eltype(linelist)[]
     for line in linelist
         line_center = line.wl * 1e8
-        if !any(λs[begin] < line_center < λs[end] for λs in wl_ranges)
+        if !any((λs[begin]-max_distance) < line_center < (λs[end]+max_distance) for λs in wl_ranges)
             continue
         end
 

src/synthesize.jl

@@ -88,7 +88,8 @@ solution = synthesize(atm, linelist, A_X, 5000, 5100)
 - `bezier_radiative_transfer` (default: false): Use the radiative transfer scheme.  This is for 
    testing purposes only.
 - `molecular_cross_sections` (default: `[]`): A vector of precomputed molecular cross-sections. See 
-   [`MolecularCrossSection`](@ref) for how to generate these.
+   [`MolecularCrossSection`](@ref) 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.
 - `use_chemical_equilibrium_from` (default: `nothing`): Takes another solution returned by 
    `synthesize`. 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`,
@@ -152,7 +153,12 @@ function synthesize(atm::ModelAtmosphere, linelist, A_X::AbstractVector{<:Real},
     # frequencies at which to calculate the continuum, as a single vector
     sorted_cntmνs = c_cgs ./ reverse(cntmλs) 
 
-    # sort linelist and remove lines far from the synthesis region 
+    # sort linelist and remove lines far from the synthesis region
+    # first just the ones needed for α5 (fall back to default if they aren't provided)
+    if !bezier_radiative_transfer
+        linelist5 = get_alpha_5000_linelist(linelist)
+    end
+    # now the ones for the synthesis
     linelist = filter_linelist(linelist, wl_ranges, line_buffer)
 
     if length(A_X) != MAX_ATOMIC_NUMBER || (A_X[1] != 12)
@@ -204,6 +210,15 @@ function synthesize(atm::ModelAtmosphere, linelist, A_X::AbstractVector{<:Real},
     #vector of continuum-absorption interpolators
     α_cntm = last.(triples) 
 
+    # line contributions to α5
+    if ! bezier_radiative_transfer
+        α_cntm_5 = [_->a for a in copy(α5)] # lambda per layer
+        line_absorption!(view(α5, :, 1), linelist5, [5e-5:1:5e-5], get_temps(atm), nₑs, number_densities,
+                         partition_funcs, vmic*1e5, α_cntm_5; cutoff_threshold=line_cutoff_threshold)
+        interpolate_molecular_cross_sections!(view(α5, :, 1), molecular_cross_sections, [5e-5],
+                                              get_temps(atm), vmic, number_densities)
+    end
+
     source_fn = blackbody.((l->l.temp).(atm.layers), all_λs')
     cntm = nothing
     if return_cntm
@@ -225,10 +240,9 @@ function synthesize(atm::ModelAtmosphere, linelist, A_X::AbstractVector{<:Real},
         end
     end
 
-    line_absorption!(α, linelist, wl_ranges, [layer.temp for layer in atm.layers], nₑs,
-        number_densities, partition_funcs, vmic*1e5, α_cntm, cutoff_threshold=line_cutoff_threshold;
-        verbose=verbose)
-    interpolate_molecular_cross_sections!(α, molecular_cross_sections, get_temps(atm), vmic, number_densities)
+    line_absorption!(α, linelist, wl_ranges, get_temps(atm), nₑs, number_densities, partition_funcs,
+                     vmic*1e5, α_cntm; cutoff_threshold=line_cutoff_threshold, verbose=verbose)
+    interpolate_molecular_cross_sections!(α, molecular_cross_sections, all_λs, get_temps(atm), vmic, number_densities)
 
     flux, intensity = if bezier_radiative_transfer
         RadiativeTransfer.BezierTransfer.radiative_transfer(atm, α, source_fn, n_mu_points)
@@ -266,7 +280,7 @@ construct_wavelength_ranges(wls::AbstractRange) = [wls]
 
 Return a new linelist containing only lines within the provided wavelength ranges.
 """
-function filter_linelist(linelist, wl_ranges, line_buffer)
+function filter_linelist(linelist, wl_ranges, line_buffer; warn_empty=true)
     # this could be made faster by using a binary search (e.g. searchsortedfirst/last)
 
     #sort the lines if necessary
@@ -285,12 +299,46 @@ function filter_linelist(linelist, wl_ranges, line_buffer)
         return false
     end
 
-    if nlines_before != 0 && length(linelist) == 0
+    if nlines_before != 0 && length(linelist) == 0 && warn_empty
         @warn "The provided linelist was not empty, but none of the lines were within the provided wavelength range."
     end
     linelist
 end
 
+"""
+    get_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`](@ref))
+"""
+function get_alpha_5000_linelist(linelist)
+    # start by getting the lines in the provided linelist which effect the synthesis at 5000 Å
+    # use a 21 Å line buffer, which 1 Å bigger than the coverage of the fallback linelist.
+    linelist5 = filter_linelist(linelist, [5e-5:1:5e-5], 21e-8; warn_empty=false)
+    # if there aren't any, use the built-in one
+    if length(linelist5) == 0
+        _alpha_5000_default_linelist
+    # if there are some, but they don't actually cross over 5000 Å, use the built-in one where they 
+    # aren't present
+    elseif linelist5[1].wl   > 5e-5
+        ll = filter(_alpha_5000_default_linelist) do line
+            line.wl < linelist5[1].wl
+        end
+        [ll ; linelist5]
+    elseif linelist5[end].wl < 5e-5
+        ll = filter(_alpha_5000_default_linelist) do line
+            line.wl > linelist5[end].wl
+        end
+        [linelist5 ; ll]
+    else # if the built-in lines span 5000 Å, use them
+        linelist5
+    end
+end
+
 """
     format_A_X(default_metals_H, default_alpha_H, abundances; kwargs... )
 

test/autodiff.jl

@@ -19,19 +19,24 @@
 
     @testset "autodiff just one abundance"  begin
         atm = read_model_atmosphere("data/sun.mod")
-        linelist = [Korg.Line(5000e-8, 0.0, Korg.species"C I", 0.0)]
+        linelist = [Korg.Line(6000e-8, 0.0, Korg.species"C I", 0.0)]
+        # If this line is super weak (or removed), the test will fail due to numerics in the 
+        # abundances and continuum opacities. It used to be a fake line at 5000, but the fact that 
+        # that's the reference wavelength caused instability in the finie differnces calculation.
         function f(A_C)
             A_X = format_A_X(Dict("C"=>A_C))
-            synthesize(atm, linelist, A_X, 5000, 5000).flux[1]
+            synthesize(atm, linelist, A_X, 6000, 6000).flux[1]
         end
         @test FiniteDiff.finite_difference_derivative(f, 0.0) ≈ ForwardDiff.derivative(f, 0.0) rtol=1e-4
     end
 
     @testset "line params" begin
         atm = read_model_atmosphere("data/sun.mod")
         function f(loggf)
-            linelist = [Korg.Line(5000e-8, loggf, Korg.species"Na I", 0.0)]
-            synthesize(atm, linelist, format_A_X(), 5000, 5000).flux[1]
+            linelist = [Korg.Line(6000e-8, loggf, Korg.species"Na I", 0.0)]
+            # This used to be a fake line at 5000, but the fact that that's the reference wavelength
+            # caused instability in the finie differnces calculation.
+            synthesize(atm, linelist, format_A_X(), 6000, 6000).flux[1]
         end
         @test FiniteDiff.finite_difference_derivative(f, 0.0) ≈ ForwardDiff.derivative(f, 0.0) rtol=1e-4
     end

test/molecular_cross_sections.jl

@@ -27,7 +27,7 @@
     filename = tempname()
     Korg.save_molecular_cross_section(filename, table)
     deserialized_table = Korg.read_molecular_cross_section(filename)
-    @test table.itp == deserialized_table.itp
+    @test table.itp.itp.coefs == deserialized_table.itp.itp.coefs
     @test all(table.wls .== deserialized_table.wls)
     @test table.species == deserialized_table.species
 end

test/qfactors.jl

@@ -35,7 +35,7 @@
     msk[1:100] .= false
 
     Q = Korg.Qfactor(synth_flux, synth_wls, apowls, LSF_model; obs_mask=msk)
-    @test Q ≈ 812.2698584824295
+    @test Q ≈ 810.9267657701387 rtol=1e-4
     SNR = flux ./ obs_err
     RMS_SNR = sqrt(mean(SNR[msk].^2))
     Q_prec = Korg.RV_prec_from_Q(Q, RMS_SNR, count(msk))