Skip to content

Commit 74f8efe

Browse files
authored
Merge pull request #34 from icweaver/G16
Added G16 model
2 parents 01be6da + 2072118 commit 74f8efe

File tree

10 files changed

+2006
-8
lines changed

10 files changed

+2006
-8
lines changed

docs/plots.jl

Lines changed: 31 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
using Plots, LaTeXStrings
2-
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
2+
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
33

44
dir = joinpath(@__DIR__, "src", "assets")
55

@@ -119,6 +119,19 @@ xlabel!(L"x\ \left[\mu m ^{-1}\right]")
119119
ylabel!(L"A(x)/A(V)")
120120
savefig(joinpath(dir, "F19_plot.svg"))
121121

122+
#--------------------------------------------------------------------------------
123+
# M14
124+
125+
w = range(0.3, 3.3, length=1000)
126+
plot()
127+
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
128+
m = m14_invum.(w, rv)
129+
plot!(w, m, label="Rv=$rv")
130+
end
131+
xlabel!(L"x\ \left[\mu m ^{-1}\right]")
132+
ylabel!(L"A(x)/A(V)")
133+
savefig(joinpath(dir, "M14_plot.svg"))
134+
122135
#--------------------------------------------------------------------------------
123136
# FM90
124137

@@ -143,14 +156,26 @@ ylabel!(L"E(\lambda - V)/E(B - V)")
143156
savefig(joinpath(dir, "FM90_plot.svg"))
144157

145158
#--------------------------------------------------------------------------------
146-
# M14
159+
# G16
147160

148-
w = range(0.3, 3.3, length=1000)
161+
# Fixed f_A = 1.0, variable Rv
162+
w = range(0.3, 10.0, length=1000)
163+
x = 1e4 ./ w
149164
plot()
150165
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
151-
m = m14_invum.(w, rv)
152-
plot!(w, m, label="Rv=$rv")
166+
m = G16(Rv=rv, f_A=1.0).(x)
167+
plot!(w, m, label="Rv=$rv", legendtitle="f_A = 1.0", legend=:topleft)
153168
end
154169
xlabel!(L"x\ \left[\mu m ^{-1}\right]")
155170
ylabel!(L"A(x)/A(V)")
156-
savefig(joinpath(dir, "M14_plot.svg"))
171+
savefig(joinpath(dir, "G16_fixed_f_A_plot.svg"))
172+
173+
# Fixed Rv = 3.1, variable f_A
174+
plot()
175+
for f_A in [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
176+
m = G16(Rv=3.1, f_A=f_A).(x)
177+
plot!(w, m, label="f_A=$f_A", legendtitle="Rv = 3.1", legend=:topleft)
178+
end
179+
xlabel!(L"x\ \left[\mu m ^{-1}\right]")
180+
ylabel!(L"A(x)/A(V)")
181+
savefig(joinpath(dir, "G16_fixed_Rv_plot.svg"))

docs/src/assets/G16_fixed_Rv_plot.svg

Lines changed: 910 additions & 0 deletions
Loading

docs/src/assets/G16_fixed_f_A_plot.svg

Lines changed: 801 additions & 0 deletions
Loading

docs/src/assets/g16.bib

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
@ARTICLE{2016ApJ...826..104G,
2+
author = {{Gordon}, Karl D. and {Fouesneau}, Morgan and {Arab}, Heddy and
3+
{Tchernyshyov}, Kirill and {Weisz}, Daniel R. and
4+
{Dalcanton}, Julianne J. and {Williams}, Benjamin F. and
5+
{Bell}, Eric F. and {Bianchi}, Luciana and {Boyer}, Martha and
6+
{Choi}, Yumi and {Dolphin}, Andrew and {Girardi}, L{\'e}o and
7+
{Hogg}, David W. and {Kalirai}, Jason S. and {Kapala}, Maria and
8+
{Lewis}, Alexia R. and {Rix}, Hans-Walter and {Sandstrom}, Karin and
9+
{Skillman}, Evan D.},
10+
title = "{The Panchromatic Hubble Andromeda Treasury. XV. The BEAST: Bayesian Extinction and Stellar Tool}",
11+
journal = {\apj},
12+
keywords = {dust, extinction, galaxies: individual: M31, methods: data analysis, methods: statistical, stars: fundamental parameters, Astrophysics - Astrophysics of Galaxies},
13+
year = 2016,
14+
month = aug,
15+
volume = {826},
16+
number = {2},
17+
eid = {104},
18+
pages = {104},
19+
doi = {10.3847/0004-637X/826/2/104},
20+
archivePrefix = {arXiv},
21+
eprint = {1606.06182},
22+
primaryClass = {astro-ph.GA},
23+
adsurl = {https://ui.adsabs.harvard.edu/abs/2016ApJ...826..104G},
24+
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
25+
}

docs/src/color_laws.md

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -227,3 +227,14 @@ DustExtinction.bounds
227227
```@docs
228228
FM90
229229
```
230+
231+
### Mixture Extinction Laws
232+
233+
#### Gordon et al. (2016)
234+
235+
![](assets/G16_fixed_f_A_plot.svg)
236+
![](assets/G16_fixed_Rv_plot.svg)
237+
238+
```@docs
239+
G16
240+
```

docs/src/index.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@ There are various citations relevant to this work. Please be considerate when us
3838
| [`VCG04`](@ref) | [Valencic, Clayton, & Gordon (2004)](https://ui.adsabs.harvard.edu/abs/2004ApJ...616..912V) | [download](assets/vcg04.bib) |
3939
| [`GCC09`](@ref) | [Gordon, Cartledge, & Clayton (2009)](https://ui.adsabs.harvard.edu/abs/2009ApJ...705.1320G) | [download](assets/gcc09.bib) |
4040
| [`FM90`](@ref) | [Fitzpatrick & Massa (1990)](https://ui.adsabs.harvard.edu/abs/1990ApJS...72..163F) | [download](assets/fm90.bib) |
41+
| [`G16`](@ref) | [Gordon et al (2016) ](https://ui.adsabs.harvard.edu/abs/2016ApJ...826..104G) | [download](assets/g16.bib) |
4142
| [`SFD98Map`](@ref) | [Schlegel, Finkbeiner and Davis (1998)](https://ui.adsabs.harvard.edu/abs/1998ApJ...500..525S) | [download](assets/sfd98.bib) |
4243
| [`F99`](@ref) | [Fitzpatrick (1999)](https://ui.adsabs.harvard.edu/abs/1999PASP..111...63F) | [download](assets/f99.bib) |
4344
| [`F04`](@ref) | [Fitzpatrick (2004)](https://ui.adsabs.harvard.edu/abs/2004ASPC..309...33F) | [download](assets/f04.bib) |

src/DustExtinction.jl

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,9 @@ export redden,
1919
M14,
2020
# Fittable laws
2121
FM90,
22+
# Mixture laws
23+
G16,
24+
G03_SMCBar,
2225
# Dust maps
2326
SFD98Map,
2427
ebv_galactic
@@ -139,6 +142,7 @@ include("deprecate.jl")
139142
include("color_laws.jl")
140143
include("dust_maps.jl")
141144
include("fittable_laws.jl")
145+
include("mixture_laws.jl")
142146

143147
# --------------------------------------------------------------------------------
144148
# Here be codegen!
@@ -148,7 +152,7 @@ include("fittable_laws.jl")
148152
# at which point adding `(l::ExtinctionLaw)(wave::Quantity)` is possible, until then
149153
# using this code-gen does the trick but requires manually editing
150154
# instead of providing support for all <: ExtinctionLaw
151-
for law in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, F99, F04, F19, M14]
155+
for law in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, G16, G03_SMCBar, F99, F04, F19, M14]
152156
(l::law)(wavelength::Quantity) = l(ustrip(u"Å", wavelength)) * u"mag"
153157
end
154158

src/mixture_laws.jl

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
using DustExtinction: _curve_F99_method
2+
3+
const g03_obsdata_x = (
4+
0.455, 0.606, 0.800, 1.235, 1.538, 1.818, 2.273, 2.703,
5+
3.375, 3.625, 3.875, 4.125, 4.375, 4.625, 4.875, 5.125,
6+
5.375, 5.625, 5.875, 6.125, 6.375, 6.625, 6.875, 7.125,
7+
7.375, 7.625, 7.875, 8.125, 8.375, 8.625
8+
)
9+
const g03_obsdata_axav = (
10+
0.110, 0.169, 0.250, 0.567, 0.801, 1.000, 1.374, 1.672,
11+
2.000, 2.220, 2.428, 2.661, 2.947, 3.161, 3.293, 3.489,
12+
3.637, 3.866, 4.013, 4.243, 4.472, 4.776, 5.000, 5.272,
13+
5.575, 5.795, 6.074, 6.297, 6.436, 6.992
14+
)
15+
16+
const g03_obsdata_tolerance = 6e-2
17+
18+
const g03_C1, g03_C2, g03_C3, g03_C4 = -4.959, 2.264, 0.389, 0.461
19+
const g03_x0 = 4.6
20+
const g03_gamma = 1.0
21+
22+
const g03_optnir_axav_x = 1.0 ./ (
23+
2.198, 1.65, 1.25, 0.81, 0.65, 0.55, 0.44, 0.37
24+
)
25+
# values at 2.198 and 1.25 changed to provide smooth interpolation
26+
# as noted in Gordon et al. (2016, ApJ, 826, 104)
27+
const g03_optnir_axav_y = (0.11, 0.169, 0.25, 0.567, 0.801, 1.00, 1.374, 1.672)
28+
29+
"""
30+
G03_SMCBar(;Rv=2.74) <Internal function>
31+
32+
Gordon et al. (2003) SMCBar Average Extinction Curve.
33+
34+
The observed A(lambda)/A(V) values at 2.198 and 1.25 microns were changed to
35+
provide smooth interpolation as noted in Gordon et al. (2016, ApJ, 826, 104)
36+
37+
# Reference
38+
[Gordon et al. (2003)](https://ui.adsabs.harvard.edu/abs/2003ApJ...594..279G/)
39+
"""
40+
@with_kw struct G03_SMCBar <: DustExtinction.ExtinctionLaw
41+
Rv::Float64 = 2.74
42+
obsdata_x = g03_obsdata_x
43+
obsdata_axav = g03_obsdata_axav
44+
obsdata_tolerance::Float64 = g03_obsdata_tolerance
45+
end
46+
47+
#G03_SMCBar(Rv) = G03_SMCBar(promote(Rv)...)
48+
#G03_SMCBar(Rv=2.74) = G03_SMCBar(Rv)
49+
50+
bounds(::Type{<:G03_SMCBar}) = (1000.0, 33333.3)
51+
52+
function (law::G03_SMCBar)(wave::T) where T
53+
checkbounds(law, wave) || return zero(float(T))
54+
x = aa_to_invum(wave)
55+
return g03_invum(x, law.Rv)
56+
end
57+
58+
59+
function g03_invum(x::Real, Rv::Real)
60+
if !(0.3 <= x <= 10.0)
61+
error("out of bounds of G03_SMCBar, support is over $(bounds(G03_SMCBar)) angstrom")
62+
end
63+
64+
# return A(x)/A(V)
65+
return _curve_F99_method(
66+
x,
67+
Rv,
68+
g03_C1,
69+
g03_C2,
70+
g03_C3,
71+
g03_C4,
72+
g03_x0,
73+
g03_gamma,
74+
g03_optnir_axav_x,
75+
g03_optnir_axav_y,
76+
)
77+
end
78+
79+
"""
80+
G16(;Rv=3.1, f_A=1.0)
81+
82+
Gordon et al. (2016) Milky Way, LMC, & SMC R(V) and f_A dependent model
83+
84+
Returns E(B-V) in magnitudes at the given wavelength relative to the
85+
extinction. This is mixture model between the F99 R(V) dependent model
86+
(component A) and the [`G03_SMCBar`](@ref) model (component B) The default support is
87+
[1000, 33333] Å. Outside of that range this will return 0. Rv is the selective
88+
extinction and is valid over [2, 6]. A typical value for the Milky Way is 3.1.
89+
90+
# References
91+
[Gordon et al. (2016)](https://ui.adsabs.harvard.edu/abs/2016ApJ...826..104G/)
92+
"""
93+
@with_kw struct G16{T<:Number} <: ExtinctionLaw @deftype T
94+
Rv::Float64 = 3.1
95+
f_A = 1.0
96+
end
97+
98+
#G16(Rv, f_A) = G16(promote(Rv, f_A)...)
99+
#G16(Rv=3.1, f_A=3.1) = G16(Rv, f_A)
100+
101+
bounds(::Type{<:G16}) = (1000.0, 33333.3)
102+
103+
function (law::G16)(wave::T) where T
104+
checkbounds(law, wave) || return zero(float(T))
105+
x = aa_to_invum(wave)
106+
return g16_invum(x, law.Rv, law.f_A)
107+
end
108+
109+
"""
110+
DustExtinction.g16_invum(x, Rv)
111+
The algorithm used for the [`G16`](@ref) extinction law, given inverse microns
112+
and Rv. For more information, seek the original paper.
113+
"""
114+
function g16_invum(x::Real, Rv::Real, f_A::Number)
115+
if !(0.3 <= x <= 10.0)
116+
error("out of bounds of G16, support is over $(bounds(G16)) angstrom")
117+
end
118+
119+
# get the A component extinction model
120+
alav_A = f99_invum(x, Rv)
121+
122+
# get the B component extinction model
123+
alav_B = g03_invum(x, Rv)
124+
125+
# create the mixture model
126+
alav = @. f_A * alav_A + (1.0 - f_A) * alav_B
127+
128+
return alav
129+
end

test/mixture_laws.jl

Lines changed: 91 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,91 @@
1+
using Measurements
2+
using DustExtinction: aa_to_invum,
3+
g03_invum,
4+
g16_invum,
5+
G03_SMCBar,
6+
G16
7+
8+
@testset "G16" begin
9+
10+
function get_fA_1_ref()
11+
# From Fitzpatrick (1999) Table 3
12+
x_invum = [0.377, 0.820, 1.667, 1.828, 2.141, 2.433, 3.704, 3.846]
13+
ref_values = [0.265, 0.829, 2.688, 3.055, 3.806, 4.315, 6.265, 6.591] ./ 3.1
14+
return (x_invum, ref_values)
15+
end
16+
17+
# Test Rv=3.1, fA=1.0 mixture
18+
x_invum, ref_values = get_fA_1_ref()
19+
model = G16(Rv=3.1, f_A=1.0)
20+
wave = 1e4 ./ x_invum
21+
@test model.(wave) ref_values rtol = 2e-3
22+
23+
# Test RV=2.74, f_A=0 mixture
24+
tmodel = G16(Rv=2.74, f_A=0.0)
25+
gmodel = DustExtinction.G03_SMCBar()
26+
x_invum = gmodel.obsdata_x
27+
wave = 1e4 ./ x_invum
28+
ref_values = gmodel.obsdata_axav
29+
tolerance = gmodel.obsdata_tolerance
30+
@test collect(tmodel.(wave)) collect(ref_values) rtol = tolerance
31+
32+
# Test out of bounds
33+
model = G16(Rv=3.1, f_A=1.0)
34+
bad_waves = [100, 4e4]
35+
@test model.(bad_waves) == zeros(length(bad_waves))
36+
@test_throws ErrorException g16_invum(aa_to_invum(bad_waves[1]), 3.1, 1.0)
37+
@test_throws ErrorException g16_invum(aa_to_invum(bad_waves[2]), 3.1, 1.0)
38+
39+
# uncertainties
40+
x_invum, ref_values = get_fA_1_ref()
41+
wave = 1e4 ./ x_invum
42+
model = G16(Rv=3.1, f_A=1.0)
43+
noise = rand(length(wave)) .* 0.01
44+
wave_unc = wave noise
45+
reddening = map(w -> @uncertain(model(w)), wave_unc)
46+
@test Measurements.value.(reddening) ref_values rtol = 1e-3
47+
48+
# Unitful
49+
model = G16(Rv=3.1, f_A=1.0)
50+
x_invum, ref_values = get_fA_1_ref()
51+
wave = 1e4 ./ x_invum
52+
wave_u = wave * u"angstrom"
53+
reddening = @inferred map(model, wave_u)
54+
@test eltype(reddening) <: Gain
55+
@test ustrip.(reddening) ref_values rtol = 0.016
56+
57+
end
58+
59+
@testset "G03_SMCBar" begin
60+
61+
# Test output
62+
model = G03_SMCBar()
63+
x, y, tolerance = model.obsdata_x, model.obsdata_axav, model.obsdata_tolerance
64+
w = 1e4 ./ x
65+
@test model.(collect(w)) collect(y) rtol = tolerance
66+
67+
# Test out of bounds
68+
model = G03_SMCBar(Rv=2.74)
69+
bad_waves = [100, 4e4]
70+
@test model.(bad_waves) == zeros(length(bad_waves))
71+
@test_throws ErrorException g03_invum(aa_to_invum(bad_waves[1]), 3.1)
72+
@test_throws ErrorException g03_invum(aa_to_invum(bad_waves[2]), 3.1)
73+
74+
# uncertainties
75+
model = G03_SMCBar(Rv=2.74)
76+
x, y, tolerance = model.obsdata_x, model.obsdata_axav, model.obsdata_tolerance
77+
wave = 1e4 ./ x
78+
noise = rand(length(wave)) .* 0.01
79+
wave_unc = wave noise
80+
reddening = map(w -> @uncertain(model(w)), wave_unc)
81+
@test Measurements.value.(reddening) collect(y) rtol = tolerance
82+
83+
# Unitful
84+
model = G03_SMCBar(Rv=2.74)
85+
x, y, tolerance = model.obsdata_x, model.obsdata_axav, model.obsdata_tolerance
86+
wave = 1e4 ./ x
87+
wave_u = collect(wave) * u"angstrom"
88+
reddening = @inferred map(model, wave_u)
89+
@test eltype(reddening) <: Gain
90+
@test ustrip.(reddening) collect(y) rtol = tolerance
91+
end

test/runtests.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,9 +8,10 @@ include("deprecate.jl")
88
include("color_laws.jl")
99
include("dust_maps.jl")
1010
include("fittable_laws.jl")
11+
include("mixture_laws.jl")
1112

1213
@testset "interfaces" begin
13-
for LAW in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, F99, F04, F19, M14]
14+
for LAW in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, G16, F99, F04, F19, M14]
1415
@test bounds(LAW) == bounds(LAW())
1516
@test checkbounds(LAW, 1000) == checkbounds(LAW(), 1000)
1617
low, high = bounds(LAW)

0 commit comments

Comments
 (0)