diff --git a/Project.toml b/Project.toml index 6535f8c..5f13b50 100644 --- a/Project.toml +++ b/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" diff --git a/docs/plots.jl b/docs/plots.jl index 98776fe..7228686 100644 --- a/docs/plots.jl +++ b/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 diff --git a/docs/src/assets/F04_plot.svg b/docs/src/assets/F04_plot.svg new file mode 100644 index 0000000..3ca69a0 --- /dev/null +++ b/docs/src/assets/F04_plot.svg @@ -0,0 +1,798 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +2 + + +4 + + +6 + + +8 + + +10 + + +0 + + +2 + + +4 + + +6 + + +8 + + + + + + + + + + + + + + + + +Rv=2.0 + + + +Rv=3.1 + + + +Rv=4.0 + + + +Rv=5.0 + + + +Rv=6.0 + + diff --git a/docs/src/assets/f04.bib b/docs/src/assets/f04.bib new file mode 100644 index 0000000..616ac22 --- /dev/null +++ b/docs/src/assets/f04.bib @@ -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} +} diff --git a/docs/src/color_laws.md b/docs/src/color_laws.md index bf7c3dd..876dabd 100644 --- a/docs/src/color_laws.md +++ b/docs/src/color_laws.md @@ -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 diff --git a/docs/src/index.md b/docs/src/index.md index b990768..9a4c519 100644 --- a/docs/src/index.md +++ b/docs/src/index.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 diff --git a/src/DustExtinction.jl b/src/DustExtinction.jl index 66a6d94..eb3761c 100644 --- a/src/DustExtinction.jl +++ b/src/DustExtinction.jl @@ -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 diff --git a/src/color_laws.jl b/src/color_laws.jl index 5ca9774..071def1 100644 --- a/src/color_laws.jl +++ b/src/color_laws.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 diff --git a/test/color_laws.jl b/test/color_laws.jl index 2aba77f..c9b8800 100644 --- a/test/color_laws.jl +++ b/test/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 diff --git a/test/runtests.jl b/test/runtests.jl index 778ad45..491aae0 100644 --- a/test/runtests.jl +++ b/test/runtests.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