Added F99 model

Project.toml

@@ -5,6 +5,7 @@ version = "0.7.0"
 
 [deps]
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
+Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

docs/plots.jl

@@ -1,5 +1,5 @@
 using Plots, LaTeXStrings
-import DustExtinction: ccm89_ca, ccm89_cb, od94_ca, od94_cb, cal00_invum, ccm89_invum, vcg04_invum, gcc09_invum, FM90
+import DustExtinction: ccm89_ca, ccm89_cb, od94_ca, od94_cb, cal00_invum, ccm89_invum, vcg04_invum, gcc09_invum, FM90, f99_invum, F99
 
 dir = joinpath(@__DIR__, "src", "assets")
 
@@ -102,3 +102,16 @@ plot!(w, m4, label = "FUV rise term")
 xlabel!(L"\mu m ^{-1}")
 ylabel!(L"E(\lambda - V)/E(B - V)")
 savefig(joinpath(dir, "FM90_plot.svg"))
+
+#--------------------------------------------------------------------------------
+# F99
+
+w = range(0.3, 10.0, length=1000)
+plot()
+for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
+  m = f99_invum.(w, rv)
+  plot!(w, m, label="Rv=$rv")
+end
+xlabel!(L"x\ \left[\mu m ^{-1}\right]")
+ylabel!(L"A(x)/A(V)")
+savefig(joinpath(dir, "F99_plot.svg"))

docs/src/assets/f99.bib

@@ -0,0 +1,17 @@
+@ARTICLE{1999PASP..111...63F,
+       author = {{Fitzpatrick}, Edward L.},
+        title = "{Correcting for the Effects of Interstellar Extinction}",
+      journal = {\pasp},
+     keywords = {ISM: DUST, EXTINCTION, Astrophysics},
+         year = 1999,
+        month = jan,
+       volume = {111},
+       number = {755},
+        pages = {63-75},
+          doi = {10.1086/316293},
+archivePrefix = {arXiv},
+       eprint = {astro-ph/9809387},
+ primaryClass = {astro-ph},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/1999PASP..111...63F},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}

docs/src/color_laws.md

@@ -133,6 +133,7 @@ and is loosely associated with the size of the dust grains in the interstellar m
 - [`CAL00`](@ref)
 - [`VCG04`](@ref)
 - [`GCC09`](@ref)
+- [`F99`](@ref)
 
 ### Clayton, Cardelli and Mathis (1989)
 
@@ -174,6 +175,14 @@ VCG04
 GCC09
 ```
 
+### Fitzpatrick (1999)
+
+![](assets/F99_plot.svg)
+
+```@docs
+F99
+```
+
 ## API/Reference
 
 ```@docs

docs/src/index.md

@@ -39,6 +39,7 @@ There are various citations relevant to this work. Please be considerate when us
 |  [`GCC09`](@ref)   | [Gordon, Cartledge, & Clayton (2009)](https://ui.adsabs.harvard.edu/abs/2009ApJ...705.1320G)   | [download](assets/gcc09.bib) |
 |  [`FM90`](@ref)   | [Fitzpatrick & Massa (1990)](https://ui.adsabs.harvard.edu/abs/1990ApJS...72..163F)   | [download](assets/fm90.bib) |
 | [`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) |
 
 ## Index
 

src/DustExtinction.jl

@@ -13,6 +13,7 @@ export redden,
        OD94,
        GCC09,
        VCG04,
+       F99,
        # Fittable laws
        FM90,
        # Dust maps
@@ -144,7 +145,7 @@ include("fittable_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]
+for law in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, F99]
     (l::law)(wavelength::Quantity) = l(ustrip(u"Å", wavelength)) * u"mag"
 end
 

src/color_laws.jl

@@ -1,3 +1,5 @@
+using Dierckx
+
 # Convenience function for wavelength conversion
 @inline aa_to_invum(wave::Real) = 10000 / wave
 
@@ -228,3 +230,141 @@ function gcc09_invum(x::Real, Rv::Real)
 
     return a + b / Rv
 end
+
+# Shape models used by F99 and F04
+function _curve_F99_method(
+    x::Real,
+    Rv::Real,
+    c1::Real,
+    c2::Real,
+    c3::Real,
+    c4::Real,
+    x0::Real,
+    gamma::Real,
+    optnir_axav_x::Vector{<:Real},
+    optnir_axav_y::Vector{<:Real},
+    )
+
+    # x value above which FM90 parametrization used
+    x_cutval_uv = 10000.0 / 2700.0
+
+    # required UV points for spline interpolation
+    x_splineval_uv = 10000.0 ./ [2700.0, 2600.0]
+
+    # add in required spline points, otherwise just spline points
+    if x >= x_cutval_uv
+        xuv = vcat(x_splineval_uv, x)
+    else
+        xuv = x_splineval_uv
+    end
+
+    # FM90 model and values
+    fm90_model = FM90(c1=c1, c2=c2, c3=c3, c4=c4, x0=x0, gamma=gamma)
+    # evaluate model and get results in A(x)/A(V)
+    axav_fm90 = fm90_model.(10000.0 ./ xuv) / Rv .+ 1.0 # Expects xuv in Å
+
+    # ignore the spline points
+    if x >= x_cutval_uv
+        axav = last(axav_fm90)
+    else
+        # **Optical Portion**
+        # using cubic spline anchored in UV, optical, and IR
+
+        # optical/NIR points in input x
+
+        # save spline points
+        y_splineval_uv = axav_fm90
+
+        # spline points
+        x_splineval_optir = optnir_axav_x
+
+        # determine optical/IR values at spline points
+        y_splineval_optir = optnir_axav_y
+        spline_x = vcat(x_splineval_optir, x_splineval_uv)
+        spline_y = vcat(y_splineval_optir, y_splineval_uv)
+        spl = Spline1D(spline_x, spline_y, k=3)
+        axav = spl(x)
+    end
+
+    # return A(x)/A(V)
+    return axav
+end
+
+"""
+    F99(;Rv=3.1)
+
+Fitzpatrick (1999) dust law.
+
+# References
+[Fitzpatrick (1999)](https://ui.adsabs.harvard.edu/abs/1999PASP..111...63F/)
+"""
+@with_kw struct F99 <: ExtinctionLaw
+    Rv::Float64 = 3.1
+end
+
+function (law::F99)(wave::T) where T
+    checkbounds(law, wave) || return zero(float(T))
+    x = aa_to_invum(wave)
+    return f99_invum(x, law.Rv)
+end
+
+bounds(::Type{F99}) = (1000.0, 33333.3)
+
+"""
+    DustExtinction.f99_invum(x, Rv)
+
+The algorithm used for the [`F99`](@ref) extinction law, given inverse microns and Rv. For more information, seek the original paper.
+"""
+function f99_invum(x::Real, Rv::Real)
+    if !(0.3 <= x <= 10.0)
+        error("out of bounds of F99, support is over $(bounds(F99)) angstrom")
+    end
+
+    # constant terms
+    c3 = 3.23
+    c4 = 0.41
+    x0 = 4.596
+    gamma = 0.99
+
+    # terms depending on Rv
+    c2 = -0.824 + 4.717 / Rv
+    # original F99 c1-c2 correlation
+    c1 = 2.030 - 3.007 * c2
+
+    # spline points
+    optnir_axav_x = 10000.0 ./
+        [26500.0, 12200.0, 6000.0, 5470.0, 4670.0, 4110.0]
+
+    # determine optical/IR values at spline points
+    #    Final optical spline point has a leading "-1.208" in Table 4
+    #    of F99, but that does not reproduce Table 3.
+    #    Additional indication that this is not correct is from
+    #    fm_unred.pro
+    #    which is based on FMRCURVE.pro distributed by Fitzpatrick.
+    #    --> confirmation needed?
+    #
+    #    Also, fm_unred.pro has different coeff and # of terms,
+    #    but later work does not include these terms
+    #    --> check with Fitzpatrick?
+    opt_axebv_y = [
+        -0.426 + 1.0044 * Rv,
+        -0.050 + 1.0016 * Rv,
+         0.701 + 1.0016 * Rv,
+         1.208 + 1.0032 * Rv - 0.00033 * (Rv^2),
+    ]
+    nir_axebv_y = [0.265, 0.829] * Rv / 3.1
+    optnir_axebv_y = vcat(nir_axebv_y, opt_axebv_y) / Rv
+
+    return _curve_F99_method(
+            x,
+            Rv,
+            c1,
+            c2,
+            c3,
+            c4,
+            x0,
+            gamma,
+            optnir_axav_x,
+            optnir_axebv_y,
+        )
+end

test/color_laws.jl

@@ -5,7 +5,8 @@ using DustExtinction: ccm89_invum,
                       gcc09_invum,
                       aa_to_invum,
                       ccm89_ca,
-                      ccm89_cb
+                      ccm89_cb,
+                      f99_invum
 
 @testset "helper" begin
     @test aa_to_invum(10000) ≈ 1
@@ -248,3 +249,40 @@ end
         @test ustrip.(reddening) ≈ ref_values[rv] rtol = 0.016
     end
 end
+
+@testset "F99" begin
+    # From Fitzpatrick (1999) Table 3
+
+    x_inv_microns = [0.377, 0.820, 1.667, 1.828, 2.141, 2.433, 3.704, 3.846]
+    wave = 1e4 ./ x_inv_microns
+
+    ref_values = Dict(
+        3.1 => [0.265, 0.829, 2.688, 3.055, 3.806, 4.315, 6.265, 6.591] ./ 3.1,
+    )
+
+    # test defaults
+    @test F99().(wave) ≈ ref_values[3.1] rtol = 0.016
+
+    for rv in keys(ref_values)
+        law = F99(Rv = rv)
+        output = @inferred map(law, wave)
+        @test output ≈ ref_values[rv] rtol = 0.016
+
+        bad_xs = [0.2, 11.0]
+        @test @inferred(map(law, bad_xs)) == zeros(length(bad_xs))
+        @test_throws ErrorException f99_invum(bad_xs[1], rv)
+        @test_throws ErrorException f99_invum(bad_xs[2], rv)
+
+        # uncertainties
+        #noise = randn(length(wave)) .* 0.01
+        #wave_unc = wave .± noise
+        #reddening = @inferred map(law, wave_unc)
+        #@test Measurements.value.(reddening) ≈ ref_values[rv] rtol = 0.016
+
+        # Unitful
+        wave_u = wave * u"angstrom"
+        reddening = @inferred map(law, wave_u)
+        @test eltype(reddening) <: Gain
+        @test ustrip.(reddening) ≈ ref_values[rv] rtol = 0.016
+    end
+end

test/runtests.jl

@@ -10,7 +10,7 @@ include("dust_maps.jl")
 include("fittable_laws.jl")
 
 @testset "interfaces" begin
-    for LAW in [CCM89, OD94, CAL00, GCC09, VCG04, FM90]
+    for LAW in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, F99]
         @test bounds(LAW) == bounds(LAW())
         @test checkbounds(LAW, 1000) == checkbounds(LAW(), 1000)
 	low, high = bounds(LAW)