Merge pull request #29 from siddharthlal25/p92

Added P92 model

Author: icweaver

docs/src/assets/bibtex/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

@@ -243,6 +243,16 @@ lplot(FM90) # hide
 FM90
 ```
 
+### Pei (1992)
+
+```@example plotting
+lplot(P92) # hide
+```
+
+```@docs
+P92
+```
+
 ## Mixture Extinction Laws
 
 ### Gordon et al. (2003)

docs/src/index.md

@@ -48,6 +48,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/bibtex/sfd98.bib) |
 | [`F99`](@ref)      | [Fitzpatrick (1999)](https://ui.adsabs.harvard.edu/abs/1999PASP..111...63F)                    | [download](assets/bibtex/f99.bib)   |
 | [`F04`](@ref)      | [Fitzpatrick (2004)](https://ui.adsabs.harvard.edu/abs/2004ASPC..309...33F)                    | [download](assets/bibtex/f04.bib)   |
+| [`P92`](@ref)      | [Pei (1992)](https://ui.adsabs.harvard.edu/abs/1992ApJ...395..130P)                            | [download](assets/bibtex/P92.bib)   |
 | [`F19`](@ref)      | [Fitzpatrick (2019)](https://ui.adsabs.harvard.edu/abs/2019ApJ...886..108F)                    | [download](assets/bibtex/f19.bib)   |
 | [`M14`](@ref)      | [Maiz Apellaniz et al. (2014)](https://ui.adsabs.harvard.edu/abs/2014A%26A...564A..63M)        | [download](assets/bibtex/m14.bib)   |
 

docs/src/plotting.jl

@@ -70,6 +70,36 @@ function lplot(law::Type{<:FM90}; args...)
     fig
 end
 
+"""
+    lplot(law::P92; args...)
+
+Fittable law series plot with automatic axis labels.
+"""
+function lplot(law::Type{P92}; args...)
+    # m1
+    fig, ax, p = lplot(law(); label="total")
+
+    m2 = P92(FUV_amp=0, NUV_amp=0, SIL1_amp=0, SIL2_amp=0, FIR_amp=0)
+    Makie.lines!(ax, m2; label="BKG only")
+
+    m3 = P92(NUV_amp=0, SIL1_amp=0, SIL2_amp=0, FIR_amp=0)
+    Makie.lines!(ax, m3; label="BKG+FUV only")
+
+    m4 = P92(FUV_amp=0, SIL1_amp=0, SIL2_amp=0, FIR_amp=0)
+    Makie.lines!(ax, m4; label="BKG+NUV only")
+
+    m5 = P92(FUV_amp=0, NUV_amp=0, SIL2_amp=0)
+    Makie.lines!(ax, m5; label="BKG+FIR+SIL1 only")
+
+    m6 = P92(FUV_amp=0, NUV_amp=0, SIL1_amp=0)
+    Makie.lines!(ax, m6; label="BKG+FIR+SIL2 only")
+
+    m7 = P92(FUV_amp=0, NUV_amp=0, SIL1_amp=0, SIL2_amp=0)
+    Makie.lines!(ax, m7; label="BKG+FIR only")
+
+    fig
+end
+
 """
     mplot(law::Type{G16}, Rvs, f_A::Real; args...)
 

src/DustExtinction.jl

@@ -20,6 +20,7 @@ export redden,
        M14,
        # Fittable laws
        FM90,
+       P92,
        # Mixture laws
        G16,
        G03_SMCBar,

src/fittable_laws.jl

@@ -80,3 +80,162 @@ function (law::FM90)(wave::T) where T <: Real
 
     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.08,  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.561019978746464
+
+julia> P92(FUV_b = 2.0).([1000, 2000, 3000])
+3-element Vector{Float64}:
+ 5.17952309549434
+ 2.7581232249728607
+ 1.8081781540687367
+```
+
+## Default Parameter Values
+
+| Term  | lambda (μm) | A                     | b     | n   |
+| :---: | :---------: | :-------------------: | :---: | :-: |
+| BKG   | 0.047       | 218.57142857142858    | 90.0  | 2.0 |
+| FUV   | 0.08        | 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)
+"""
+Base.@kwdef struct P92{T<:Number} <: ExtinctionLaw
+    BKG_amp::T = 165.0 * (1 / 3.08 + 1)
+    BKG_lambda::T = 0.047
+    BKG_b::T = 90.0
+    BKG_n::T = 2.0
+    FUV_amp::T = 14.0 * (1 / 3.08 + 1)
+    FUV_lambda::T = 0.08
+    FUV_b::T = 4.0
+    FUV_n::T = 6.5
+    NUV_amp::T = 0.045 * (1 / 3.08 + 1)
+    NUV_lambda::T = 0.22
+    NUV_b::T = -1.95
+    NUV_n::T = 2.0
+    SIL1_amp::T = 0.002 * (1 / 3.08 + 1)
+    SIL1_lambda::T = 9.7
+    SIL1_b::T = -1.95
+    SIL1_n::T = 2.0
+    SIL2_amp::T = 0.002 * (1 / 3.08 + 1)
+    SIL2_lambda::T = 18.0
+    SIL2_b::T = -1.80
+    SIL2_n::T = 2.0
+    FIR_amp::T = 0.012 * (1 / 3.08 + 1)
+    FIR_lambda::T = 25.0
+    FIR_b::T = 0.00
+    FIR_n::T = 2.0
+    function 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)
+
+        BKG_amp ≥ 0 || error("`BKG_amp` must be ≥ 0, got ", BKG_amp)
+        FUV_amp ≥ 0 || error("`FUV_amp` must be ≥ 0, got ", FUV_amp)
+        0.06 ≤ FUV_lambda ≤ 0.08 || error("`FUV_lambda` must be in between [0.06, 0.08], got ", FUV_lambda)
+        NUV_amp ≥ 0 || error("`NUV_amp` must be ≥ 0, got", NUV_amp)
+        0.20 ≤ NUV_lambda ≤ 0.24 || error("`NUV_lambda` must be in between [0.20, 0.24], got ", NUV_lambda)
+        SIL1_amp ≥ 0 || error("`SIL1_amp` must be ≥ 0, got ", SIL1_amp)
+        7 ≤ SIL1_lambda ≤ 13 || error("`SIL1_lambda` must be in between [7, 13], got ", SIL1_lambda)
+        SIL2_amp ≥ 0 || error("`SIL2_amp` must be ≥ 0, got ", SIL2_amp)
+        15 ≤ SIL2_lambda ≤ 21 || error("`SIL2_lambda` must be in between [15, 21], got ", SIL2_lambda)
+        FIR_amp ≥ 0 || error("`FIR_amp` must be ≥ 0, got ", FIR_amp)
+        20 ≤ FIR_lambda ≤ 30 || error("`FIR_lambda` must be in between [20, 30], got ", FIR_lambda)
+
+        params = 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)
+
+        return new{eltype(params)}(params...)
+    end
+end
+
+bounds(::Type{<:P92}) = (10, 10_000_000)
+
+function (law::P92)(wave::T) where T <: Real
+    checkbounds(law, wave) || return zero(float(T))
+
+    x = aa_to_invum(wave)
+    lam = 1 / 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
+    l_norm_n = l_norm^n
+    return amplitude / (l_norm_n + inv(l_norm_n) + b)
+end

test/fittable_laws.jl

@@ -34,3 +34,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

@@ -13,7 +13,7 @@ using Test, Measurements, Unitful, UnitfulAstro, Random
     VERSION ≥ v"1.9" && include("makie_recipes.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)