Added P92 model

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
 

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}
+}

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)

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) |
 

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
 

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

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

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)