merged in master/M14 (01be6da)

Author: icweaver

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, f19_invum, FM90, G16
+import DustExtinction: ccm89_ca, ccm89_cb, od94_ca, od94_cb, cal00_invum, ccm89_invum, vcg04_invum, gcc09_invum, f99_invum, f04_invum, f19_invum, m14_invum, FM90, G16
 
 dir = joinpath(@__DIR__, "src", "assets")
 
@@ -119,6 +119,19 @@ xlabel!(L"x\ \left[\mu m ^{-1}\right]")
 ylabel!(L"A(x)/A(V)")
 savefig(joinpath(dir, "F19_plot.svg"))
 
+#--------------------------------------------------------------------------------
+# M14
+
+w = range(0.3, 3.3, length=1000)
+plot()
+for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
+  m = m14_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, "M14_plot.svg"))
+
 #--------------------------------------------------------------------------------
 # FM90
 
@@ -166,5 +179,3 @@ end
 xlabel!(L"x\ \left[\mu m ^{-1}\right]")
 ylabel!(L"A(x)/A(V)")
 savefig(joinpath(dir, "G16_fixed_Rv_plot.svg"))
-
-

docs/src/assets/M14_plot.svg

@@ -0,0 +1,21 @@
+@ARTICLE{2014A&A...564A..63M,
+       author = {{Ma{\'\i}z Apell{\'a}niz}, J. and {Evans}, C.~J. and {Barb{\'a}}, R.~H. and
+         {Gr{\"a}fener}, G. and {Bestenlehner}, J.~M. and {Crowther}, P.~A. and
+         {Garc{\'\i}a}, M. and {Herrero}, A. and {Sana}, H. and
+         {Sim{\'o}n-D{\'\i}az}, S. and {Taylor}, W.~D. and {van Loon}, J. Th. and
+         {Vink}, J.~S. and {Walborn}, N.~R.},
+        title = "{The VLT-FLAMES Tarantula Survey. XVI. The optical and NIR extinction laws in 30 Doradus and the photometric determination of the effective temperatures of OB stars}",
+      journal = {\aap},
+     keywords = {open clusters and associations: individual: 30 Doradus, dust, extinction, Magellanic Clouds, stars: early-type, stars: fundamental parameters, Astrophysics - Astrophysics of Galaxies, Astrophysics - Cosmology and Nongalactic Astrophysics, Astrophysics - Instrumentation and Methods for Astrophysics, Astrophysics - Solar and Stellar Astrophysics},
+         year = 2014,
+        month = apr,
+       volume = {564},
+          eid = {A63},
+        pages = {A63},
+          doi = {10.1051/0004-6361/201423439},
+archivePrefix = {arXiv},
+       eprint = {1402.3062},
+ primaryClass = {astro-ph.GA},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2014A&A...564A..63M},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}

docs/src/assets/m14.bib

@@ -201,6 +201,14 @@ F04
 F19
 ```
 
+### Maiz Apellaniz et al. (2014)
+
+![](assets/m14_plot.svg)
+
+```@docs
+M14
+```
+
 ## API/Reference
 
 ```@docs

docs/src/color_laws.md

@@ -43,6 +43,7 @@ There are various citations relevant to this work. Please be considerate when us
 | [`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) |
+| [`M14`](@ref) | [Maiz Apellaniz et al. (2014)](https://ui.adsabs.harvard.edu/abs/2014A%26A...564A..63M) | [download](assets/m14.bib) |
 
 ## Index
 

docs/src/index.md

@@ -16,6 +16,7 @@ export redden,
        F99,
        F04,
        F19,
+       M14,
        # Fittable laws
        FM90,
        # Mixture laws
@@ -151,7 +152,7 @@ include("mixture_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, G16, G03_SMCBar, F99, F04, F19]
+for law in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, G16, G03_SMCBar, F99, F04, F19, M14]
     (l::law)(wavelength::Quantity) = l(ustrip(u"Å", wavelength)) * u"mag"
 end
 

src/DustExtinction.jl

@@ -512,3 +512,107 @@ function f19_invum(x::Real, Rv::Real)
 
     return a_rV
 end
+
+# Optical
+const m14_x1 = (1.0)
+const m14_xi1 = first(m14_x1)
+const m14_x2 = (1.15, 1.81984, 2.1, 2.27015, 2.7)
+const m14_x3 = (3.5, 3.9, 4.0, 4.1, 4.2)
+const m14_xi3 = last(m14_x3)
+const m14_a1v = 0.574 * m14_x1^1.61
+const m14_b1v = -0.527 * m14_x1^1.61
+const m14_x2_diff = m14_x2 .- 1.82
+const m14_a2v = @. (
+    1
+    + 0.17699 * (m14_x2_diff)
+    - 0.50447 * (m14_x2_diff)^2
+    - 0.02427 * (m14_x2_diff)^3
+    + 0.72085 * (m14_x2_diff)^4
+    + 0.01979 * (m14_x2_diff)^5
+    - 0.77530 * (m14_x2_diff)^6
+    + 0.32999 * (m14_x2_diff)^7
+    + (0.0, 0.0, -0.011, 0.0, 0.0)
+)
+const m14_b2v = @. (
+      1.41338 * (m14_x2 - 1.82)
+    + 2.28305 * (m14_x2 - 1.82)^2
+    + 1.07233 * (m14_x2 - 1.82)^3
+    - 5.38434 * (m14_x2 - 1.82)^4
+    - 0.62251 * (m14_x2 - 1.82)^5
+    + 5.30260 * (m14_x2 - 1.82)^6
+    - 2.09002 * (m14_x2 - 1.82)^7
+    + (0.0, 0.0, 0.091, 0.0, 0.0)
+)
+const m14_a3v = @. (
+      1.752
+    - 0.316 * m14_x3
+    - 0.104 / ((m14_x3 - 4.67)^2 + 0.341)
+    + (0.442, 0.341, 0.130, 0.020, 0.000)
+)
+const m14_b3v = @. (
+     -3.090
+    + 1.825 * m14_x3
+    + 1.206 / ((m14_x3 - 4.62)^2 + 0.263)
+    - (1.256, 1.021, 0.416, 0.064, 0.000)
+)
+const m14_xn = collect((m14_x1..., m14_x2..., m14_x3...))
+const m14_anv = collect((m14_a1v..., m14_a2v..., m14_a3v...))
+const m14_bnv = collect((m14_b1v..., m14_b2v..., m14_b3v...))
+const m14_a_spl = Spline1D(m14_xn, m14_anv, bc="nearest")
+const m14_b_spl = Spline1D(m14_xn, m14_bnv, bc="nearest")
+
+"""
+    M14(;Rv=3.1)
+
+Maiz Apellaniz et al (2014) Milky Way & LMC R(V) dependent model.
+
+Returns E(B-V) in magnitudes at the given wavelength relative to the
+extinction. The published UV extinction curve is identical to Clayton,
+Cardelli, and Mathis (1989, CCM). Forcing the optical section to match smoothly
+with CCM introduces a non-physical feature at high values of R5495 around 3.9
+inverse microns; see section 5 in Maiz Apellaniz et al. (2014) for more
+discussion. For that reason, we provide the M14 model only through 3.3 inverse
+microns, the limit of the optical in CCM. 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. R5495 = A(5485)/E(4405-5495) Spectral
+equivalent to photometric R(V).
+
+# References
+[Maiz Apellaniz et al. (2014)](https://ui.adsabs.harvard.edu/abs/2014A%26A...564A..63M/)
+"""
+@with_kw struct M14 <: ExtinctionLaw
+    Rv::Float64 = 3.1
+end
+
+function (law::M14)(wave::T) where T
+    checkbounds(law, wave) || return zero(float(T))
+    x = aa_to_invum(wave)
+    return m14_invum(x, law.Rv)
+end
+
+bounds(::Type{M14}) = (3030.3, 33333.3)
+
+"""
+    DustExtinction.m14_invum(x, Rv)
+
+The algorithm used for the [`M14`](@ref) extinction law, given inverse microns and Rv. For more information, seek the original paper.
+"""
+function m14_invum(x::Real, Rv::Real)
+    if !(0.3 <= x <= 3.3)
+        error("out of bounds of M14, support is over $(bounds(M14)) angstrom")
+    end
+
+    a = zero(x)
+    b = zero(x)
+    if x < m14_xi1
+        # Infrared
+        a = 0.574 * x^1.61
+        b = -0.527 * x^1.61
+    elseif x < m14_xi3
+        # Optical
+        a = m14_a_spl(x)
+        b = m14_b_spl(x)
+    end
+
+    return a + b / Rv
+end

src/color_laws.jl

@@ -8,7 +8,8 @@ using DustExtinction: ccm89_invum,
                       ccm89_cb,
                       f99_invum,
                       f04_invum,
-                      f19_invum
+                      f19_invum,
+                      m14_invum
 
 @testset "helper" begin
     @test aa_to_invum(10000) ≈ 1
@@ -363,3 +364,40 @@ end
         @test ustrip.(reddening) ≈ ref_values[rv] rtol = 0.016
     end
 end
+
+@testset "M14" begin
+    # using R5495 = 3.1
+
+    x_inv_microns = [0.5, 2.0]
+    wave = 1e4 ./ x_inv_microns
+
+    ref_values = Dict(
+        3.1 => [0.1323, 1.141]
+    )
+
+    # test defaults
+    @test M14().(wave) ≈ ref_values[3.1] rtol = 0.016
+
+    for rv in keys(ref_values)
+        law = M14(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 m14_invum(aa_to_invum(bad_waves[1]), rv)
+        @test_throws ErrorException m14_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

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