Merge pull request #42 from Morisset/devel

1.1.22

docs/CHANGES/1.1.20-1.1.21.txt

@@ -0,0 +1,2 @@
+Correct bug in utils/manage_atomic_data/extract_flt : bytes and str confusion.
+

docs/CHANGES/1.1.21-1.1.22.txt

@@ -0,0 +1,2 @@
+Add Calzetti 2000 extinction laws
+

pyneb/atomic_data/o_iii_atom_TZ17.dat

@@ -18,7 +18,7 @@ Aij
 *** O III transition probabilities
 *** ATOM oxygen
 *** SPECTRUM 3
-*** N_LEVELS 6
+*** N_LEVELS 14
 *** GSCONFIG unknown
 *** SOURCE3 Tayal, S. S. & Zatsarinny, O. 2017, Astrophysical Journal, 850, 147
 *** NOTE3 A-values from all levels

pyneb/atomic_data/o_iii_coll_TZ17.dat

@@ -20305,6 +20305,6 @@
 *** NOTE1    Collision strengths
 *** ATOM oxygen
 *** SPECTRUM 3
-*** N_LEVELS 5
+*** N_LEVELS 202
 *** GSCONFIG p2 
 *** T_UNIT K

pyneb/core/pynebcore.py

@@ -46,7 +46,9 @@
     import h5py
 if config.INSTALLED['ai4neb']:
     from ai4neb import manage_RM
-
+
+
+
 # Change the profiler to 'cpu', 'mem' or None to profile the execution of Atom.
 profiler = None
 #profiler = 'cpu'
@@ -1211,23 +1213,29 @@ class Atom(object):
 
 
     @profile
-    def __init__(self, elem=None, spec=None, atom=None, OmegaInterp='linear', noExtrapol = False, NLevels=None):
+    def __init__(self, elem=None, spec=None, atom=None, OmegaInterp='linear', noExtrapol = False, NLevels=None,
+                 pumpingSED=None):
         """
         Atom constructor
         
         Parameters:
             elem:          symbol of the selected element
             spec:          ionization stage in spectroscopic notation (I = 1, II = 2, etc.)
             atom:          ion (e.g. 'O3').
-            OmegaInterp:   option "kind" from scipy.interpolate.interp1d method: 
+            OmegaInterp [linear]:   option "kind" from scipy.interpolate.interp1d method: 
                             'linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic', 'previous', 'next', 
                             where 'zero', 'slinear', 'quadratic' and 'cubic' refer to a spline interpolation of 
                             zeroth, first, second or third order; 'previous' and 'next' simply return the 
                             previous or next value of the point. 
                             "Cheb" works only for fits files for historical reasons.
-            noExtrapol:    if set to False (default), Omega will be extrapolated above and below
+            noExtrapol [False]:     if set to False (default), Omega will be extrapolated above and below
                             the highest and lowest temperatures where it is defined. If set to True
                             a NaN will be return.
+            NLevels [None]: If specified, set the maximum number of levels the atom will consider. 
+                            The actual value may be lower, depending on the number of levels
+                            available in the coll and atom data.
+            pumpingSED [None]: If specified, fluorescence pumping is performed in the level population estimation.
+                            It is a function of the wavelength (Angstrom) givien the flux as 
             
         **Usage:**
             O3 = pn.Atom('O',3)
@@ -1282,6 +1290,7 @@ def __init__(self, elem=None, spec=None, atom=None, OmegaInterp='linear', noExtr
         self.calling = 'Atom ' + self.atom
         self.log_.message('Making atom object for {0} {1}'.format(self.elem, self.spec), calling=self.calling)
         self.NLevels = NLevels
+        self.pumpingSED = pumpingSED
         dataFile = atomicData.getDataFile(self.atom, data_type='atom')
         if dataFile is None:
             self.atomFileType = None
@@ -1369,6 +1378,11 @@ def __init__(self, elem=None, spec=None, atom=None, OmegaInterp='linear', noExtr
         self.energy_eV = CST.RYD_EV * self.energy_Ryd
 
         self._A = self.getA() # index = quantum number - 1
+        self._B = np.zeros_like(self._A)
+        for i in range(self.atomNLevels): #upper
+            for j in range(i): #lower
+                self._B[i,j] = self._A[i,j] / (8 * np.pi * CST.HPLANCK * (self.getEnergy(i+1, unit='cm-1')-self.getEnergy(j+1, unit='cm-1'))**3) 
+                self._B[j,i] = self._B[i,j] * self.getStatWeight(i+1) / self.getStatWeight(j+1)
         self._Energy = self.getEnergy() # Angstrom^-1
         self._StatWeight = self.getStatWeight()
         if self.NLevels > 0:
@@ -1404,6 +1418,21 @@ def __init__(self, elem=None, spec=None, atom=None, OmegaInterp='linear', noExtr
                                 }
 
 
+    def getB(self, lev_1= -1, lev_2= -1):
+        """
+        Return the value of B(i,j) or B(j,i)
+
+        Parameters:
+            lev_1:  first level
+            lev_2:  end level
+        
+        **Usage:**
+        
+            O3.getB(5,2)
+
+        """
+        return self._B[lev_1-1, lev_2-1]
+
     def getOmega(self, tem, lev_i= -1, lev_j= -1, wave= -1):
         """
         Return interpolated value of the collision strength value at the given temperature 
@@ -1787,11 +1816,18 @@ def getPopulations(self, tem, den, product=True, NLevels=None):
             den_rav = den.ravel()
             q = self.getCollRates(tem_rav, n_level)
             A = self._A[:n_level, :n_level]
+            if self.pumpingSED is None:
+                FB = self._B[:n_level, :n_level] * 0.0
+            else:
+                FB = self._B[:n_level, :n_level] * self.pumpingSED(self.wave_Ang[:n_level, :n_level])
+                for i in range(n_level):
+                    FB[i,i] = 0.0
             pop_result = np.zeros(res_shape_rav1)
             coeff_matrix = np.ones(res_shape_rav2)
             sum_q_up = np.zeros(res_shape_rav1)
             sum_q_down = np.zeros(res_shape_rav1)
             sum_A = A.sum(axis=1)
+            sum_FB = FB.sum(axis=1)
             n_tem = tem_rav.size
             # Following line changed 29/11/2012. It made the code crash when atom_nlevels diff coll_nlevels
             #Atem = np.outer(self._A, np.ones(n_tem)).reshape(n_level, n_level, n_tem)
@@ -1806,13 +1842,13 @@ def getPopulations(self, tem, den, product=True, NLevels=None):
             for row in range(1, n_level):
                 # upper right half            
                 for col in range(row + 1, n_level):
-                    coeff_matrix[row, col] = den_rav * q[col, row] + A[col, row]
+                    coeff_matrix[row, col] = den_rav * q[col, row] + A[col, row] + FB[col, row]
                 # lower left half
                 for col in range(row):
-                    coeff_matrix[row, col] = den_rav * q[col, row]
+                    coeff_matrix[row, col] = den_rav * q[col, row] + FB[col, row]
                 # diagonal
-                coeff_matrix[row, row] = -(den_rav * (sum_q_up[row] + sum_q_down[row]) + sum_A[row])
-
+                coeff_matrix[row, row] = -(den_rav * (sum_q_up[row] + sum_q_down[row]) + sum_A[row] + sum_FB[row])
+                
             vect = np.zeros(n_level)
             vect[0] = 1.
 
@@ -1827,8 +1863,7 @@ def getPopulations(self, tem, den, product=True, NLevels=None):
             pop = np.squeeze(pop_result.reshape(res_shape1))
 
         return pop
-
-
+
     def getLowDensRatio(self, lev_i1=-1, lev_i2=-1, lev_j1=-1, lev_j2=-1, 
                         wave1=-1, wave2=-1, to_eval=None, tem=1e4, lowden = 1e-20):
 

pyneb/extinction/red_corr.py

@@ -77,6 +77,7 @@ def my_X(wave, params = [5000., 1., 2., 3.]):
         #self._laws_dict['F99-like IDL'] = self._F99_like_IDL
         self._laws_dict['F99'] = self._F99
         self._laws_dict['F88 F99 LMC'] = self._FM88_F99_LMC
+        self._laws_dict['Cal00'] = self._Cal00
 
         self.UserFunction = UserFunction
         self.UserParams = None
@@ -332,7 +333,7 @@ def plot(self, w_inf=1000., w_sup=10000., laws=None, ax=None, **kwargs):
         if ax is None:
             f, ax = plt.subplots()
         colors = ['r', 'g', 'b', 'y', 'm', 'c']
-        styles = ['-', '--', '-:', ':']
+        styles = ['-', '--', '-.', ':']
         if not pn.config.INSTALLED['plt']:
             pn.log_.error('matplotlib.pyplot not available for plotting', calling=self.calling)
         old_E_BV = self.E_BV
@@ -918,7 +919,33 @@ def _FM88_F99_LMC(self, wave):
         self.FitzParams = [x0, gamma, c1, c2, c3, c4]
         return self._F99_like(wave)
 
-
+    def _Cal00(self, wave):
+
+        """
+        This function returns the extinction curve proposed by Calzetti et al. 2000 
+        for actively star-forming galaxies.
+
+        reference:
+        The Dust Content and Opacity of Actively Star-forming Galaxies 
+        Calzetti, D., Armus, L., Bohlin, R. C., Kinney, A. L., Koornneef, J., & Storchi-Bergmann, T.
+        2000, Astrophysical Journal, 533, 682-695
+        http://adsabs.harvard.edu/abs/2000ApJ...533..682C
+        """
+
+        x = 1e4 / np.asarray([wave]) # inv microns
+
+        Xx = np.zeros_like(x)
+
+        tt = (x > 1/0.63) & (x <= 1/0.12)
+        res = 2.659 * (((0.011 * x[tt] - 0.198) * x[tt] + 1.509) * x[tt] - 2.156)
+        Xx[tt] = res + self.R_V
+
+        tt = (x > 1/2.2) & (x <= 1/0.63)
+        res = 2.659 * (-1.857 + 1.040 * x[tt])
+        Xx[tt] = res + self.R_V
+
+        return np.squeeze(Xx)
+
     def _zeros(self, wave):
         """
         No correction, return 0.0

pyneb/version.py

@@ -1,2 +1,3 @@
 # PyNeb version
-__version__ = '1.1.21'
+__version__ = '1.1.22'
+