interpolated metal-poor atmospheres with negative optical depths

[segasai]

Hi,

I don't know if it's worth reporting, but I looked at korg failures over large grid that are not 'chemistry convergence' issues.
And this
thing stand out:

ERROR: LoadError: DomainError with -7.986061283680538e-6:
log was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
    @ Base.Math ./special/log.jl:301
  [3] log
    @ ./special/log.jl:267 [inlined]
  [4] _broadcast_getindex_evalf
    @ ./broadcast.jl:709 [inlined]
  [5] _broadcast_getindex
    @ ./broadcast.jl:682 [inlined]
  [6] getindex
    @ ./broadcast.jl:636 [inlined]
  [7] macro expansion
    @ ./broadcast.jl:1004 [inlined]
  [8] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [9] copyto!
    @ ./broadcast.jl:1003 [inlined]
 [10] copyto!
    @ ./broadcast.jl:956 [inlined]
 [11] copy
    @ ./broadcast.jl:928 [inlined]
 [12] materialize
    @ ./broadcast.jl:903 [inlined]
 [13] radiative_transfer(α::Matrix{Float64}, S::Matrix{Float64}, spatial_coord::Vector{Float64}, μ_points::Int64, spherical::Bool; include_inward_rays::Bool, α_ref::Vector{Float64}, τ_ref::Vector{Float64}, I_scheme::String, τ_scheme::String)
    @ Korg.RadiativeTransfer ~/.julia/packages/Korg/6TfTZ/src/RadiativeTransfer/RadiativeTransfer.jl:110
 [14] radiative_transfer
    @ ~/.julia/packages/Korg/6TfTZ/src/RadiativeTransfer/RadiativeTransfer.jl:73 [inlined]
 [15] #radiative_transfer#1
    @ ~/.julia/packages/Korg/6TfTZ/src/RadiativeTransfer/RadiativeTransfer.jl:61 [inlined]
 [16] synthesize(::Korg.PlanarAtmosphere{Float64, Float64, Float64, Float64, Float64}, ::Vector{Korg.Line{Float64, Float64, Float64, Float64, Float64, Float64}}, ::Vector{Float64}, ::Int64, ::Vararg{Int64}; vmic::Int64, line_buffer::Int64, cntm_step::Float64, air_wavelengths::Bool, hydrogen_lines::Bool, use_MHD_for_hydrogen_lines::Bool, hydrogen_line_window_size::Int64, mu_values::Int64, line_cutoff_threshold::Float64, electron_number_density_warn_threshold::Float64, electron_number_density_warn_min_value::Float64, return_cntm::Bool, I_scheme::String, tau_scheme::String, ionization_energies::Dict{UInt8, Vector{Float64}}, partition_funcs::Dict{Korg.Species, Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, log_equilibrium_constants::Dict{Korg.Species, Union{Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Korg.var"#logK#57"{Korg.Species, Dict{Korg.Species, Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}}}}, molecular_cross_sections::Vector{Any}, use_chemical_equilibrium_from::Nothing, verbose::Bool)
    @ Korg ~/.julia/packages/Korg/6TfTZ/src/synthesize.jl:224
 [17] top-level scope
 

This is caused by this log

Korg.jl/src/RadiativeTransfer/RadiativeTransfer.jl

Line 110 in 7516c46

log_τ_ref = log.(τ_ref)

The reproducer (requires the linelist that i previously shared)

using Korg
using Serialization

lines=open(deserialize, "vald_all_3000_10000.julia.bin")
teff=(4118.689742870629)
logg=(5.291771519929171)
X_H=[12. , 10.93 , -2.84274731, -2.51274731, -1.19274731,
4.49725269, 3.88725269, 5.35591968, 0.66725269, 4.53591968,
2.27725269, 4.22591968, 2.47725269, 4.20591968, 1.46725269,
3.83591968, 1.60725269, 2.87591968, 1.18725269, 3.00591968,
-0.72274731, 1.59591968, 0.10725269, 1.74725269, 1.49725269,
3.55725269, 1.02725269, 2.33725269, 0.31725269, 0.70725269,
-1.01274731, -0.31274731, -1.60274731, -0.56274731, -1.33274731,
-0.64274731, -1.29274731, -0.97274731, -1.68274731, -1.31274731,
-2.47274731, -1.97274731, -8.89274731, -2.05274731, -2.77274731,
-2.23274731, -2.95274731, -2.12274731, -2.29274731, -1.89274731,
-2.89274731, -1.70274731, -2.38274731, -1.65274731, -2.82274731,
-1.72274731, -2.76274731, -2.19274731, -3.31274731, -2.44274731,
-8.89274731, -2.89274731, -3.37274731, -2.78274731, -3.61274731,
-2.75274731, -3.38274731, -2.96274731, -3.89274731, -2.81274731,
-3.83274731, -3.01274731, -4.06274731, -2.78274731, -3.66274731,
-2.64274731, -2.51274731, -2.25274731, -2.88274731, -2.76274731,
-2.99274731, -1.89274731, -3.24274731, -8.89274731, -8.89274731,
-8.89274731, -8.89274731, -8.89274731, -8.89274731, -3.83274731,
                                                                     -8.89274731, -4.41274731]
marcs = Korg.interpolate_marcs((teff),
                                   (logg),
                                   X_H,
                                   clamp_abundances=true)

minlam=3400
maxlam=5000
spec = Korg.synthesize(marcs,
                       lines,X_H,minlam,maxlam,vmic=1,
                        line_buffer=100,
                           hydrogen_line_window_size=400)

[ajwheeler]

Thanks for reporting this. It turns out to be a problem arising in the metal-poor atmosphere interpolator, where for a subset of stellar parameters, some of the optical depths at 5000 Angstrom go negative.

Here's a quick plot showing where the problem arises.

code to generate

teffs = 7875 : -50 : 2800
loggs = 5.1 : 0.05 : 5.25
m_Hs = -5 : 0.5 : -2.5

@time all_atms = interpolate_marcs.(teffs, loggs', reshape(m_Hs, 1, 1, :), 0.4, 0);

f(a) = sum(Korg.get_tau_5000s(a) .< 0)
all_negative = f.(all_atms);

for (i, m_H) in enumerate(m_Hs)
    ps = teffs .=> loggs'
    scatter(first.(ps)[:], last.(ps)[:] .+ (m_H+2.5)*1e-2, c=all_negative[:, :, i][:], s=2, vmin=0, vmax=10)
end
colorbar(label="number of layers with τ_5000 < 0")
gca().invert_xaxis()
gca().invert_yaxis()
xlabel("Teff [K]")
ylabel("log(g) + ([m/H]+2)/100")

Ideally this would never happen, but for now I will make sure to at least throw an error at atmosphere interpolation. I will leave this issue open until this is fully fixed.

[segasai]

Thanks for looking into this.
I don't know the insides of the atmosphere's interpolation in korg, but if it's done in tau, presumably one can switch the interpolation to log(tau) that would naturally enforce the non-negativity.

[ajwheeler]

I will look into doing this. Currently the reference optical depths are directly interpolated (some of these other are log or otherwise transformed), but I will need to check what this does to accuracy in other parts of parameter space.