Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

throw an error if interpolated atmosphere has negative optical depths #381

Merged
merged 1 commit into from
Dec 12, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 13 additions & 1 deletion src/atmosphere.jl
Original file line number Diff line number Diff line change
Expand Up @@ -219,6 +219,13 @@ function _get_cool_dwarfs_atm_itp()
_cool_dwarfs_atm_itp
end

struct AtmosphereInterpolationError <: Exception
msg::String
end
function Base.showerror(io::IO, e::AtmosphereInterpolationError)
print(io, "Chemical equilibrium failed: ", e.msg)
end

"""
interpolate_marcs(Teff, logg, A_X; kwargs...)
interpolate_marcs(Teff, logg, m_H=0, alpha_m=0, C_m=0; kwargs...)
Expand Down Expand Up @@ -292,7 +299,7 @@ function interpolate_marcs(Teff, logg, m_H=0, alpha_m=0, C_m=0; spherical=logg <
archives=(_sdss_marcs_atmospheres, _get_cool_dwarfs_atm_itp(),
_low_Z_marcs_atmospheres))
# cool dwarfs
if Teff <= 4000 && logg >= 3.5 && m_H >= -2.5 && resampled_cubic_for_cool_dwarfs
atm = if Teff <= 4000 && logg >= 3.5 && m_H >= -2.5 && resampled_cubic_for_cool_dwarfs
itp, nlayers = archives[2]
atm_quants = itp(1:nlayers, 1:5, Teff, logg, m_H, alpha_m, C_m)
PlanarAtmosphere(PlanarAtmosphereLayer.(atm_quants[:, 4],
Expand Down Expand Up @@ -337,6 +344,11 @@ function interpolate_marcs(Teff, logg, m_H=0, alpha_m=0, C_m=0; spherical=logg <
exp.(atm_quants[nanmask, 3])))
end
end

if any(get_tau_5000s(atm) .< 0)
throw(AtmosphereInterpolationError("Interpolated atmosphere has negative optical depths and is not reliable. See https://github.com/ajwheeler/Korg.jl/issues/378 for details."))
end
atm
end
# handle the case where Teff, logg, and [m/H] are integers. As long as not all (interpolated) params
# are passed in as integers, there's no problem.
Expand Down
2 changes: 2 additions & 0 deletions test/atmosphere.jl
Original file line number Diff line number Diff line change
Expand Up @@ -131,6 +131,8 @@
@test assert_atmospheres_close(atm2, atm3; tol=1e-10)

@test_throws ArgumentError interpolate_marcs(5000, 4.5, -3, 0)

@test_throws Korg.AtmosphereInterpolationError interpolate_marcs(4100, 5.25, -4, 0.4, 0)
end

@testset "integer arguments" begin
Expand Down
38 changes: 19 additions & 19 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,26 +5,26 @@ using Korg, Test, Logging, HDF5, ForwardDiff, FiniteDiff, Aqua
# tools for testing: assert_allclose and assert_allclose_grid
include("utilities.jl")

#include("wavelengths.jl")
#include("molecular_cross_sections.jl")
#include("cubic_splines.jl")
#include("transfer.jl")
#include("species.jl")
#include("interval.jl")
#include("continuum_absorption.jl") # test this after the "Interval" testset
#include("partition_funcs.jl")
#include("statmech.jl")
#include("linelist.jl")
include("wavelengths.jl")
include("molecular_cross_sections.jl")
include("cubic_splines.jl")
include("transfer.jl")
include("species.jl")
include("interval.jl")
include("continuum_absorption.jl") # test this after the "Interval" testset
include("partition_funcs.jl")
include("statmech.jl")
include("linelist.jl")
include("fit.jl") # slow
#include("autodiff.jl") # slow
#include("autodiffable_conv.jl")
#include("atmosphere.jl") # slow
#include("abundances.jl")
#include("synthesize.jl")
#include("prune_linelist.jl")
#include("utils.jl")
#include("line_absorption.jl")
#include("qfactors.jl")
include("autodiff.jl") # slow
include("autodiffable_conv.jl")
include("atmosphere.jl") # slow
include("abundances.jl")
include("synthesize.jl")
include("prune_linelist.jl")
include("utils.jl")
include("line_absorption.jl")
include("qfactors.jl")
@testset "Aqua automated checks" begin
# see https://github.com/JuliaTesting/Aqua.jl/issues/77 for why I'm doing it this way,
# (basically the default bahavior is different and this is what we want to avoid errors in
Expand Down
Loading