diff --git a/docs/plots.jl b/docs/plots.jl index 4736322..53fb052 100644 --- a/docs/plots.jl +++ b/docs/plots.jl @@ -1,5 +1,6 @@ using Plots, LaTeXStrings -import DustExtinction: ccm89_ca, ccm89_cb, od94_ca, od94_cb, cal00_invum, ccm89_invum, vcg04_invum, gcc09_invum, f99_invum, f04_invum, f19_invum, m14_invum, FM90, G16 + +import DustExtinction: ccm89_ca, ccm89_cb, od94_ca, od94_cb, cal00_invum, ccm89_invum, vcg04_invum, gcc09_invum, f99_invum, f04_invum, f19_invum, m14_invum, FM90, G16, P92 dir = joinpath(@__DIR__, "src", "assets") @@ -155,6 +156,43 @@ xlabel!(L"\mu m ^{-1}") ylabel!(L"E(\lambda - V)/E(B - V)") savefig(joinpath(dir, "FM90_plot.svg")) +#-------------------------------------------------------------------------------- +# P92 + +# plotting x in micrometers +# wave numbers generated in mu-1 +w = exp10.(range(-3, 3, step = 0.001)) + +# wave generated in angstrom +x = 1e4 ./ w +# size was modified to accomodate the tally table +plot(size = (700, 500), xaxis = :log10, yaxis = :log10, xticks = ([0.001, 0.01, 0.1, 1, 10, 100, 1000]), yticks = ([0.001, 0.01, 0.1, 1, 10])) + +m1 = P92().(x) +plot!(1 ./w, m1, label = "total") + +m2 = P92(FUV_amp = 0, NUV_amp = 0, SIL1_amp = 0, SIL2_amp = 0, FIR_amp = 0).(x) +plot!(1 ./w, m2, label = "BKG only") + +m3 = P92(NUV_amp=0.0, SIL1_amp=0.0, SIL2_amp=0.0, FIR_amp=0.0).(x) +plot!(1 ./w, m3, label = "BKG+FUV only") + +m4 = P92(FUV_amp=0.0, SIL1_amp=0.0, SIL2_amp=0.0, FIR_amp=0.0).(x) +plot!(1 ./w, m4, label = "BKG+NUV only") + +m5 = P92(FUV_amp = 0, NUV_amp = 0, SIL2_amp = 0).(x) +plot!(1 ./w, m5, label = "BKG+FIR+SIL1 only") + +m6 = P92(FUV_amp=0., NUV_amp=0.0, SIL1_amp=0.0).(x) +plot!(1 ./w, m6, label = "BKG+FIR+SIL2 only") + +m7 = P92(FUV_amp=0., NUV_amp=0.0, SIL1_amp=0.0, SIL2_amp=0.0).(x) +plot!(1 ./w, m7, label = "BKG+FIR only") + +xlabel!(L"x\ \left[\mu m\right]") +ylabel!(L"A(X)/A(V)") +savefig(joinpath(dir, "P92_plot.svg")) + #-------------------------------------------------------------------------------- # G16 diff --git a/docs/src/assets/P92.bib b/docs/src/assets/P92.bib new file mode 100644 index 0000000..958d982 --- /dev/null +++ b/docs/src/assets/P92.bib @@ -0,0 +1,13 @@ +@ARTICLE{1992ApJ...395..130P, + author = {{Pei}, Yichuan C.}, + title = "{Interstellar Dust from the Milky Way to the Magellanic Clouds}", + journal = {\apj}, + keywords = {Cosmic Dust, Intergalactic Media, Interstellar Extinction, Interstellar Matter, Magellanic Clouds, Milky Way Galaxy, Chemical Evolution, Far Ultraviolet Radiation, Kramers-Kronig Formula, Astrophysics, GALAXIES: INTERGALACTIC MEDIUM, GALAXIES: INTERSTELLAR MATTER, GALAXIES: MAGELLANIC CLOUDS, ISM: DUST, EXTINCTION}, + year = 1992, + month = aug, + volume = {395}, + pages = {130}, + doi = {10.1086/171637}, + adsurl = {https://ui.adsabs.harvard.edu/abs/1992ApJ...395..130P}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} diff --git a/docs/src/assets/P92_plot.svg b/docs/src/assets/P92_plot.svg new file mode 100644 index 0000000..a81317f --- /dev/null +++ b/docs/src/assets/P92_plot.svg @@ -0,0 +1,4344 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/color_laws.md b/docs/src/color_laws.md index 1378bfe..7620225 100644 --- a/docs/src/color_laws.md +++ b/docs/src/color_laws.md @@ -228,6 +228,14 @@ DustExtinction.bounds FM90 ``` +#### Pei (1992) + +![](assets/P92_plot.svg) + +```@docs +P92 +``` + ### Mixture Extinction Laws #### Gordon et al. (2016) diff --git a/docs/src/index.md b/docs/src/index.md index f2ae54f..515f9f6 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -42,6 +42,7 @@ There are various citations relevant to this work. Please be considerate when us | [`SFD98Map`](@ref) | [Schlegel, Finkbeiner and Davis (1998)](https://ui.adsabs.harvard.edu/abs/1998ApJ...500..525S) | [download](assets/sfd98.bib) | | [`F99`](@ref) | [Fitzpatrick (1999)](https://ui.adsabs.harvard.edu/abs/1999PASP..111...63F) | [download](assets/f99.bib) | | [`F04`](@ref) | [Fitzpatrick (2004)](https://ui.adsabs.harvard.edu/abs/2004ASPC..309...33F) | [download](assets/f04.bib) | +| [`P92`](@ref) | [Pei (1992)](https://ui.adsabs.harvard.edu/abs/1992ApJ...395..130P) | [download](assets/P92.bib) | | [`F19`](@ref) | [Fitzpatrick (2019)](https://ui.adsabs.harvard.edu/abs/2019ApJ...886..108F) | [download](assets/f19.bib) | | [`M14`](@ref) | [Maiz Apellaniz et al. (2014)](https://ui.adsabs.harvard.edu/abs/2014A%26A...564A..63M) | [download](assets/m14.bib) | diff --git a/src/DustExtinction.jl b/src/DustExtinction.jl index f504989..518a3ec 100644 --- a/src/DustExtinction.jl +++ b/src/DustExtinction.jl @@ -19,6 +19,7 @@ export redden, M14, # Fittable laws FM90, + P92, # Mixture laws G16, G03_SMCBar, @@ -36,7 +37,7 @@ The abstract super-type for dust extinction laws. See the extended help (`??Dust ## Interface Here's how to make a new extinction law, called `MyLaw` -* Create your struct. We strongly recommend using `Parameters.jl` to facilitate creating keyword argument constructors if your model is parameterized, which allows convenient usage with [`redden`](@ref) and [`deredden`](@ref). +* Create your struct. We strongly recommend using `Parameters.jl` to facilitate creating keyword argument constructors if your model is parameterized, which allows convenient usage with [`redden`](@ref) and [`deredden`](@ref). ```julia struct MyLaw <: DustExtinction.ExtinctionLaw end ``` @@ -152,7 +153,8 @@ include("mixture_laws.jl") # at which point adding `(l::ExtinctionLaw)(wave::Quantity)` is possible, until then # using this code-gen does the trick but requires manually editing # instead of providing support for all <: ExtinctionLaw -for law in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, G16, G03_SMCBar, F99, F04, F19, M14] + +for law in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, G16, G03_SMCBar, F99, F04, F19, M14, P92] (l::law)(wavelength::Quantity) = l(ustrip(u"Å", wavelength)) * u"mag" end diff --git a/src/fittable_laws.jl b/src/fittable_laws.jl index bde19d0..afde40c 100644 --- a/src/fittable_laws.jl +++ b/src/fittable_laws.jl @@ -44,7 +44,7 @@ the overall normalization, possibly changing the expected behavior of reddening ## References [Fitzpatrick & Massa (1990)](https://ui.adsabs.harvard.edu/abs/1990ApJS...72..163F) """ -@with_kw struct FM90{T<:Number} <: ExtinctionLaw @deftype T +@with_kw struct FM90{T<:Number} <: ExtinctionLaw @deftype T c1 = 0.10 c2 = 0.70 c3 = 3.23 @@ -77,3 +77,172 @@ function (law::FM90)(wave::T) where T return exvebv end + +""" + P92(BKG_amp=218.57, + BKG_lambda=0.047, + BKG_b=90.0, + BKG_n=2.0, + FUV_amp=18.54, + FUV_lambda=0.07, + FUV_b=4.0, + FUV_n=6.5, + NUV_amp=0.0596, + NUV_lambda=0.22, + NUV_b=-1.95, + NUV_n=2.0, + SIL1_amp=0.0026, + SIL1_lambda=9.7, + SIL1_b=-1.95, + SIL1_n=2.0, + SIL2_amp=0.0026, + SIL2_lambda=18.0, + SIL2_b=-1.8, + SIL2_n=2.0, + FIR_amp=0.0159, + FIR_lambda=25.0, + FIR_b=0.0, + FIR_n=2.0) + +Pei (1992) generative extinction model applicable from the extreme UV to far-IR. + +## Parameters + +Background Terms +* `BKG_amp` - amplitude +* `BKG_lambda` - central wavelength +* `BKG_b` - b coefficient +* `BKG_n` - n coefficient + +Far-Ultraviolet Terms +* `FUV_amp` - amplitude +* `FUV_lambda` - central wavelength +* `FUV_b` - b coefficent +* `FUV_n` - n coefficient + +Near-Ultraviolet (2175 Å) Terms +* `NUV_amp` - amplitude +* `NUV_lambda` - central wavelength +* `NUV_b` - b coefficent +* `NUV_n` - n coefficient + +1st Silicate Feature (~10 micron) Terms +* `SIL1_amp` - amplitude +* `SIL1_lambda` - central wavelength +* `SIL1_b` - b coefficent +* `SIL1_n` - n coefficient + +2nd Silicate Feature (~18 micron) Terms +* `SIL2_amp` - amplitude +* `SIL2_lambda` - central wavelength +* `SIL2_b` - b coefficient +* `SIL2_n` - n coefficient + +Far-Infrared Terms +* `FIR_amp` - amplitude +* `FIR_lambda` - central wavelength +* `FIR_b` - b coefficent +* `FIR_n` - n coefficient + +If `λ` is a `Unitful.Quantity` it will be automatically converted to Å and the returned value will be `UnitfulAstro.mag`. + +## Examples +```jldoctest +julia> model = P92(); + +julia> model(1500) +2.396291891812002 + +julia> P92(FUV_b = 2.0).([1000, 2000, 3000]) +3-element Array{Float64,1}: + 3.8390886792306187 + 2.7304534614548697 + 1.806181164464396 + +``` + +## Default Parameter Values + +|Term |lambda|A|b|n| +|:---:|:---:|:---:|:---:|:---:| +|BKG |0.047 |218.57142857142858 |90 |2 | +|FUV |0.07 |18.545454545454547 |4.0 |6.5| +|NUV |0.22 |0.05961038961038961 |-1.95|2.0| +|SIL1 |9.7 |0.0026493506493506496|-1.95|2.0| +|SIL2 |18.0 |0.0026493506493506496|-1.8 |2.0| +|FIR |25.0 |0.015896103896103898 |0.0 |2.0| + + +## References +[Pei (1992)](https://ui.adsabs.harvard.edu/abs/1992ApJ...395..130P) +""" +@with_kw struct P92{T<:Number} <: ExtinctionLaw @deftype T + BKG_amp = 165.0 * (1 / 3.08 + 1) + BKG_lambda = 0.047 + BKG_b = 90.0 + BKG_n = 2.0 + FUV_amp = 14.0 * (1 / 3.08 + 1) + FUV_lambda = 0.07 + FUV_b = 4.0 + FUV_n = 6.5 + NUV_amp = 0.045 * (1 / 3.08 + 1) + NUV_lambda = 0.22 + NUV_b = -1.95 + NUV_n = 2.0 + SIL1_amp = 0.002 * (1 / 3.08 + 1) + SIL1_lambda = 9.7 + SIL1_b = -1.95 + SIL1_n = 2.0 + SIL2_amp = 0.002 * (1 / 3.08 + 1) + SIL2_lambda = 18.0 + SIL2_b = -1.80 + SIL2_n = 2.0 + FIR_amp = 0.012 * (1 / 3.08 + 1) + FIR_lambda = 25.0 + FIR_b = 0.00 + FIR_n = 2.0 + @assert BKG_amp ≥ 0 "`BKG_amp` must be ≥ 0, got $BKG_amp" + @assert FUV_amp ≥ 0 "`FUV_amp` must be ≥ 0, got $FUV_amp" + @assert 0.06 ≤ FUV_lambda ≤ 0.08 "`FUV_lambda` must be in between [0.06, 0.08], got $FUV_lambda" + @assert NUV_amp ≥ 0 "`NUV_amp` must be ≥ 0, got $NUV_amp" + @assert 0.20 ≤ NUV_lambda ≤ 0.24 "`NUV_lambda` must be in between [0.20, 0.24], got $NUV_lambda" + @assert SIL1_amp ≥ 0 "`SIL1_amp` must be ≥ 0, got $SIL1_amp" + @assert 7 ≤ SIL1_lambda ≤ 13 "`SIL1_lambda` must be in between [7, 13], got $SIL1_lambda" + @assert SIL2_amp ≥ 0 "`SIL2_amp` must be ≥ 0, got $SIL2_amp" + @assert 15 ≤ SIL2_lambda ≤ 21 "`SIL2_lambda` must be in between [15, 21], got $SIL2_lambda" + @assert FIR_amp ≥ 0 "`FIR_amp` must be ≥ 0, got $FIR_amp" + @assert 20 ≤ FIR_lambda ≤ 30 "`FIR_lambda` must be in between [20, 30], got $FIR_lambda" +end + +P92(BKG_amp, BKG_lambda, BKG_b, BKG_n, FUV_amp, FUV_lambda, FUV_b, FUV_n, + NUV_amp, NUV_lambda, NUV_b, NUV_n, SIL1_amp, SIL1_lambda, SIL1_b, + SIL1_n, SIL2_amp, SIL2_lambda, SIL2_b, SIL2_n, FIR_amp, FIR_lambda, + FIR_b, FIR_n) = + P92(promote(BKG_amp, BKG_lambda, BKG_b, BKG_n, FUV_amp, FUV_lambda, FUV_b, FUV_n, + NUV_amp, NUV_lambda, NUV_b, NUV_n, SIL1_amp, SIL1_lambda, SIL1_b, + SIL1_n, SIL2_amp, SIL2_lambda, SIL2_b, SIL2_n, FIR_amp, FIR_lambda, + FIR_b, FIR_n)...) + +bounds(::Type{<:P92}) = (10, 10000000) + +function (law::P92)(wave::T) where T + checkbounds(law, wave) || return zero(float(T)) + + x = aa_to_invum(wave) + lam = 1.0 / x # wavelength is in microns + + axav = _p92_single_term(lam, law.BKG_amp, law.BKG_lambda, law.BKG_b, law.BKG_n) + axav += _p92_single_term(lam, law.FUV_amp, law.FUV_lambda, law.FUV_b, law.FUV_n) + axav += _p92_single_term(lam, law.NUV_amp, law.NUV_lambda, law.NUV_b, law.NUV_n) + axav += _p92_single_term(lam, law.SIL1_amp, law.SIL1_lambda, law.SIL1_b, law.SIL1_n) + axav += _p92_single_term(lam, law.SIL2_amp, law.SIL2_lambda, law.SIL2_b, law.SIL2_n) + axav += _p92_single_term(lam, law.FIR_amp, law.FIR_lambda, law.FIR_b, law.FIR_n) + + return axav +end + +# function for calculating a single P92 term +function _p92_single_term(x::Real, amplitude::Real, cen_wave::Real, b::Real, n::Real) + l_norm = x / cen_wave + return amplitude / (l_norm^n + inv(l_norm^n) + b) +end diff --git a/test/fittable_laws.jl b/test/fittable_laws.jl index f3e6117..c1d1389 100644 --- a/test/fittable_laws.jl +++ b/test/fittable_laws.jl @@ -31,3 +31,46 @@ @test ustrip.(reddening) ≈ ref_values rtol = 1e-4 @test ustrip.(reddening1) ≈ ref_values rtol = 1e-4 end + +@testset "P92" begin + + x_inv_microns = [0.21, 0.29, 0.45, 0.61, 0.80, 1.11, 1.43, 1.82, 2.27, 2.50, 2.91, 3.65, 4.00, 4.17, 4.35, + 4.57, 4.76, 5.00, 5.26, 5.56, 5.88, 6.25, 6.71, 7.18, 7.60, 8.00, 8.50, 9.00, 9.50, 10.00] + + wave = 1e4 ./ x_inv_microns + + MW_exvebv = [-3.02, -2.91, -2.76, -2.58, -2.23, -1.60, -0.78, 0.00, 1.00, 1.30, 1.80, 3.10, 4.19, 4.90, 5.77, + 6.57, 6.23, 5.52, 4.90, 4.65, 4.60, 4.73, 4.99, 5.36, 5.91, 6.55, 7.45, 8.45, 9.80, 11.30] + + Rv = 3.08 + ref_values = MW_exvebv ./ Rv .+ 1 + + model = P92() + model1 = P92(FUV_b = 4) + + # Test out of bounds + bad_waves = [9, 1e8] + @test model.(bad_waves) == zeros(length(bad_waves)) + @test model1.(bad_waves) == zeros(length(bad_waves)) + + # testing main part + @test model.(wave) ≈ ref_values rtol = 0.25 atol = 0.01 + @test model1.(wave) ≈ ref_values rtol = 0.25 atol = 0.01 + + # uncertainties + noise = randn(length(wave)) .* 0.01 + wave_unc = wave .± noise + reddening = @inferred broadcast(model, wave_unc) + reddening1 = @inferred broadcast(model1, wave_unc) + @test Measurements.value.(reddening) ≈ ref_values rtol = 0.25 atol = 0.01 + @test Measurements.value.(reddening1) ≈ ref_values rtol = 0.25 atol = 0.01 + + # Unitful + wave_u = wave * u"angstrom" + reddening = @inferred broadcast(model, wave_u) + reddening1 = @inferred broadcast(model1, wave_u) + @test eltype(reddening) <: Gain + @test eltype(reddening1) <: Gain + @test ustrip.(reddening) ≈ ref_values rtol = 0.25 atol = 0.01 + @test ustrip.(reddening1) ≈ ref_values rtol = 0.25 atol = 0.01 +end diff --git a/test/runtests.jl b/test/runtests.jl index d1cfd05..3491478 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,7 +11,7 @@ include("fittable_laws.jl") include("mixture_laws.jl") @testset "interfaces" begin - for LAW in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, G16, F99, F04, F19, M14] + for LAW in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, G16, F99, F04, F19, M14, P92] @test bounds(LAW) == bounds(LAW()) @test checkbounds(LAW, 1000) == checkbounds(LAW(), 1000) low, high = bounds(LAW)