Merge pull request #32 from icweaver/F19

Author: mileslucas

Project.toml

@@ -1,10 +1,11 @@
 name = "DustExtinction"
 uuid = "fb44c06c-c62f-5397-83f5-69249e0a3c8e"
 license = "MIT"
-version = "0.9.0"
+version = "0.10.0"
 
 [deps]
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"

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, f99_invum, f04_invum, FM90
+import DustExtinction: ccm89_ca, ccm89_cb, od94_ca, od94_cb, cal00_invum, ccm89_invum, vcg04_invum, gcc09_invum, f99_invum, f04_invum, f19_invum, FM90
 
 dir = joinpath(@__DIR__, "src", "assets")
 
@@ -80,6 +80,19 @@ xlabel!(L"\mu m ^{-1}")
 ylabel!("E(B-V)")
 savefig(joinpath(dir, "vcg04_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"))
+
 #--------------------------------------------------------------------------------
 # F04
 
@@ -93,6 +106,19 @@ xlabel!(L"x\ \left[\mu m ^{-1}\right]")
 ylabel!(L"A(x)/A(V)")
 savefig(joinpath(dir, "F04_plot.svg"))
 
+#--------------------------------------------------------------------------------
+# F19
+
+w = range(0.3, 8.7, length=1000)
+plot()
+for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
+  m = f19_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, "F19_plot.svg"))
+
 #--------------------------------------------------------------------------------
 # FM90
 
@@ -115,16 +141,3 @@ 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/F19_plot.svg

@@ -0,0 +1,19 @@
+@ARTICLE{2019ApJ...886..108F,
+       author = {{Fitzpatrick}, E.~L. and {Massa}, Derck and {Gordon}, Karl D. and
+         {Bohlin}, Ralph and {Clayton}, Geoffrey C.},
+        title = "{An Analysis of the Shapes of Interstellar Extinction Curves. VII. Milky Way Spectrophotometric Optical-through-ultraviolet Extinction and Its R-dependence}",
+      journal = {\apj},
+     keywords = {Interstellar dust, Interstellar dust extinction, Interstellar dust processes, Interstellar absorption, Interstellar reddening, Interstellar scattering, Astrophysics - Astrophysics of Galaxies, Astrophysics - Solar and Stellar Astrophysics},
+         year = 2019,
+        month = dec,
+       volume = {886},
+       number = {2},
+          eid = {108},
+        pages = {108},
+          doi = {10.3847/1538-4357/ab4c3a},
+archivePrefix = {arXiv},
+       eprint = {1910.08852},
+ primaryClass = {astro-ph.GA},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2019ApJ...886..108F},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}

docs/src/assets/f19.bib

@@ -135,6 +135,7 @@ and is loosely associated with the size of the dust grains in the interstellar m
 - [`GCC09`](@ref)
 - [`F99`](@ref)
 - [`F04`](@ref)
+- [`F19`](@ref)
 
 ### Clayton, Cardelli and Mathis (1989)
 
@@ -192,6 +193,14 @@ F99
 F04
 ```
 
+### Fitzpatrick (2019)
+
+![](assets/F19_plot.svg)
+
+```@docs
+F19
+```
+
 ## API/Reference
 
 ```@docs

docs/src/color_laws.md

@@ -41,6 +41,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) |
+| [`F19`](@ref) | [Fitzpatrick (2019)](https://ui.adsabs.harvard.edu/abs/2019ApJ...886..108F) | [download](assets/f19.bib) |
 
 ## Index
 

docs/src/index.md

@@ -15,6 +15,7 @@ export redden,
        VCG04,
        F99,
        F04,
+       F19,
        # Fittable laws
        FM90,
        # Dust maps
@@ -146,7 +147,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, F99, F04]
+for law in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, F99, F04, F19]
     (l::law)(wavelength::Quantity) = l(ustrip(u"Å", wavelength)) * u"mag"
 end
 
@@ -163,6 +164,13 @@ function __init__()
     "https://sncosmo.github.io/data/dust/SFD_dust_4096_sgp.fits"],
     ["50b6aaad0b880762d0fd081177802dcc17c39d7044a410dd5649e2dfd0503e97",
     "84891a59054adab44a7be54051e4dcf0e66e3f13eee0d845ce3739242f553b83"]))
+    register(DataDep("F19",
+    """
+    Tabulated data from Fitzpatrick et al (2019) extinction model calculation
+    Paper: https://ui.adsabs.harvard.edu/abs/2019ApJ...886..108F
+    """,
+    ["https://raw.githubusercontent.com/karllark/dust_extinction/master/dust_extinction/data/F19_tabulated.dat"],
+    ["34011693065bba25396085823ab1df43c058336a95da7d4fdb567f03103b35bc"]))
 end
 
 end # module

src/DustExtinction.jl

@@ -1,4 +1,4 @@
-using Dierckx
+using Dierckx, DelimitedFiles
 
 # Convenience function for wavelength conversion
 @inline aa_to_invum(wave::Real) = 10000 / wave
@@ -452,3 +452,63 @@ function f04_invum(x::Real, Rv::Real)
             optnir_axebv_y,
         )
 end
+
+"""
+    F19(;Rv=3.1)
+
+Fitzpatrick (2019) 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 [1149, 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.
+
+Fitzpatrick, Massa, Gordon et al. (2019, ApJ, 886, 108) model. Based on a
+sample of stars observed spectroscopically in the optical with HST/STIS.
+
+# References
+[Fitzpatrick (2019)](https://ui.adsabs.harvard.edu/abs/2019ApJ...886..108F/)
+"""
+@with_kw struct F19 <: ExtinctionLaw
+    Rv::Float64 = 3.1
+end
+
+function (law::F19)(wave::T) where T
+    checkbounds(law, wave) || return zero(float(T))
+    x = aa_to_invum(wave)
+    return f19_invum(x, law.Rv)
+end
+
+bounds(::Type{F19}) = (1149.4, 33333.3)
+
+"""
+    DustExtinction.f19_invum(x, Rv)
+
+The algorithm used for the [`F19`](@ref) extinction law, given inverse microns
+and Rv. For more information, seek the original paper.
+"""
+function f19_invum(x::Real, Rv::Real)
+    # read and unpack tabulated data
+    data_x, data_k, data_delta_k, data_sigma_k = let data = readdlm(joinpath(datadep"F19", "F19_tabulated.dat"), skipstart=1)
+        (data[:, i] for i in 1:4)
+    end
+
+    if !(0.3 <= x <= 8.7)
+        error("out of bounds of F19, support is over $(bounds(F19)) angstrom")
+    end
+
+    # compute E(lambda-55)/E(B-55) on the tabulated x points
+    k_rV_tab_x = @. data_k + data_delta_k * (Rv - 3.10) * 0.990
+
+    # setup spline interpolation
+    spl = Spline1D(collect(data_x), collect(k_rV_tab_x), k=3)
+
+    # use spline interpolation to evaluate the curve for the input x values
+    k_rV = spl(x)
+
+    # convert to A(x)/A(55) from E(x-55)/E(44-55)
+    a_rV = k_rV / Rv + 1.0
+
+    return a_rV
+end

src/color_laws.jl

@@ -7,7 +7,8 @@ using DustExtinction: ccm89_invum,
                       ccm89_ca,
                       ccm89_cb,
                       f99_invum,
-                      f04_invum
+                      f04_invum,
+                      f19_invum
 
 @testset "helper" begin
     @test aa_to_invum(10000) ≈ 1
@@ -324,3 +325,41 @@ end
         @test ustrip.(reddening) ≈ ref_values[rv] rtol = 0.016
     end
 end
+
+@testset "F19" begin
+    # x values from Fitzpatrick et al. (2000) Table 3
+
+    x_inv_microns = [1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
+    wave = 1e4 ./ x_inv_microns
+
+    # convert from E(x-V)/E(B-V) to A(x)/A(V)
+    ref_values = Dict(
+        3.1 => @. ([-1.757, -0.629, 0.438, 2.090, 4.139, 5.704, 4.904, 5.684, 7.150] + 3.1) / 3.1
+    )
+
+    # test defaults
+    @test F19().(wave) ≈ ref_values[3.1] rtol = 0.016
+
+    for rv in keys(ref_values)
+        law = F19(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 f19_invum(aa_to_invum(bad_waves[1]), rv)
+        @test_throws ErrorException f19_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/color_laws.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, F99, F04]
+    for LAW in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, F99, F04, F19]
         @test bounds(LAW) == bounds(LAW())
         @test checkbounds(LAW, 1000) == checkbounds(LAW(), 1000)
         low, high = bounds(LAW)