Skip to content

GitLab

C
cosmic_dustbox
Commit 438d4f66 authored by Martin Glatzle's avatar Martin Glatzle
Browse files
Options

Update size dists.

parent e567d46f
Hide whitespace changes
Inline Side-by-side
Showing with 152 additions and 24 deletions
cosmic_dustbox/data/sdist/WD01-2001ApJ...548..296W-Table1.txt 0 → 100644
View file @ 438d4f66
# Table 1 of 2001ApJ...548..296W which gives size distribution parameters for different lines of sight and assumptions.
# R_V 10^5 b_C b_C Case \alpha_g \beta_g a_{t,g} [m] a_{c,g} [m] C_g \alpha_s \beta_s a_{t,s} [m] a_{c,s} [m] C_s
3.1 0.0 0.0 A -2.25 -0.0648 0.00000000745 0.000000606 9.94e-11 -1.48 -9.34 0.000000172 0.0000001 1.02e-12
3.1 1.0 0.00001 A -2.17 -0.0382 0.00000000373 0.000000586 3.79e-10 -1.46 -10.3 0.000000174 0.0000001 1.09e-12
3.1 2.0 0.00002 A -2.04 -0.111 0.00000000828 0.000000543 5.57e-11 -1.43 -11.7 0.000000173 0.0000001 1.27e-12
3.1 3.0 0.00003 A -1.91 -0.125 0.00000000837 0.000000499 4.15e-11 -1.41 -11.5 0.000000171 0.0000001 1.33e-12
3.1 4.0 0.00004 A -1.84 -0.132 0.00000000898 0.000000489 2.90e-11 -2.1 -0.114 0.000000169 0.0000001 1.26e-13
3.1 5.0 0.00005 A -1.72 -0.322 0.0000000254 0.000000438 3.20e-12 -2.1 -0.0407 0.000000166 0.0000001 1.27e-13
3.1 6.0 0.00006 A -1.54 -0.165 0.0000000107 0.000000428 9.99e-12 -2.21 0.3 0.000000164 0.0000001 1.0e-13
4.0 0.0 0.0 A -2.26 -0.199 0.0000000241 0.000000861 5.47e-12 -2.03 0.668 0.000000189 0.0000001 5.2e-14
4.0 1.0 0.00001 A -2.16 -0.0862 0.00000000867 0.000000803 4.58e-11 -2.05 0.832 0.000000188 0.0000001 4.81e-14
4.0 2.0 0.00002 A -2.01 -0.0973 0.00000000811 0.000000696 3.96e-11 -2.06 0.995 0.000000185 0.0000001 4.7e-14
4.0 3.0 0.00003 A -1.83 -0.175 0.0000000117 0.000000604 1.42e-11 -2.08 1.29 0.000000184 0.0000001 4.26e-14
4.0 4.0 0.00004 A -1.64 -0.247 0.0000000152 0.000000536 5.83e-12 -2.09 1.58 0.000000183 0.0000001 3.94e-14
5.5 0.0 0.0 A -2.35 -0.668 0.000000148 0.00000196 4.82e-14 -1.57 1.1 0.000000198 0.0000001 4.24e-14
5.5 1.0 0.00001 A -2.12 -0.67 0.0000000686 0.00000135 3.65e-13 -1.57 1.25 0.000000197 0.0000001 4.00e-14
5.5 2.0 0.00002 A -1.94 -0.853 0.0000000786 0.000000921 2.57e-13 -1.55 1.33 0.000000195 0.0000001 4.05e-14
5.5 3.0 0.00003 A -1.61 -0.722 0.0000000418 0.00000072 7.58e-13 -1.59 2.12 0.000000193 0.0000001 3.2e-14
cosmic_dustbox/sdist.py
View file @ 438d4f66
###############################################################################
import astropy.units as _u
import numpy as _np
import os as _os
from scipy.special import erf as _erf
###############################################################################
......@@ -169,9 +171,15 @@ class PowerLaw(SizeDist):
return
class Log_Normal(SizeDist):
class DoubleLogNormal(SizeDist):
"""
Sum of two log normal distributions as in Eq. (2) of [1]_.
def __init__(self, sizeMin, sizeMax, rho, sgma, bc, b, a0):
References
----------
.. [1] Weingartner & Draine 2001ApJ...548..296W
"""
def __init__(self, sizeMin, sizeMax, rho, sgma, bc, a0):
# mass of carbon atom
m_C = 12.0107*_u.u
......@@ -181,9 +189,9 @@ class Log_Normal(SizeDist):
denominator = (
(2*_np.pi**2)**1.5 * rho * a0[i]**3 * sgma *
(1 +
_np.erf((3 * sgma/_np.sqrt(2)) +
(_np.log((a0[i]/3.5/_u.angstrom).decompose().value) /
(sgma * _np.sqrt(2)))
_erf((3 * sgma/_np.sqrt(2)) +
(_np.log((a0[i]/3.5/_u.angstrom).decompose().value) /
(sgma * _np.sqrt(2)))
)
)
)
......@@ -195,7 +203,7 @@ class Log_Normal(SizeDist):
_B = [B_val(i) for i in range(2)]
def f(a):
return sum([_B(i)/a * _np.exp(-0.5*(
return sum([_B[i]/a * _np.exp(-0.5*(
_np.log((a/a0[i]).decompose().value)/sgma)**2)
for i in range(2)])
......@@ -203,30 +211,95 @@ class Log_Normal(SizeDist):
return
class WD01_dst(SizeDist):
class WD01ExpCutoff(SizeDist):
"""
Power law distribution with exponential cutoff as in Eqs. (4) and (5) of
[1]_.
def __init__(self, sizeMin, sizeMax, a_t, beta, a_c, alpha, C):
References
----------
.. [1] Weingartner & Draine 2001ApJ...548..296W
"""
def __init__(self, sizeMin, sizeMax, alpha, beta, a_t, a_c, C):
if beta >= 0:
def F(a):
return 1 + (beta * a) / a_t
else:
def F(a):
return 1 / (1 - (beta * a) / a_t)
if beta >= 0:
def F(a):
return 1 + (beta * a) / a_t
else:
def F(a):
return 1 / (1 - (beta * a) / a_t)
def exp_func(a):
r = _np.ones_like(a.value)
ind = _np.where(a > a_t)
r[ind] = _np.exp(-(((a[ind] - a_t)/a_c).decompose().value)**3)
return r
def exp_func(a):
r = _np.ones_like(a.value)
ind = _np.where(a > a_t)
r[ind] = _np.exp(-(((a[ind] - a_t)/a_c).decompose().value)**3)
return r
def f(a):
return C*(a / a_t)**alpha/a \
* F(a) * exp_func(a)
def f(a):
return C*(a / a_t)**alpha/a \
* F(a) * exp_func(a)
super().__init__(sizeMin, sizeMax, f)
return
super().__init__(sizeMin, sizeMax, f)
return
def WD01(RV, bc, case):
# pandas dataframe would be far better suited for this but don't know if
# pulling in another dependency makes sense just for this
data = _np.genfromtxt(
_os.path.join(
_os.path.dirname(__file__),
'data/sdist/WD01-2001ApJ...548..296W-Table1.txt'
),
dtype=None,
encoding=None
)
# doing == float comparison here, might come back to bite us at some point
params = data[
[
j for j, x in enumerate(data)
if (x[0] == RV and x[1] == bc and x[3] == case)
]
]
if len(params) > 1:
# this should never happen
raise ValueError("Could not uniquely identify parameter set.")
elif len(params) == 0:
raise ValueError("Could not identify parameter set.")
params = params[0]
sizeMin = 3.5*_u.angstrom
sizeMax = 10*_u.micron
rho = 2.24*_u.g/_u.cm**3
s_gra = DoubleLogNormal(
sizeMin,
sizeMax,
rho,
0.4,
[0.75*params[2], 0.25*params[2]],
[3.5*_u.angstrom, 30*_u.angstrom]
)
l_gra = WD01ExpCutoff(
sizeMin,
sizeMax,
params[4],
params[5],
params[6]*_u.m,
params[7]*_u.m,
params[8]
)
sil = WD01ExpCutoff(
sizeMin,
sizeMax,
params[9],
params[10],
params[11]*_u.m,
params[12]*_u.m,
params[13]
)
return s_gra, l_gra, sil
###############################################################################
if __name__ == "__main__":
......
cosmic_dustbox/tests/test_sdist.py
View file @ 438d4f66
......@@ -79,3 +79,39 @@ class TestSdist(TestCase):
np.isclose(b(self.__class__._sizes), r, rtol=1e-10, atol=0)
))
return
class TestWD01(TestCase):
def test_smoke_test(self):
sdist.WD01(3.1, 0.0, 'A')
return
def test_wrong_params(self):
with self.assertRaises(ValueError) as context:
sdist.WD01(3.1, 8.0, 'A')
with self.assertRaises(ValueError) as context:
sdist.WD01(3.1, 8.0, 'C')
return
def test_out_of_bounds(self):
sg, lg, sil = sdist.WD01(3.1, 0.0, 'A')
for sd in [sg, lg, sil]:
self.assertEqual(
sd(np.array([1])*u.angstrom)[0],
0.0/u.angstrom
)
self.assertEqual(
sd(np.array([11])*u.micron)[0],
0.0/u.angstrom
)
sg, lg, sil = sdist.WD01(3.1, 6.0, 'A')
for sd in [sg, lg, sil]:
self.assertEqual(
sd(np.array([1])*u.angstrom)[0],
0.0/u.angstrom
)
self.assertEqual(
sd(np.array([11])*u.micron)[0],
0.0/u.angstrom
)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment