...
 
Commits (6)
......@@ -25,8 +25,8 @@ grain sizes between 1 nm and 1 micron:
>>> sd_car, sd_sil = cb.sdist.WD01(3.1, 6.0, 'A')
>>> sizes = np.logspace(-9, -6, num=10)*u.m
>>> sd_car(sizes)
<Quantity [7.39939133e+00, 9.02469859e-02, 1.94996300e-02, 1.07158763e-03,
1.18483502e-04, 1.30850238e-05, 1.24648409e-06, 9.43835731e-08,
<Quantity [3.94727090e+01, 2.60581728e-01, 7.53742633e-02, 1.57870904e-03,
1.18599378e-04, 1.30850244e-05, 1.24648409e-06, 9.43835731e-08,
2.41782877e-09, 2.43258206e-15] 1 / m>
>>> sd_sil(sizes)
<Quantity [7.86318728e+000, 6.70676535e-001, 5.73425242e-002,
......
import os as _os
data_dir = _os.path.join(
_os.path.dirname(__file__), 'data'
)
from scipy import interpolate as _interp
import astropy.units as _u
import numpy as _np
import os as _os
class Crefin(object):
......@@ -104,9 +103,9 @@ class Crefin(object):
return _np.broadcast_to(interpolation(lam), (len(a), len(lam)))
else:
rinterpolation = _interp.interp2d(lam, a, n.real,
bounds_error=bounds_error)
bounds_error=bounds_error)
iinterpolation = _interp.interp2d(lam, a, n.imag,
bounds_error=bounds_error)
bounds_error=bounds_error)
# scipy.interpolate.interp2d implicitly sorts its `x` and `y`
# arguments and returns an array that's sorted accordingly. Here we
......@@ -127,9 +126,9 @@ class Crefin(object):
else:
def f(a, lam):
n = rinterpolation(lam, a) + 1j*iinterpolation(lam, a)
if _np.prod(r.shape)*_np.prod(lam.shape) > 1:
if _np.prod(a.shape)*_np.prod(lam.shape) > 1:
rowind_unsorted = _np.argsort(lam)
colind_unsorted = _np.argsort(r)
colind_unsorted = _np.argsort(a)
return n[rowind_unsorted][:, colind_unsorted]
else:
return _np.atleast_2d(n)
......@@ -175,21 +174,3 @@ class SGPAHCrefin(Crefin):
a = None
n = _np.array(n)
return cls.fromData(a, lam, n, bounds_error=True, a_bounds_error=False)
cpaD03 = SGPAHCrefin.fromFiles([
_os.path.join(_os.path.dirname(__file__),
'data/crefin/callindex.out_CpaD03_0.01'),
_os.path.join(_os.path.dirname(__file__),
'data/crefin/callindex.out_CpaD03_0.10')
])
cpeD03 = SGPAHCrefin.fromFiles([
_os.path.join(_os.path.dirname(__file__),
'data/crefin/callindex.out_CpeD03_0.01'),
_os.path.join(_os.path.dirname(__file__),
'data/crefin/callindex.out_CpeD03_0.10')
])
silD03 = SGPAHCrefin.fromFiles([
_os.path.join(_os.path.dirname(__file__),
'data/crefin/callindex.out_silD03')
])
from cosmic_dustbox import sdist
from cosmic_dustbox import sdist as _sdist
import astropy.units as _u
sdists = {
'gra': sdist.PowerLaw(
'gra': _sdist.PowerLaw(
50*_u.angstrom, 0.25*_u.micron, 10**-25.13*_u.cm**2.5, -3.5),
'sil': sdist.PowerLaw(
'sil': _sdist.PowerLaw(
50*_u.angstrom, 0.25*_u.micron, 10**-25.11*_u.cm**2.5, -3.5)
}
###############################################################################
from cosmic_dustbox.crefin import SGPAHCrefin as _SGPAHCrefin
from cosmic_dustbox.config import data_dir as _data_dir
import os as _os
###############################################################################
crefins = {
'cpa': _SGPAHCrefin.fromFiles([
_os.path.join(_data_dir,
'crefin/callindex.out_CpaD03_0.01'),
_os.path.join(_data_dir,
'crefin/callindex.out_CpaD03_0.10')
]),
'cpe': _SGPAHCrefin.fromFiles([
_os.path.join(_data_dir,
'crefin/callindex.out_CpeD03_0.01'),
_os.path.join(_data_dir,
'crefin/callindex.out_CpeD03_0.10')
]),
'sil': _SGPAHCrefin.fromFiles([
_os.path.join(_data_dir,
'crefin/callindex.out_silD03')
])
}
import astropy.units as _u
import numpy as _np
from cosmic_dustbox.sdist import LogNormal as _LogNormal
from cosmic_dustbox.sdist import WD01ExpCutoff as _WD01ExpCutoff
from cosmic_dustbox.config import data_dir as _data_dir
import os as _os
def sdists(RV, bc, case):
"""
Get size distributions as defined in [1]_.
Specifically check Table 1 of [1]_.
Parameters
----------
RV : float
R_V value for which
bc : float
Fraction of carbon atoms locked in very small grains in units of 1e5.
case : string
Case 'A' or case 'B' parameters for size distributions. See [1]_ for
definitions.
Returns
-------
car : SizeDist
Size distribution for carbonaceous grains.
sil : SizeDist
Size distribution for silicate grains.
References
----------
.. [1] Weingartner & Draine 2001ApJ...548..296W
"""
# 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(
_data_dir,
'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_car = \
_LogNormal(
sizeMin,
sizeMax,
rho,
0.4,
0.75*params[2],
3.5*_u.angstrom
) + \
_LogNormal(
sizeMin,
sizeMax,
rho,
0.4,
0.25*params[2],
30*_u.angstrom
)
l_car = _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]
)
car = s_car + l_car
return car, sil
from unittest import TestCase
from .. import WD01
import numpy as np
import astropy.units as u
class Testsdists(TestCase):
def test_smoke_test(self):
WD01.sdists(3.1, 0.0, 'A')
return
def test_wrong_params(self):
with self.assertRaises(ValueError):
WD01.sdists(3.1, 8.0, 'A')
with self.assertRaises(ValueError):
WD01.sdists(3.1, 8.0, 'C')
return
def test_out_of_bounds(self):
for sd in WD01.sdists(3.1, 0.0, 'A'):
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
)
for sd in WD01.sdists(3.1, 6.0, 'A'):
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
)
return
......@@ -238,73 +238,6 @@ class WD01ExpCutoff(SizeDist):
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_car = LogNormal(
sizeMin,
sizeMax,
rho,
0.4,
0.75*params[2],
3.5*_u.angstrom
) + \
LogNormal(
sizeMin,
sizeMax,
rho,
0.4,
0.25*params[2],
30*_u.angstrom
)
l_car = 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_car + l_car, sil
###############################################################################
if __name__ == "__main__":
import doctest
......
from unittest import TestCase
from .. import crefin
from .. import config
import numpy as np
import astropy.units as u
import os
......@@ -16,8 +17,9 @@ class TestFromData(TestCase):
cls.a2 = np.array([1e-9, 1e-8])
cls.lam2 = np.linspace(1, 10, num=20)
cls.n2 = (
np.linspace(0, 8, num=40)
+ 1j*np.linspace(0, 9, num=40)).reshape(2, 20)
np.linspace(0, 8, num=40)
+ 1j*np.linspace(0, 9, num=40)
).reshape(2, 20)
return
@classmethod
......@@ -25,6 +27,9 @@ class TestFromData(TestCase):
return
def test_no_a(self):
"""
Crefin that has no dependence on grain size.
"""
c = self.__class__
a = crefin.Crefin.fromData(c.a1, c.lam1, c.n1)
r = a(
......@@ -39,38 +44,37 @@ class TestFromData(TestCase):
r[0, 0], complex(0, 0))
self.assertAlmostEqual(
r[1, 0], complex(0, 0))
# check bounds error
with self.assertRaises(ValueError) as context:
# check bounds error in wavelength
with self.assertRaises(ValueError):
a(
np.array([1])*u.micron,
np.array([11])*u.micron
)
return
def test_normal(self):
def test_with_a(self):
c = self.__class__
a = crefin.Crefin.fromData(c.a2, c.lam2, c.n2)
crefin.Crefin.fromData(c.a2, c.lam2, c.n2)
return
def test_single_element(self):
c = self.__class__
a = crefin.Crefin.fromData(c.a2, c.lam2, c.n2)
n = a(
np.array([1e-9])*u.micron,
np.array([8])*u.micron
)
self.assertTrue
a(np.array([1e-9])*u.micron, np.array([8])*u.micron)
return
def test_single_element_no_bounds(self):
c = self.__class__
a = crefin.Crefin.fromData(c.a2, c.lam2, c.n2)
a = crefin.Crefin.fromData(c.a2, c.lam2, c.n2,
bounds_error=False)
a(np.array([1e-9])*u.micron, np.array([8])*u.micron)
return
def test_normal_no_bounds(self):
c = self.__class__
a = crefin.Crefin.fromData(c.a2, c.lam2, c.n2,
bounds_error=False)
a(np.array([1e-9])*u.micron, np.array([8])*u.micron)
return
......@@ -79,7 +83,7 @@ class TestSGPAHCrefin(TestCase):
@classmethod
def setUpClass(cls):
cls.basepath = os.path.join(
os.path.dirname(os.path.abspath(__file__)), '../data/crefin')
config.data_dir, 'crefin')
return
@classmethod
......@@ -207,21 +211,13 @@ class TestSGPAHCrefin(TestCase):
self.assertAlmostEqual(r.imag, nim_expect[i][j])
return
def test_fromFiles_cpar(self):
paths = [
os.path.join(self.__class__.basepath, 'callindex.out_CpaD03_0.01'),
os.path.join(self.__class__.basepath, 'callindex.out_CpaD03_0.10')
]
a = crefin.SGPAHCrefin.fromFiles(paths)
return
def test_fromFiles_diff_wavelengths(self):
paths = [
os.path.join(self.__class__.basepath, 'callindex.out_CpeD03_0.01'),
os.path.join(self.__class__.basepath, 'callindex.out_CpaD03_0.10')
]
with self.assertRaises(ValueError) as context:
a = crefin.SGPAHCrefin.fromFiles(paths)
with self.assertRaises(ValueError):
crefin.SGPAHCrefin.fromFiles(paths)
return
def test_fromFiles_same_size(self):
......@@ -229,6 +225,6 @@ class TestSGPAHCrefin(TestCase):
os.path.join(self.__class__.basepath, 'callindex.out_CpeD03_0.01'),
os.path.join(self.__class__.basepath, 'callindex.out_CpeD03_0.01')
]
with self.assertRaises(ValueError) as context:
a = crefin.SGPAHCrefin.fromFiles(paths)
with self.assertRaises(ValueError):
crefin.SGPAHCrefin.fromFiles(paths)
return
......@@ -126,37 +126,3 @@ class TestWD01ExpCutoff(TestCase):
sd(np.array([s.value])*s.unit)[0].value, r.value
)
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):
for sd in sdist.WD01(3.1, 0.0, 'A'):
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
)
for sd in sdist.WD01(3.1, 6.0, 'A'):
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
)
from setuptools import setup
import metadata as md
setup(
name=md.name,
version=md.version,
......