...
 
Commits (36)
"""
Toolbox for cosmic dust species and mixtures and their properties.
"""
from . import sdist
from . import crefin
from scipy import interpolate as _interp
import astropy.units as _u
import numpy as _np
import os as _os
class Crefin(object):
"""
Base class used to represent the complex refractive index of a certain
grain material possibly along a certain axis.
This class implements a `__call__` method which computes the crefin on a
grid of grain size and photon energy (or equivalent).
Parameters
----------
f : callable
A callable object that takes as first argument a 1-D array of grain
sizes in micron and as second argument a 1-D array of wavelengths in
micron. It shall return a 2-D array of the complex refractive index
evaluated at the given sizes and wavelengths. The first axis of this
array shall correspond to grain size and the second axis to
wavelength. Note that this axis ordering is reversed wrt to the output
of `np.meshgrid`.
Attributes
----------
_f : callable
Directly taken from parameters.
"""
def __init__(self, f):
"""
See class docstring.
"""
self._f = f
return
def __call__(self, a, wave):
"""
Compute the complex refractive index as a function of grain size and
wavelength or equivalent.
Parameters
----------
a : array-valued Quantity [L]
1-D array of grain sizes as which to evaluate crefin.
wave : array-valued Quantity
1-D array of photon wavelengths or equivalent at which to evaluate
crefin.
Returns
-------
: array
2-D array of crefin values. The first axis corresponds to grain
size and the second to photon wavelength.
"""
return self._f(
a.to(_u.micron, equivalencies=_u.spectral()).value,
wave.to(_u.micron, equivalencies=_u.spectral()).value)
@classmethod
def fromData(cls, a, lam, n, bounds_error=True, a_bounds_error=False):
"""
Linearly interpolate complex refractive index from given data.
Create an object which interpolates the complex refractive index in
grain size and wavelength from the given data in its `__call__` method.
Parameters
----------
a, lam : array
1-D arrays of grain sizes and wavelengths in micron. Must have
length >= 2. `a` can also be None if there is no size dependence or
it is not known.
n : array
2-D array of the complex refractive index. Note: its first index is
grain size and its second index is wavelength, so if len(a) == l
and len(lam) == k, n.shape == (l, k). If a is None, `n` shall be
1-D with the same length as `lam`.
bounds_error : bool, optional
Whether to raise out of bounds errors when interpolating. If False,
extrapolation is used. See also `a_bounds_error`.
a_bounds_error : bool, optional
If `bounds_error` is True, this can be used to specify if bounds
errors on grain size should be ignored by setting it False. This is
useful because often `n` is not known at many grain sizes. Note
that if the bounds in grain size are to be ignored, the values of
`n` for the largest/smallest grain size will be used outside the
bounds instead of extrapolating. If `bounds_error` is False,
`a_bounds_error` will be ignored, i.e. in this case extrapolation
is used as for the wavelengths. If `a` is None, it is assumed that
there is no size dependence of `n`.
Returns
-------
: cls
An instance of the calling class the `__call__` method of which
interpolates the crefin on the given data.
"""
if a is None:
interpolation = _interp.interp1d(lam, n, bounds_error=bounds_error)
def f(a, lam):
return _np.broadcast_to(interpolation(lam), (len(a), len(lam)))
else:
rinterpolation = _interp.interp2d(lam, a, n.real,
bounds_error=bounds_error)
iinterpolation = _interp.interp2d(lam, a, n.imag,
bounds_error=bounds_error)
# scipy.interpolate.interp2d implicitly sorts its `x` and `y`
# arguments and returns an array that's sorted accordingly. Here we
# shuffle the rows and columns of the returned array around so that
# their row-(column-)order corresponds to the order of the input
# arrays. This is done for convenience, however it is inefficient
# and might be improved in the future.
if bounds_error and not a_bounds_error:
def f(aa, lam):
r = _np.clip(aa, _np.amin(a), _np.amax(a))
n = rinterpolation(lam, r) + 1j*iinterpolation(lam, r)
if _np.prod(r.shape)*_np.prod(lam.shape) > 1:
rowind_unsorted = _np.argsort(r)
colind_unsorted = _np.argsort(lam)
return n[rowind_unsorted][:, colind_unsorted]
else:
return _np.atleast_2d(n)
else:
def f(a, lam):
n = rinterpolation(lam, a) + 1j*iinterpolation(lam, a)
if _np.prod(r.shape)*_np.prod(lam.shape) > 1:
rowind_unsorted = _np.argsort(lam)
colind_unsorted = _np.argsort(r)
return n[rowind_unsorted][:, colind_unsorted]
else:
return _np.atleast_2d(n)
return cls(f)
class SGPAHCrefin(Crefin):
@classmethod
def parseCrefinFile(cls, path):
with open(path) as f:
data = _np.loadtxt(f, skiprows=5)
f.seek(0)
next(f)
size = float(f.readline().split('=')[0])
return size, data
@classmethod
def fromFiles(cls, paths):
a = []
n = []
for j, path in enumerate(paths):
size, data = cls.parseCrefinFile(path)
if size in a:
raise ValueError("Grain size "
+ str(size) + "provided twice.")
else:
a.append(size)
if j == 0:
lam = data[:, 0]
else:
if not _np.array_equal(lam, data[:, 0]):
raise ValueError(
"Wavelengths in file " + path
+ "differ from those in the other files visited "
"so far.")
nt = data[:, 3] + 1 + 1j*data[:, 4]
n.append(nt)
if len(a) > 1:
a = _np.array(a)
elif len(a) == 1:
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')
])
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
# 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
###############################################################################
import numpy as _np
import astropy.units as u
###############################################################################
class GrainSpecies(object):
"""
Grain species base class.
"""
def __init__(self):
return
def scatt(self, sizes, wave, **kwargs):
return self._scatt(
sizes.to(_u.micron),
wave.to(_u.micron, equivalencies=_u.spectral()),
**kwargs
)
class SphericalGS(GrainSpecies):
"""
Base class for spherical grain species.
"""
def __init__(self):
return
def volume(self, sizes):
return 4*_np.pi/3 * size**3
def geomCS(self, sizes):
return _np.pi * sizes**2
def surfaceArea(self, sizes):
return 4*_np.pi * sizes**2
class HomSphericalGS(SphericalGS):
"""
Base class for homogeneous spherical grain species.
Represents grain species for which Mie theory and/or anomalous diffraction
theory can be used to compute optical properties.
"""
def __init__(self):
return
def mass(self, sizes):
return self.volume(sizes)*self.massDensity
class IsoHomSphericalGS(HomSphericalGS):
"""
Base class for isotropic and homogeneous spherical grain species.
For these grain species, Mie theory and anomalous diffraction theory can be
directly employed.
"""
def __init__(self, crefin, mieSolver):
self.crefin = crefin
return
def _scatt(self, sizes, wave):
return
class cAxisSphericalGS(HomSphericalGS):
"""
Base class for homogeneous spherical grain species composed of anisotropic
material.
Represents grain species that are spherical but composed of a material that
has crystal axes so that its properties depend on orientation. It is
assumed that they can be treated in an "effective" way using Mie
theory/anomalous diffraction theory by averaging properties along different
axes appropriately. This is then supposed to represent the average
properties of an ensemble of randomly oriented grains.
The classic example of this are graphite grains where the crystal axis is
perpendicular to the so-called basal plane and properties are computed
using the "1/3 - 2/3 approximation", which has been shown to work well
[1]_.
References
----------
.. [1] Draine & Malhotra 1993ApJ...414..632D
"""
def __init__(self, crefins, axesWeights, mieSolver):
self.crefins = crefins
self.axesWeights = axesWeights
self.solver = mieSolver
return
def _scatt(self, sizes, wave):
return
from cosmic_dustbox import sdist
from cosmic_dustbox import Gspcs as gs
import astropy.units as _u
import os as _os
import numpy as _np
import sys
import pandas
Path = _os.path.dirname(_os.path.dirname(_os.path.realpath(__file__))) + \
'\\lib\\paramsWD01a.csv'
sys.path.insert(0, Path)
class WD01(object):
with open(Path)as f:
R_V = [] ; bc=[] ; alpha_g = [] ; beta_g = [] ; a_tg = [] ; a_cg = []
C_g=[] ; alpha_s = []; beta_s = []; a_ts = [] ; a_cs = [] ; C_s =[]
df = pandas.read_csv(f)
R_V.append(df['R_V'])
bc.append(df['b_C'])
alpha_g.append(df['\\alpha_g'])
beta_g.append(df['\\beta_g'])
a_tg.append(df['a_{t,g} [m]'])
a_cg.append(df['a_{c,g} [m]'])
C_g.append(df['C_g'])
alpha_s.append(df['\\alpha_s'])
beta_s.append(df['\\beta_s'])
a_ts.append(df['a_{t,s} [m]'])
a_cs.append(df['a_{c,s} [m]'])
C_s.append(df['C_s'])
rho_g = 2.24 *_u.g/_u.cm**3
b = _np.array([0.75 * bc[1],0.25 * bc[2]])
#
def __init__(self, Indx):
def load_file(Indx, st):
sdists = {
'gra': sdist.Log_Normal(3.5*_u.Angstrom, _np.inf, rho_g,0.4,bc,b,a0) + \
sdist.WD01_dst(3.5*_u.Angstrom, _np.inf , a_tg[Indx],
beta_g[Indx],a_cg[Indx], alpha_g[Indx],C_g[Indx]) , \
'sil': sdist.WD01(
sdist.WD01_dst(3.5*_u.Angstrom, _np.inf,a_ts[Indx], beta_s[Indx], \
a_cs[Indx],alpha_s[Indx],C_s[Indx])
}
return sdist[st]
self.sdist_gr = load_file(Indx,'gra')
self.sdist_sil = load_file(Indx,'sil')
return
#
\ No newline at end of file
###############################################################################
import astropy.units as _u
import numpy as _np
import os as _os
from scipy.special import erf as _erf
###############################################################################
......@@ -83,7 +86,7 @@ class SizeDist(object):
Overloading of ``+`` operator. Can add an instance of ``SizeDist`` (or
subclass thereof), any callable or a scalar to ``self.func``. No check
on units is performed. Returns an instance of ``self.__class__``, the
on units is performed. Returns an instance of ``SizeDist``, the
``func`` attribute of which is defined to return the corresponding sum.
If ``other`` is an instance of ``SizeDist`` (or subclass), take maximum
......@@ -121,15 +124,15 @@ class SizeDist(object):
def func(sizes):
return self.func(sizes) + other.func(sizes)
return self.__class__(sizeMin, sizeMax, func)
return SizeDist(sizeMin, sizeMax, func)
elif callable(other):
def func(sizes):
return self.func(sizes) + other(sizes)
return self.__class__(self.sizeMin, self.sizeMax, func)
return SizeDist(self.sizeMin, self.sizeMax, func)
else:
def func(sizes):
return other + self.func(sizes)
return self.__class__(self.sizeMin, self.sizeMax, func)
return SizeDist(self.sizeMin, self.sizeMax, func)
# make addition commutative
__radd__ = __add__
......@@ -164,10 +167,140 @@ class PowerLaw(SizeDist):
def __init__(self, sizeMin, sizeMax, power, C):
def f(a):
return C*a**power
super(self).__init__(sizeMin, sizeMax, f)
super().__init__(sizeMin, sizeMax, f)
return
class DoubleLogNormal(SizeDist):
"""
Sum of two log normal distributions as in Eq. (2) of [1]_.
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
def B_val(i):
nominator = (3 * _np.exp(-4.5 * sgma**2) * bc[i] * m_C)
denominator = (
(2*_np.pi**2)**1.5 * rho * a0[i]**3 * sgma *
(1 +
_erf((3 * sgma/_np.sqrt(2)) +
(_np.log((a0[i]/3.5/_u.angstrom).decompose().value) /
(sgma * _np.sqrt(2)))
)
)
)
# FIXME: what do we do if the denominator is zero
if denominator != 0:
B = (nominator/denominator).decompose().value
return B
_B = [B_val(i) for i in range(2)]
def f(a):
return sum([_B[i]/a * _np.exp(-0.5*(
_np.log((a/a0[i]).decompose().value)/sgma)**2)
for i in range(2)])
super().__init__(sizeMin, sizeMax, f)
return
class WD01ExpCutoff(SizeDist):
"""
Power law distribution with exponential cutoff as in Eqs. (4) and (5) of
[1]_.
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)
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)
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_car = DoubleLogNormal(
sizeMin,
sizeMax,
rho,
0.4,
[0.75*params[2], 0.25*params[2]],
[3.5*_u.angstrom, 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
import numpy as np
import astropy.units as u
import os
class TestFromData(TestCase):
@classmethod
def setUpClass(cls):
cls.a1 = None
cls.lam1 = np.linspace(1, 10, num=10)
cls.n1 = np.linspace(0, 10, num=10) + 1j*np.linspace(0, 10, num=10)
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)
return
@classmethod
def tearDownClass(cls):
return
def test_no_a(self):
c = self.__class__
a = crefin.Crefin.fromData(c.a1, c.lam1, c.n1)
r = a(
np.array([1, 2])*u.micron,
np.array([1, 5, 9])*u.micron
)
self.assertTrue(
r.shape[0] == 2)
self.assertTrue(
r.shape[1] == 3)
self.assertAlmostEqual(
r[0, 0], complex(0, 0))
self.assertAlmostEqual(
r[1, 0], complex(0, 0))
# check bounds error
with self.assertRaises(ValueError) as context:
a(
np.array([1])*u.micron,
np.array([11])*u.micron
)
return
def test_normal(self):
c = self.__class__
a = 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
return
def test_single_element_no_bounds(self):
c = self.__class__
a = crefin.Crefin.fromData(c.a2, c.lam2, c.n2)
return
def test_normal_no_bounds(self):
c = self.__class__
a = crefin.Crefin.fromData(c.a2, c.lam2, c.n2,
bounds_error=False)
return
class TestSGPAHCrefin(TestCase):
@classmethod
def setUpClass(cls):
cls.basepath = os.path.join(
os.path.dirname(os.path.abspath(__file__)), '../data/crefin')
return
@classmethod
def tearDownClass(cls):
return
def test_parseCrefinFile_asil(self):
path = os.path.join(self.__class__.basepath, 'callindex.out_silD03')
size, data = crefin.SGPAHCrefin.parseCrefinFile(path)
self.assertAlmostEqual(0.1, size)
self.assertEqual(data.shape[0], 837)
self.assertEqual(data.shape[1], 5)
return
def test_parseCrefinFile_cpar(self):
path = os.path.join(self.__class__.basepath,
'callindex.out_CpaD03_0.01')
size, data = crefin.SGPAHCrefin.parseCrefinFile(path)
self.assertAlmostEqual(0.01, size)
self.assertEqual(data.shape[0], 386)
self.assertEqual(data.shape[1], 5)
return
def test_parseCrefinFile_cperp(self):
path = os.path.join(self.__class__.basepath,
'callindex.out_CpeD03_0.01')
size, data = crefin.SGPAHCrefin.parseCrefinFile(path)
self.assertAlmostEqual(0.01, size)
self.assertEqual(data.shape[0], 383)
self.assertEqual(data.shape[1], 5)
return
def test_fromFiles_asil(self):
path = os.path.join(self.__class__.basepath, 'callindex.out_silD03')
c = crefin.SGPAHCrefin.fromFiles([path])
a = np.array([0.1])*u.micron
wave = np.array([1.239840e5, 5.481167e2, 6.199200e-5])*u.micron
nre_expect = [2.4350+1, 2.4040+1, -1.9450e-6+1]
nim_expect = [1.1190e-3, 8.2550e-2, 1.7830e-8]
n = c(a, wave)
self.assertEqual(n.shape[0], len(a))
self.assertEqual(n.shape[1], len(wave))
for j, r in enumerate(n[0, :]):
self.assertAlmostEqual(r.real, nre_expect[j])
self.assertAlmostEqual(r.imag, nim_expect[j])
return
def test_fromFiles_asil_multiple_a(self):
path = os.path.join(self.__class__.basepath, 'callindex.out_silD03')
c = crefin.SGPAHCrefin.fromFiles([path])
a = np.array([1.0, 0.1, 0.01])*u.micron
wave = np.array([1.239840e5, 5.481167e2, 6.199200e-5])*u.micron
nre_expect = [2.4350+1, 2.4040+1, -1.9450e-6+1]
nim_expect = [1.1190e-3, 8.2550e-2, 1.7830e-8]
n = c(a, wave)
self.assertEqual(n.shape[0], len(a))
self.assertEqual(n.shape[1], len(wave))
for i in range(len(a)):
for j, r in enumerate(n[i, :]):
self.assertAlmostEqual(r.real, nre_expect[j])
self.assertAlmostEqual(r.imag, nim_expect[j])
return
def test_fromFiles_cperp(self):
paths = [
os.path.join(self.__class__.basepath, 'callindex.out_CpeD03_0.01'),
os.path.join(self.__class__.basepath, 'callindex.out_CpeD03_0.10')
]
c = crefin.SGPAHCrefin.fromFiles(paths)
a = np.array([1.0, 0.1, 0.01, 0.001])*u.micron
wave = np.array([1.239840E+05, 6.322489E+00, 6.199200E-05])*u.micron
nre_expect = np.array([
[1.0596E+03, 4.5430E+00, 1.5405E-10],
[1.0596E+03, 4.5430E+00, 1.5405E-10],
[3.4137E+02, 4.5568E+00, 1.5405E-10],
[3.4137E+02, 4.5568E+00, 1.5405E-10],
])+1
nim_expect = [
[1.0636E+03, 4.2997E+00, 1.4471E-17],
[1.0636E+03, 4.2997E+00, 1.4471E-17],
[3.4149E+02, 4.3122E+00, 1.3963E-16],
[3.4149E+02, 4.3122E+00, 1.3963E-16],
]
n = c(a, wave)
self.assertEqual(n.shape[0], len(a))
self.assertEqual(n.shape[1], len(wave))
for i in [0, 1]:
for j, r in enumerate(n[i, :]):
self.assertAlmostEqual(r.real, nre_expect[i][j],
msg='Failed at i='+str(i)+' j='+str(j))
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')
]
c = crefin.SGPAHCrefin.fromFiles(paths)
a = np.array([1.0, 0.1, 0.01, 0.001])*u.micron
wave = np.array([1.239840E+05, 6.199200E+00, 6.199200E-05])*u.micron
nre_expect = np.array([
[1.5051E+02, 1.1266E+00, 2.3165E-10],
[1.5051E+02, 1.1266E+00, 2.3165E-10],
[1.4329E+02, 1.0261E+00, 2.3165E-10],
[1.4329E+02, 1.0261E+00, 2.3165E-10],
])+1
nim_expect = [
[1.5156E+02, 2.3038E-02, 1.4161E-17],
[1.5156E+02, 2.3038E-02, 1.4161E-17],
[1.4434E+02, 3.0850E+00, 1.5614E-17],
[1.4434E+02, 3.0850E+00, 1.5614E-17],
]
n = c(a, wave)
self.assertEqual(n.shape[0], len(a))
self.assertEqual(n.shape[1], len(wave))
for i in [0, 1]:
for j, r in enumerate(n[i, :]):
self.assertAlmostEqual(r.real, nre_expect[i][j],
msg='Failed at i='+str(i)+' j='+str(j))
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)
return
def test_fromFiles_same_size(self):
paths = [
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)
return
......@@ -63,7 +63,9 @@ class TestSdist(TestCase):
r = a(self.__class__._sizes)
r[np.where(r != 0)] = r[np.where(r != 0)] + scalar
self.assertTrue(
np.all(b(self.__class__._sizes) == r))
np.all(
np.isclose(b(self.__class__._sizes), r, rtol=1e-10, atol=0)
))
return
def test_add_scalar_quantity_float(self):
......@@ -73,5 +75,41 @@ class TestSdist(TestCase):
r = a(self.__class__._sizes)
r[np.where(r != 0)] = r[np.where(r != 0)] + scalar
self.assertTrue(
np.all(b(self.__class__._sizes) == r))
np.all(
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):
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
)
......@@ -3,7 +3,7 @@
\usetikzlibrary{math}
\begin{document}
\begin{tikzpicture}
\draw[thick] (0, 0) rectangle (10, 10);
\draw (0, 0) rectangle (10, 10);
% red circles
% import numpy as np
......
......@@ -5,10 +5,10 @@
\usetikzlibrary{positioning}
\begin{document}
\begin{tikzpicture}
\draw[thick] (0, 0) rectangle (5, 5);
\draw[thick] (5, 0) rectangle (10, 5);
\draw[thick] (0, 5) rectangle (5, 10);
\draw[thick] (5, 5) rectangle (10, 10);
\draw (0, 0) rectangle (5, 5);
\draw (5, 0) rectangle (10, 5);
\draw (0, 5) rectangle (5, 10);
\draw (5, 5) rectangle (10, 10);
% red circles
% import numpy as np
......
......@@ -2,20 +2,22 @@ from setuptools import setup
with open('README.md') as fh:
long_desc = fh.read()
name = 'cosmic_dustbox'
setup(
name='cosmic_dustbox',
version='0.1.dev',
name=name,
version='0.1.dev0',
author_email='findessp@yandex.ru',
description='Toolbox for cosmic dust.',
long_description=long_desc,
long_description_content_type='text/markdown',
url='https://github.com/findesgh/cosmic_dustbox',
url='https://github.com/findesgh/'+name,
license='GPL-3.0',
packages=['cosmic_dustbox'],
packages=[name],
install_requires=[
'numpy',
'astropy'
'astropy',
'scipy'
],
tests_require=[
'pytest'
......