diff --git a/docs/CHANGES/1.1.20-1.1.21.txt b/docs/CHANGES/1.1.20-1.1.21.txt new file mode 100644 index 00000000..47a2b02a --- /dev/null +++ b/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. + diff --git a/docs/CHANGES/1.1.21-1.1.22.txt b/docs/CHANGES/1.1.21-1.1.22.txt new file mode 100644 index 00000000..b57667b2 --- /dev/null +++ b/docs/CHANGES/1.1.21-1.1.22.txt @@ -0,0 +1,2 @@ +Add Calzetti 2000 extinction laws + diff --git a/pyneb/atomic_data/o_iii_atom_TZ17.dat b/pyneb/atomic_data/o_iii_atom_TZ17.dat index 3e34c486..67bf1fe3 100644 --- a/pyneb/atomic_data/o_iii_atom_TZ17.dat +++ b/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 diff --git a/pyneb/atomic_data/o_iii_coll_TZ17.dat b/pyneb/atomic_data/o_iii_coll_TZ17.dat index 6793d862..8e142145 100644 --- a/pyneb/atomic_data/o_iii_coll_TZ17.dat +++ b/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 diff --git a/pyneb/core/pynebcore.py b/pyneb/core/pynebcore.py index 115c4f30..f9aa1d5d 100644 --- a/pyneb/core/pynebcore.py +++ b/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,7 +1213,8 @@ 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 @@ -1219,15 +1222,20 @@ def __init__(self, elem=None, spec=None, atom=None, OmegaInterp='linear', noExtr 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): diff --git a/pyneb/extinction/red_corr.py b/pyneb/extinction/red_corr.py index 5509d3e3..64551147 100644 --- a/pyneb/extinction/red_corr.py +++ b/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 diff --git a/pyneb/version.py b/pyneb/version.py index a058468d..56f4af50 100644 --- a/pyneb/version.py +++ b/pyneb/version.py @@ -1,2 +1,3 @@ # PyNeb version -__version__ = '1.1.21' +__version__ = '1.1.22' +