...
 
Commits (4)
###############################################################################
import os as _os
###############################################################################
data_dir = _os.path.join(
_os.path.dirname(__file__), 'data'
......
......@@ -136,7 +136,10 @@ class Crefin(object):
class SGPAHCrefin(Crefin):
"""
Crefin subclass specific to the SGPAH model and the data format in which
the complex refractive indices for it are published.
"""
@classmethod
def parseCrefinFile(cls, path):
with open(path) as f:
......
###############################################################################
import numpy as _np
import astropy.units as u
import astropy.units as _u
import warnings as _warnings
###############################################################################
......@@ -12,11 +13,7 @@ class GrainSpecies(object):
return
def scatt(self, sizes, wave, **kwargs):
return self._scatt(
sizes.to(_u.micron),
wave.to(_u.micron, equivalencies=_u.spectral()),
**kwargs
)
return self._scatt(sizes, wave, **kwargs)
class SphericalGS(GrainSpecies):
......@@ -59,10 +56,13 @@ class IsoHomSphericalGS(HomSphericalGS):
"""
def __init__(self, crefin, mieSolver):
self.crefin = crefin
self.mieSolver = mieSolver
return
def _scatt(self, sizes, wave):
return
def _scatt(self, sizes, wave, **kwargs):
if kwargs:
_warnings.warn("Kwargs provided, these are not used, however!")
return self.mieSolver(sizes, wave, self.crefin(sizes, wave))
class cAxisSphericalGS(HomSphericalGS):
......@@ -80,17 +80,30 @@ class cAxisSphericalGS(HomSphericalGS):
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]_.
([1]_).
References
----------
.. [1] Draine & Malhotra 1993ApJ...414..632D
"""
def __init__(self, crefins, axesWeights, mieSolver):
self.crefins = crefins
super().__init__()
self.axesWeights = axesWeights
self.solver = mieSolver
self.crefins = crefins
self.mieSolver = mieSolver
return
def _scatt(self, sizes, wave):
return
def _scatt(self, sizes, wave, **kwargs):
if kwargs:
_warnings.warn("Kwargs provided, these are not used, however!")
# TODO: very nasty implementation, there's surely a nicer way
gen = (
self.mieSolver(sizes, wave, x(sizes, wave)) for x in self.crefins
)
sabs = sscat = g = 0
for weight, scatt in zip(self.axesWeights, gen):
sabs += weight*scatt[0]
sscat += weight*scatt[1]
g += weight*scatt[2]
return sabs, sscat, g
###############################################################################
from cosmic_dustbox.crefin import SGPAHCrefin as _SGPAHCrefin
from cosmic_dustbox.config import data_dir as _data_dir
import os as _os
from libpymie.miewrapper import MieWrapper as _MieWrapper
from cosmic_dustbox.config import data_dir as _data_dir
from cosmic_dustbox.crefin import SGPAHCrefin as _SGPAHCrefin
from cosmic_dustbox.gspecies import cAxisSphericalGS as _cAxisSphericalGS
from cosmic_dustbox.gspecies import IsoHomSphericalGS as _IsoHomSphericalGS
###############################################################################
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')
])
}
def crefins(name):
"""
Create SGPAHCrefin instances from name.
This function creates SGPAHCrefin instances that correspond to the complex
refractive indices published alongside [1]_ and [2]_.
Parameters
----------
name : string
Identifier. Can be 'cpa' for graphite parallel to the optical axis,
'cpe' for graphite perpendicular to the optical axis or 'sil' for
astronomical silicate.
Returns
-------
c : SGPAHCrefin
Complex refractive index object.
References
----------
.. [1] Draine 2003ApJ...598.1017D
.. [2] Draine 2003ApJ...598.1026D
"""
if name == 'cpa':
c = _SGPAHCrefin.fromFiles([
_os.path.join(_data_dir,
'crefin/callindex.out_CpaD03_0.01'),
_os.path.join(_data_dir,
'crefin/callindex.out_CpaD03_0.10')
])
elif name == 'cpe':
c = _SGPAHCrefin.fromFiles([
_os.path.join(_data_dir,
'crefin/callindex.out_CpeD03_0.01'),
_os.path.join(_data_dir,
'crefin/callindex.out_CpeD03_0.10')
])
elif name == 'sil':
c = _SGPAHCrefin.fromFiles([
_os.path.join(_data_dir,
'crefin/callindex.out_silD03')
])
else:
raise ValueError("Unknown name for crefin: "+str(name))
return c
def gspecies(name):
"""
Create GrainSpecies instances from name.
"""
if name == 'sil':
g = _IsoHomSphericalGS(
crefins(name), _MieWrapper(implementation='cloudy')
)
elif name == 'gra':
g = _cAxisSphericalGS(
[crefins('cpa'), crefins('cpe')],
[1/3, 2/3],
_MieWrapper(implementation='cloudy')
)
else:
raise ValueError("Unknown name for grain species: "+str(name))
return g
###############################################################################
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):
......
###############################################################################
from unittest import TestCase
from ddt import ddt, data
from .. import D03
import numpy as np
import os
from libpymie.tests.test_miewrapper import data_path
import astropy.units as u
###############################################################################
@ddt
class Testsdists(TestCase):
@data(
'cpa',
'cpe',
'sil',
)
def test_init(self, name):
D03.crefins(name)
return
@ddt
class Testgspecies(TestCase):
@data(
'sil',
'gra',
)
def test_init(self, name):
D03.gspecies(name)
return
def test_Fig241(self):
cs_sil_ref = np.load(os.path.join(data_path, 'cabs_sil.npy'))*u.cm**2
cs_gra_ref = np.load(os.path.join(data_path, 'cabs_gra.npy'))*u.cm**2
inv_wave = np.logspace(-3, 2, num=1000)/u.micron
sizes = np.array([0.01, 0.1, 1.0])*u.micron
cs_sil = D03.gspecies('sil').scatt(sizes, inv_wave)[0]
cs_gra = D03.gspecies('gra').scatt(sizes, inv_wave)[0]
self.assertTrue(
np.allclose(cs_sil, cs_sil_ref, atol=0, rtol=1e-7)
)
self.assertTrue(
np.allclose(cs_gra, cs_gra_ref, atol=0, rtol=1e-7)
)
return
###############################################################################
from unittest import TestCase
from .. import WD01
import numpy as np
import astropy.units as u
###############################################################################
class Testsdists(TestCase):
......
......@@ -15,7 +15,8 @@ setup(
install_requires=[
'numpy',
'astropy',
'scipy'
'scipy',
'libpymie'
],
tests_require=[
'pytest'
......