Commit e1d15de8 authored by Martin Glatzle's avatar Martin Glatzle

Merge branch 'crefin'

parents aee20c82 60846885
from scipy import interpolate as _interp
import astropy.units as _u
import numpy as _np
import os as _os
class Crefin(object):
def __init__(self, f):
self.f = f
return
def __call__(self, a, wave):
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):
"""
Interpolate complex refractive index from given data.
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
Array of the complex refractive index. Note: if `n` is 2-D, its
first index is wavelength and its second index is grain size. For
more details see the documentation of scipy.interpolate.interp2d.
bounds_error : bool, optional
Whether to raise out of bounds errors when interpolating.
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. If
`bounds_error` is False or if `a` is None, this is ignored.
Note
----
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 highly inefficient but convenient. Like this the interpolation
behaves exactly as a 2D function evaluated on a meshgrid of `x` and `y`
would.
"""
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(a, lam, n.real,
bounds_error=bounds_error)
iinterpolation = _interp.interp2d(a, lam, n.imag,
bounds_error=bounds_error)
if bounds_error and not a_bounds_error:
def f(aa, lam):
r = _np.clip(aa, _np.amin(a), _np.amax(a))
xind = _np.argsort(r)
yind = _np.argsort(lam)
return (
rinterpolation(r, lam)
+ 1j*iinterpolation(r, lam))[yind][:, xind].T
else:
def f(a, lam):
xind = _np.argsort(r)
yind = _np.argsort(lam)
return (
rinterpolation(a, lam)
+ 1j*iinterpolation(a, lam))[yind][:, xind].T
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)
n = _np.array(n).T
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.
import astropy.units as _u
import numpy as _numpy
import warnings as _warnings
import os as _os
from scipy.interpolate import interp1d
from scipy import integrate
import sys
Path = _os.path.dirname(_os.path.dirname(_os.path.realpath(__file__))) + \
'//Solver//'
sys.path.insert(0, Path)
import sdist as sd
import GrainSpecies as gs
###############################################################################
###############################################################################
class Dp(object):
def __init__(self, en):
self.en = en
return
def Sigma_abs_E(en):
sizes = _np.logspace(self.sizeMin, self.sizeMax,
num=10000)
def func1(sizes):
if issubclass(self.__class__, sd.SizeDist):
dist_gr = self.__class__(sizes).sdist_gr
dist_sil = self.__class__(sizes).sdist_sil
return dist_gr, dist_sil
def func2(sizes):
sigma_cr = gs.SpGraphite.Sigma_abs(en, sizes)[0]
sigma_sl = gs.SpAsil.Sigma_abs(en, sizes)[0]
return sigma_cr, sigma_sl
f1_gr = interp1d(sizes, func1(sizes)[0], kind='linear')
f2_sl = interp1d(sizes, func1(sizes)[1], kind='linear')
f3_gr = interp1d(sizes, func2(sizes)[0], kind='linear')
f4_sl = interp1d(sizes, func2(sizes)[1], kind='linear')
Integ1 =integrate.quad(f1_gr(Sizes) * f3_gr(Sizes),self.sizeMin,
self.sizeMax)
Integ2 = integrate.quad(f1_gr(Sizes), self.sizeMin,self.sizeMax)
Integ3 = integrate.quad(f2_sl(Sizes) * f4_sl(Sizes),self.sizeMin,
self.sizeMax)
Integ4 = integrate.quad(f2_sl(Sizes), self.sizeMin,self.sizeMax)
sigma_abs_1 = Integ1/Integ2
sigma_abs_2 = Integ3/Integ4
sigma_abs_E = sigma_abs_1 + sigma_abs_2
return sigma_abs_E
import astropy.units as _u
import numpy as _np
import warnings as _warnings
import os as _os
import miepython as miepy
from scipy import interpolate as _interp
import sys
Path = _os.path.dirname(_os.path.dirname(_os.path.realpath(__file__))) + \
'//Solver//'
sys.path.insert(0, Path)
import Solver as sv
###############################################################################
###############################################################################
class GrainSpecies(object):
def __init__(self):
return
class SpAsil(Grainspecies):
def getSilica_props(en):
en = en.to(u.micron, equivalencies=u.spectral())
if en.unit != _u.eV or en.unit != en.J :
raise('The unit of Energy should be in eV or Joule')
elif en.unit == en.J:
en = en * _u.J.to(_u.eV)
lamda = en.unit * _u.eV.to(_u.micron, \
equivalencies=_u.spectral())
path = cls.loadFiles['Silica']
silica = loadOprops(path)
assert lamda.unit == _u.micron, 'The wave length is not in micron'
Silica_d = Find_Interpolate(lamda,silica)
return Silica_d
def Sigma_abs(en,sizes):
m = getSilica_props(en,sizes)
QABS = []
sigma_abs = []
for i, size in enumerate(sizes):
QABS append.(miepy.mie(m,size))
sigma_abs.append.(_np.pi * size**2 * QABS[i])
return _np.array(sigma_abs)
def Sigma_sca(en,sizes):
m = getSilica_props(en,sizes)
QSCA = []
sigma_sca = []
for i, size in enumerate(sizes):
QSCA append.(miepy.mie(m,size))
sigma_sca.append.(_np.pi * size**2 * QSCA[i])
return _np.array(sigma_sca)
class SpGraphite(Grainspecies):
"""
"""
def getCarbon_props(en):
if en.unit != _u.eV or en.unit != en.J :
raise('The unit of Energy should be in eV or Joule')
elif en.unit == en.J:
en = en * _u.J.to(_u.eV)
lamda = en.to(_u.micron, \
equivalencies=_u.spectral())
path1 = cls.loadFiles['carbon_par_0.01']
c_par_0.01_ = loadOprops(path1)
path2 = cls.loadFiles['carbon_par_0.1']
c_par_0.1_ = loadOprops(path1)
path3 = cls.loadFiles['carbon_per_0.01']
c_per_0.01_d = loadOprops(path3)
path4 = cls.loadFiles['carbon_per_0.1']
c_per_0.1_d = loadOprops(path4)
assert lamda.unit == _u.micron, 'The wave length is not in micron'
#
c_par_0.01_d = Find_Interpolate(lamda,c_par_0.01_)
c_par_0.1_d = Find_Interpolate(lamda,c_par_0.1_)
c_per_0.01_d = Find_Interpolate(lamda,c_per_0.01_)
c_par_0.1_d = Find_Interpolate(lamda,c_per_0.1_)
#
return c_par_0.01_d , c_par_0.1_d , c_per_0.01_d , c_per_0.1_d
def Sigma_abs(en,sizes):
c_par_0.01_, c_par_0.1_, c_per_0.01_, c_per_0.1_ =
getCarbon_props(en,sizes)
QABS_par_0.01_ = [] ; QABS_par_0.1_ = [] ; QABS_per_0.01_ = []
QABS_per_0.1_ = []
sigma_par_0.01_ = [] ; sigma_par_0.1_ = [] ; sigma_per_0.01_ = []
sigma_per_0.01_ = []
for i, size in enumerate(sizes):
QABS_par_0.01_ = miepy.mie(c_par_0.01_ ,size)
QABS_par_0.1_ = miepy.mie(c_par_0.1_ ,size)
QABS_per_0.01_ = miepy.mie(c_per_0.01_ ,size)
QABS_per_0.1_ = miepy.mie(c_per_0.1_ ,size)
#
sigma_par_0.01_.append.(_np.pi * size**2 * QABS_par_0.01_[i])
sigma_par_0.1_.append.(_np.pi * size**2 * QABS_par_0.1_[i])
sigma_per_0.01_.append.(_np.pi * size**2 * QABS_per_0.01_[i])
sigma_per_0.1_.append.(_np.pi * size**2 * QABS_per_0.1_[i])
#
sigma_abs_0.01_ = ( 2 * _np.array(sigma_par_0.01_) + \
_np.array(sigma_per_0.01_) ) / 3
sigma_abs_0.1_ = ( 2 * _np.array(sigma_par_0.1_) + \
_np.array(sigma_per_0.1_) ) / 3
#
return sigma_abs_0.01_ , sigma_abs_0.1_
#
def Sigma_sca(en,sizes):
#
c_par_0.01_, c_par_0.1_, c_per_0.01_, c_per_0.1_ =
getCarbon_props(en,sizes)
QSCA_par_0.01_ = [] ; QSCA_par_0.1_ = [] ; QSCA_per_0.01_ = []
QSCA_per_0.1_ = []
sigma_par_0.01_ = [] ; sigma_par_0.1_ = [] ; sigma_per_0.01_ = []
sigma_per_0.01_ = []
for i, size in enumerate(sizes):
QSCA_par_0.01_ = miepy.mie(c_par_0.01_ ,size)
QSCA_par_0.1_ = miepy.mie(c_par_0.1_ ,size)
QSCA_per_0.01_ = miepy.mie(c_per_0.01_ ,size)
QSCA_per_0.1_ = miepy.mie(c_per_0.1_ ,size)
#
sigma_par_0.01_.append.(_np.pi * size**2 * QSCA_par_0.01_[i])
sigma_par_0.1_.append.(_np.pi * size**2 * QSCA_par_0.1_[i])
sigma_per_0.01_.append.(_np.pi * size**2 * QSCA_per_0.01_[i])
sigma_per_0.1_.append.(_np.pi * size**2 * QSCA_per_0.1_[i])
#
sigma_sca_0.01_ = ( 2 * _np.array(sigma_par_0.01_) + \
_np.array(sigma_per_0.01_) ) / 3
sigma_sca_0.1_ = ( 2 * _np.array(sigma_par_0.1_) + \
_np.array(sigma_per_0.1_) ) / 3
#
return sigma_sca_0.01_ , sigma_sca_0.1_
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
###############################################################################
......@@ -164,10 +165,69 @@ 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 Log_Normal(SizeDist):
def __init__(self, sizeMin, sizeMax, rho, sgma, bc, b, 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 +
_np.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 WD01_dst(SizeDist):
def __init__(self, sizeMin, sizeMax, a_t, beta, a_c, alpha, 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
###############################################################################
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.random.random(10) + 1j * np.random.random(10)
cls.a2 = np.array([1e-9, 1e-8])
cls.lam2 = np.linspace(1, 10, num=10)
cls.n2 = (
np.random.random(20) + 1j * np.random.random(20)).reshape(10, 2)
return
@classmethod
def tearDownClass(cls):
return
def test_no_a(self):
c = self.__class__
a = crefin.Crefin.fromData(c.a1, c.lam1, c.n1)
return
def test_normal(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],