Merge pull request #30 from icweaver/F99

Project.toml

@@ -1,17 +1,19 @@
 name = "DustExtinction"
 uuid = "fb44c06c-c62f-5397-83f5-69249e0a3c8e"
 license = "MIT"
-version = "0.7.0"
+version = "0.8.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"
 UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f"
 
 [compat]
 DataDeps = "0.7"
+Dierckx = "0.4"
 FITSIO = "0.13.0, 0.14, 0.15"
 Parameters = "0.12"
 Unitful = "0.17.0, 1"

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, f99_invum, FM90
 
 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,143 @@ function gcc09_invum(x::Real, Rv::Real)
 
     return a + b / Rv
 end
+
+# x value above which FM90 parametrization used
+const f99_x_cutval_uv = aa_to_invum(2700)
+# required UV points for spline interpolation
+const f99_x_splineval_uv = aa_to_invum.((2700, 2600))
+
+# Shape models used by F99 and F04
+function _curve_F99_method(
+    x,
+    Rv,
+    c1,
+    c2,
+    c3,
+    c4,
+    x0,
+    gamma,
+    optnir_axav_x,
+    optnir_axav_y,
+    )
+
+    # add in required spline points, otherwise just spline points
+    if x >= f99_x_cutval_uv
+        xuv = (f99_x_splineval_uv..., x)
+    else
+        xuv = f99_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(aa_to_invum(xuv)) / Rv + 1
+
+    # ignore the spline points
+    if x >= f99_x_cutval_uv
+        axav = last(axav_fm90)
+    else
+        # **Optical Portion**
+
+        # save spline points
+        y_splineval_uv = axav_fm90
+
+        # spline points
+        x_splineval_optir = (0.0, optnir_axav_x...)
+
+        # determine optical/IR values at spline points
+        y_splineval_optir = (0.0, optnir_axav_y...)
+        spline_x = (x_splineval_optir..., f99_x_splineval_uv...)
+        spline_y = (y_splineval_optir..., y_splineval_uv...)
+        spl = Spline1D(collect(spline_x), collect(spline_y), k=3)
+        axav = spl(x)
+    end
+
+    # return A(x)/A(V)
+    return axav
+end
+
+# spline points
+const f99_optnir_axav_x = aa_to_invum.((26500, 12200, 6000, 5470, 4670, 4110))
+const f99_nir_axebv_y_params = @. (0.265, 0.829) / 3.1
+
+# c1-c2 correlation terms
+const f99_c3 = 3.23
+const f99_c4 = 0.41
+const f99_x0 = 4.596
+const f99_gamma = 0.99
+
+"""
+    F99(;Rv=3.1)
+
+Fitzpatrick (1999) dust law.
+
+Returns E(B-V) in magnitudes at the given wavelength relative to the
+extinction. This model applies to the UV and optical to NIR spectral range.
+The default support is [1000, 33333] Å. Outside of that range this will return
+0. Rv is the selective extinction and is valid over [2, 6]. A typical value for
+the Milky Way is 3.1.
+
+# 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
+
+    # terms depending on Rv
+    c2 = @evalpoly (1. / Rv) -0.824 4.717
+    # original F99 c1-c2 correlation
+    c1 = @evalpoly c2 2.030 -3.007
+
+    # 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 =
+        (@evalpoly Rv -0.426 1.0044),
+        (@evalpoly Rv -0.050 1.0016),
+        (@evalpoly Rv  0.701 1.0016),
+        (@evalpoly Rv  1.208 1.0032 -0.00033)
+
+    nir_axebv_y = @. f99_nir_axebv_y_params * Rv
+    optnir_axebv_y = @. (nir_axebv_y..., opt_axebv_y...) / Rv
+
+    return _curve_F99_method(
+            x,
+            Rv,
+            c1,
+            c2,
+            f99_c3,
+            f99_c4,
+            f99_x0,
+            f99_gamma,
+            f99_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_waves = [100, 4e4]
+        @test @inferred(map(law, bad_waves)) == zeros(length(bad_waves))
+        @test_throws ErrorException f99_invum(aa_to_invum(bad_waves[1]), rv)
+        @test_throws ErrorException f99_invum(aa_to_invum(bad_waves[2]), rv)
+
+        # uncertainties
+        noise = rand(length(wave)) .* 0.01
+        wave_unc = wave .± noise
+        reddening = map(w -> @uncertain(law(w)), wave_unc)
+        @test Measurements.value.(reddening) ≈ ref_values[rv] rtol = 1e-3
+
+        # 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)