Merge branch 'Cosmic_dustbox2' into dev

cosmic_dustbox/crefin.py

+    loadFiles = { 'carbon_par_0.01': 'callindex.out_CpaD03_0.01.txt', 
+                    'carbon_par_0.1':'callindex.out_CpaD03_0.10.txt',
+                    'carbon_per_0.01':'callindex.out_CpeD03_0.01.txt' ,
+                    'carbon_per_0.1': 'callindex.out_CpeD03_0.10.txt',
+                    'silica': 'callindex.out_silD03.txt' }
+
+
+    def Find_Interpolate(lamda,array):
+        w = array[0]
+        if lamda in w:
+            ind = _np.where(w == lamda)
+            Re_n_out = array[1][ind]
+            Im_n_out = array[2][ind]
+        else:    
+            ind1 = _np.where(w >= lamda).min()
+            ind2 = _np.where(w <= lamda).max()
+            wght1 = (w[ind1] - lamda)/lamda
+            wght2= (lamda - w[ind2]) / lamda         
+            Re_n_out = (wght1* array[1][ind1] + wght2* array[1][ind2]) / 2 
+            Im_n_out = (wght1* array[2][ind1] + wght2* array[2][ind2]) / 2
+            Ref_Indx = Re_n_out + j * Im_n_out
+        return Ref_Indx
+
+
+class Crefin(object):
+
+    def __init__(self, f):
+        self.f = f
+        return
+
+    def __call__(self, a, en):
+        return self.f(a, en)
+
+
+class SGPAHCrefin(Crefin):
+
+    @classmethod
+    def loadOprops(path):
+        """
+        Load the optical properties files in CSV format.
+
+        'path': path to Refractive Index and Material properties file
+
+        returns:
+            energy (1D) [eV] at which Q values where computed
+            Real Refractive index and Imaginary Refractive index
+        """
+            rel_path = '\\lib\\' + path
+            path = os.path.dirname(os.path.realpath(__file__)) + rel_path
+
+        data = []
+        with open(path) as f:
+            for line in f:
+                if 'number of wavelengths' not in line:
+                    continue
+                else:
+                    line = f.next().split()
+                    w = []
+                    eps_1 = []
+                    eps_2 = []
+                    real_n = []
+                    Im_n = []
+                    while line:
+                        try:
+                            line = f.next().split()
+                        except StopIteration:
+                            line = []
+                        if line:
+                            w.append(float(line[0]))
+                            Re_n.append(float(line['Re(n)-1']) + 1)
+                            Im_n.append(float(line['Im(n)']))
+                    data.append([
+                        _np.array(w),
+                        _np.array(Re_n),
+                        _np.array(Im_n)])
+
+             return data
+
+
+class Crefin_Asil(SGPAHCrefin):
+
+    def __init__(self):
+        return
+
+
+class Crefin_Gra(SGPAHCrefin):
+
+    def __init__(self):
+        return

cosmic_dustbox/gpop.py

+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

cosmic_dustbox/gspecies.py

+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_

cosmic_dustbox/models/SGPAH/WD01.py

+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

cosmic_dustbox/sdist.py

 ###############################################################################
+import astropy.units as _u
 import numpy as _np
 ###############################################################################
 
@@ -168,6 +169,60 @@ class PowerLaw(SizeDist):
         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(self).__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(self).__init__(sizeMin, sizeMax, f)
+            return
+
+
 ###############################################################################
 if __name__ == "__main__":
     import doctest