Merge pull request #31 from icweaver/F04

Author: mileslucas

Project.toml

@@ -1,7 +1,7 @@
 name = "DustExtinction"
 uuid = "fb44c06c-c62f-5397-83f5-69249e0a3c8e"
 license = "MIT"
-version = "0.8.0"
+version = "0.9.0"
 
 [deps]
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"

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, FM90
+import DustExtinction: ccm89_ca, ccm89_cb, od94_ca, od94_cb, cal00_invum, ccm89_invum, vcg04_invum, gcc09_invum, f99_invum, f04_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"))
 
+#--------------------------------------------------------------------------------
+# F04
+
+w = range(0.3, 10.0, length=1000)
+plot()
+for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
+  m = f04_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, "F04_plot.svg"))
+
 #--------------------------------------------------------------------------------
 # FM90
 

docs/src/assets/F04_plot.svg

@@ -0,0 +1,17 @@
+@INPROCEEDINGS{2004ASPC..309...33F,
+       author = {{Fitzpatrick}, E.~L.},
+        title = "{Interstellar Extinction in the Milky Way Galaxy}",
+     keywords = {Astrophysics},
+    booktitle = {Astrophysics of Dust},
+         year = 2004,
+       editor = {{Witt}, Adolf N. and {Clayton}, Geoffrey C. and {Draine}, Bruce T.},
+       series = {Astronomical Society of the Pacific Conference Series},
+       volume = {309},
+        month = may,
+        pages = {33},
+archivePrefix = {arXiv},
+       eprint = {astro-ph/0401344},
+ primaryClass = {astro-ph},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2004ASPC..309...33F},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}

docs/src/assets/f04.bib

@@ -134,6 +134,7 @@ and is loosely associated with the size of the dust grains in the interstellar m
 - [`VCG04`](@ref)
 - [`GCC09`](@ref)
 - [`F99`](@ref)
+- [`F04`](@ref)
 
 ### Clayton, Cardelli and Mathis (1989)
 
@@ -183,6 +184,14 @@ GCC09
 F99
 ```
 
+### Fitzpatrick (2004)
+
+![](assets/F04_plot.svg)
+
+```@docs
+F04
+```
+
 ## API/Reference
 
 ```@docs

docs/src/color_laws.md

@@ -40,6 +40,7 @@ There are various citations relevant to this work. Please be considerate when us
 |  [`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) |
+| [`F04`](@ref) | [Fitzpatrick (2004)](https://ui.adsabs.harvard.edu/abs/2004ASPC..309...33F) | [download](assets/f04.bib) |
 
 ## Index
 

docs/src/index.md

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

src/DustExtinction.jl

@@ -298,15 +298,12 @@ 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/)
 """
@@ -324,7 +321,6 @@ 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)
@@ -370,3 +366,89 @@ function f99_invum(x::Real, Rv::Real)
             optnir_axebv_y,
         )
 end
+
+# spline points
+const f04_opt_axav_x = aa_to_invum.((6000, 5470, 4670, 4110))
+const f04_nir_axav_x = (0.50, 0.75, 1.0)
+
+# **Use NIR spline x values in FM07, clipped to K band for now
+const f04_optnir_axav_x = (f04_nir_axav_x..., f04_opt_axav_x...)
+
+# c1-c2 correlation terms
+const f04_c3 = 2.991
+const f04_c4 = 0.319
+const f04_x0 = 4.592
+const f04_gamma = 0.922
+
+"""
+    F04(;Rv=3.1)
+Fitzpatrick (2004) 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.
+Equivalent to the F99 model with an updated NIR Rv dependence
+See also Fitzpatrick & Massa (2007, ApJ, 663, 320)
+# References
+[Fitzpatrick (2004)](https://ui.adsabs.harvard.edu/abs/2004ASPC..309...33F/)
+"""
+@with_kw struct F04 <: ExtinctionLaw
+    Rv::Float64 = 3.1
+end
+
+function (law::F04)(wave::T) where T
+    checkbounds(law, wave) || return zero(float(T))
+    x = aa_to_invum(wave)
+    return f04_invum(x, law.Rv)
+end
+
+bounds(::Type{F04}) = (1000.0, 33333.3)
+
+"""
+    DustExtinction.f04_invum(x, Rv)
+The algorithm used for the [`F04`](@ref) extinction law, given inverse microns and Rv. For more information, seek the original paper.
+"""
+function f04_invum(x::Real, Rv::Real)
+    if !(0.3 <= x <= 10.0)
+        error("out of bounds of F04, support is over $(bounds(F04)) angstrom")
+    end
+
+    # original F99 Rv dependence
+    c2 = @evalpoly (1. / Rv) -0.824 4.717
+    # updated F04 C1-C2 correlation
+    c1 = @evalpoly c2 2.18 -2.91
+
+    # **Keep optical spline points from F99:
+    #    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?
+    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)
+
+    # updated NIR curve from F04, note R dependence
+    # Julia v1.0 workaround: https://github.com/JuliaAstro/DustExtinction.jl/pull/31#commitcomment-40605778
+    nir_axebv_y_coeff = (0.63 * Rv - 0.84)
+    nir_axebv_y = @. nir_axebv_y_coeff * f04_nir_axav_x^1.84
+
+    optnir_axebv_y = @. (nir_axebv_y..., opt_axebv_y...) / Rv
+
+    return _curve_F99_method(
+            x,
+            Rv,
+            c1,
+            c2,
+            f04_c3,
+            f04_c4,
+            f04_x0,
+            f04_gamma,
+            f04_optnir_axav_x,
+            optnir_axebv_y,
+        )
+end

src/color_laws.jl

@@ -6,7 +6,8 @@ using DustExtinction: ccm89_invum,
                       aa_to_invum,
                       ccm89_ca,
                       ccm89_cb,
-                      f99_invum
+                      f99_invum,
+                      f04_invum
 
 @testset "helper" begin
     @test aa_to_invum(10000) ≈ 1
@@ -286,3 +287,40 @@ end
         @test ustrip.(reddening) ≈ ref_values[rv] rtol = 0.016
     end
 end
+
+@testset "F04" 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.185, 0.772, 2.688, 3.055, 3.805, 4.315, 6.456, 6.781] ./ 3.1,
+    )
+
+    # test defaults
+    @test F04().(wave) ≈ ref_values[3.1] rtol = 0.016
+
+    for rv in keys(ref_values)
+        law = F04(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 f04_invum(aa_to_invum(bad_waves[1]), rv)
+        @test_throws ErrorException f04_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,11 +10,11 @@ include("dust_maps.jl")
 include("fittable_laws.jl")
 
 @testset "interfaces" begin
-    for LAW in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, F99]
+    for LAW in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, F99, F04]
         @test bounds(LAW) == bounds(LAW())
         @test checkbounds(LAW, 1000) == checkbounds(LAW(), 1000)
-	low, high = bounds(LAW)
-	@test all(checkbounds.(LAW, rand(low:high, 1000)))
+        low, high = bounds(LAW)
+        @test all(checkbounds.(LAW, rand(low:high, 1000)))
     end
     @test bounds(DustExtinction.ExtinctionLaw) == (0, Inf)
 end