diff --git a/colradpy/__init__.py b/colradpy/__init__.py index ae0e1e3..1ebebba 100644 --- a/colradpy/__init__.py +++ b/colradpy/__init__.py @@ -7,3 +7,5 @@ from .burgess_tully_rates import * from .split_multiplet import * from .nist_read_txt import * +from .read_FAC import * +from .convolve_EEDF import * \ No newline at end of file diff --git a/colradpy/colradpy_class.py b/colradpy/colradpy_class.py index ff47afe..329e305 100644 --- a/colradpy/colradpy_class.py +++ b/colradpy/colradpy_class.py @@ -34,6 +34,7 @@ sys.path.append('./') from colradpy.r8necip import * from colradpy.read_adf04_py3_class import * +from colradpy.read_FAC import * from colradpy.ecip_rates import * from colradpy.burgess_tully_rates import * from colradpy.split_multiplet import * @@ -160,7 +161,10 @@ def __init__(self,fil,metas=np.array([]),temp_grid=np.array([]),electron_den=np. default_pop_norm=True,temp_dens_pair=False,rate_interp_ion = 'slinear', rate_interp_recomb='log_slinear',rate_interp_col='log_slinear', rate_interp_cx='log_quadratic',use_cx=False,scale_file_ioniz=False, - ne_tau = -1): + ne_tau = -1, + atomic_data_type='adf04', ele='H', nele=0, Zele=1, + EEDF='Maxwellian', atomic_physics='incl_all', + ): """The initializing method. Sets up the nested list for data storage and starts to populate with user data as well as reading in the adf04 file @@ -238,6 +242,14 @@ def __init__(self,fil,metas=np.array([]),temp_grid=np.array([]),electron_den=np. print('Exit here fix to run') sys.exit() + # Atomic data file management + self.data['user']['atomic_data_type'] = atomic_data_type + self.data['user']['ele'] = ele + self.data['user']['nele'] = nele + self.data['user']['Zele'] = Zele + self.data['user']['EEDF'] = EEDF + self.data['user']['atomic_physics'] = atomic_physics + self.populate_data(fil) self.data['atomic']['metas'] = np.asarray(metas) @@ -325,8 +337,22 @@ def populate_data(self,fil): """ if(type(fil) == str or type(fil) == np.str_): - - self.data = self.update_dict(self.data,read_adf04(fil)) + # If atomic data is in adf04-format + if self.data['user']['atomic_data_type'] == 'adf04': + self.data = self.update_dict(self.data,read_adf04(fil)) + # If atomic data is in FAC-format + elif self.data['user']['atomic_data_type'] == 'FAC': + self.data = self.update_dict( + self.data, read_FAC( + ele=self.data['user']['ele'], + nele=self.data['user']['nele'], + Zele=self.data['user']['Zele'], + fil=fil, + EEDF=self.data['user']['EEDF'], + reacts=self.data['user']['atomic_physics'], + Te=self.data['user']['temp_grid'], # [eV] + ) + ) self.data['user']['file_loc'] = str(fil)#make this a string changed for hdfdict else: if(self.data): @@ -860,6 +886,14 @@ def populate_cr_matrix(self): np.sum(np.einsum('ij,j->ij',self.data['rates']['ioniz']['ionization'][i,:,:], self.data['user']['dens_grid']),axis=0) + # Includes autoionization rate if given + if 'autoioniz' in self.data['rates']['ioniz'].keys(): + self.data['cr_matrix']['cr'][i,i,:] -= np.sum( + self.data['rates']['ioniz']['autoioniz'][i,:] + ) + self.data['cr_matrix']['cr_loss'][i,i,:] -= np.sum( + self.data['rates']['ioniz']['autoioniz'][i,:] + ) self.data['cr_matrix']['cr'][i,0:len(self.data['atomic']['energy']),:] = \ @@ -918,6 +952,15 @@ def populate_cr_matrix(self): 0:len(self.data['atomic']['energy']),:] + \ np.einsum('ij,j->ij',self.data['rates']['ioniz']['ionization'][:,p,:], self.data['user']['dens_grid']) + + # Includes autoionization rate if given + if 'autoioniz' in self.data['rates']['ioniz'].keys(): + self.data['cr_matrix']['cr'][len(self.data['atomic']['energy'])+ p, + 0:len(self.data['atomic']['energy']),:] += ( + self.data['rates']['ioniz']['autoioniz'][:,p] + )[:,None] + + if(self.data['user']['use_cx']): for p in range(0,nsigmaplus_cx): @@ -958,6 +1001,16 @@ def populate_cr_matrix(self): self.data['cr_matrix']['cr_loss'][i,i,:,:] = self.data['cr_matrix']['cr'][i,i,:,:] - \ np.sum(np.einsum('ij,k->ijk',self.data['rates']['ioniz']['ionization'][i,:,:], self.data['user']['dens_grid']),axis=0) + + # Includes autoionization rate if given + if 'autoioniz' in self.data['rates']['ioniz'].keys(): + self.data['cr_matrix']['cr'][i,i,:,:] -= np.sum( + self.data['rates']['ioniz']['autoioniz'][i,:] + ) + self.data['cr_matrix']['cr_loss'][i,i,:,:] -= np.sum( + self.data['rates']['ioniz']['autoioniz'][i,:] + ) + #level i populating mechanisms #these are the transition rates from higher levels into the level i self.data['cr_matrix']['cr'][i,0:len(self.data['atomic']['energy']),:,:] = \ @@ -1030,6 +1083,13 @@ def populate_cr_matrix(self): np.einsum('ij,k->ijk',self.data['rates']['ioniz']['ionization'][:,p,:], self.data['user']['dens_grid']) + # Includes autoionization rate if given + if 'autoioniz' in self.data['rates']['ioniz'].keys(): + self.data['cr_matrix']['cr'][len(self.data['atomic']['energy'])+ p, + 0:len(self.data['atomic']['energy']),:,:] += ( + self.data['rates']['ioniz']['autoioniz'][:,p] + )[:,None,None] + @@ -1242,16 +1302,44 @@ def solve_quasi_static(self): self.data['rates']['ioniz']['ionization'][levels_to_keep,:], self.data['processed']['F']) + # Includes autoionization rate if given + if 'autoioniz' in self.data['rates']['ioniz'].keys(): + self.data['processed']['scd'] += np.einsum('ipl,iml->mpl', + np.einsum('ip,l->ipl', + self.data['rates']['ioniz']['autoioniz'][levels_to_keep,:], + 1/self.data['user']['dens_grid'] + ), + self.data['processed']['F'] + ) + if(self.data['processed']['driving_populations_norm']): self.data['processed']['scd'] = self.data['processed']['scd'] + \ np.einsum('ipk,ik->ipk', self.data['rates']['ioniz']['ionization'][self.data['atomic']['metas'],:], 1/(1+np.sum(self.data['processed']['pops_no_norm'][:,self.data['atomic']['metas'],:],axis=0))) + # Includes autoionization rate if given + if 'autoioniz' in self.data['rates']['ioniz'].keys(): + self.data['processed']['scd'] += np.einsum('ipl,il->ipl', + np.einsum('ip,l->ipl', + self.data['rates']['ioniz']['autoioniz'][self.data['atomic']['metas'],:], + 1/self.data['user']['dens_grid'] + ), + 1/(1+np.sum(self.data['processed']['pops_no_norm'][:,self.data['atomic']['metas'],:],axis=0)) + ) + else: self.data['processed']['scd'] = self.data['processed']['scd'] + \ self.data['rates']['ioniz']['ionization'][self.data['atomic']['metas'],:,:] + # Includes autoionization rate if given + if 'autoioniz' in self.data['rates']['ioniz'].keys(): + self.data['processed']['scd'] += ( + np.einsum('ip,l->ipl', + self.data['rates']['ioniz']['autoioniz'][self.data['atomic']['metas'],:], + 1/self.data['user']['dens_grid'] + ) + ) if(self.data['user']['use_recombination'] and self.data['user']['use_recombination_three_body']): #this is the total recombination with three body and the rates that in included in the adf04 file @@ -1373,16 +1461,45 @@ def solve_quasi_static(self): self.data['rates']['ioniz']['ionization'][levels_to_keep,:,:], self.data['processed']['F']) + # Includes autoionization rate if given + if 'autoioniz' in self.data['rates']['ioniz'].keys(): + self.data['processed']['scd'] += np.einsum('ipl,imkl->mpkl', + np.einsum('ip,l->ipl', + self.data['rates']['ioniz']['autoioniz'][levels_to_keep,:], + 1/self.data['user']['dens_grid'] + ), + self.data['processed']['F'] + ) + if(self.data['processed']['driving_populations_norm']): self.data['processed']['scd'] = self.data['processed']['scd'] + \ np.einsum('ipk,ikl->ipkl', self.data['rates']['ioniz']['ionization'][self.data['atomic']['metas'],:,:], 1/(1+np.sum(self.data['processed']['pops_no_norm'][:,self.data['atomic']['metas'],:,:],axis=0))) + # Includes autoionization rate if given + if 'autoioniz' in self.data['rates']['ioniz'].keys(): + self.data['processed']['scd'] += np.einsum('ipl,ikl->ipkl', + np.einsum('ip,l->ipl', + self.data['rates']['ioniz']['autoioniz'][self.data['atomic']['metas'],:], + 1/self.data['user']['dens_grid'] + ), + 1/(1+np.sum(self.data['processed']['pops_no_norm'][:,self.data['atomic']['metas'],:,:],axis=0)) + ) + else: self.data['processed']['scd'] = self.data['processed']['scd'] + \ self.data['rates']['ioniz']['ionization'][self.data['atomic']['metas'],:,:,None] + # Includes autoionization rate if given + if 'autoioniz' in self.data['rates']['ioniz'].keys(): + self.data['processed']['scd'] += ( + np.einsum('ip,l->ipl', + self.data['rates']['ioniz']['autoioniz'][self.data['atomic']['metas'],:], + 1/self.data['user']['dens_grid'] + ) + )[:,:,None,:] + if(self.data['user']['use_recombination'] and self.data['user']['use_recombination_three_body']): #this is the total recombination with three body and the rates that in included in the adf04 file recomb_coeff = np.einsum('ijk,l->ijkl', @@ -1535,6 +1652,16 @@ def solve_time_dependent(self): self.data['processed']['td']['scd'] = np.einsum('ipk,itkl->ptkl',self.data['rates']['ioniz']['ionization'], self.data['processed']['td']['td_pop'][0:len(self.data['rates']['ioniz']['ionization']),:,:,:]) + # Includes autoionization rate if given + if 'autoioniz' in self.data['rates']['ioniz'].keys(): + self.data['processed']['scd'] += np.einsum('ipl,itkl->ptkl', + np.einsum('ip,l->ipl', + self.data['rates']['ioniz']['autoioniz'], + 1/self.data['user']['dens_grid'] + ), + self.data['processed']['td']['td_pop'][0:len(self.data['rates']['ioniz']['ionization']),:,:,:] + ) + def split_pec_multiplet(self): """This function will solve take LS resolved PECs and split them statistically among the j levels @@ -2330,7 +2457,7 @@ def plot_pec_ratio_dens(self,pec1,pec2,temp=np.array([0]),meta = np.array([0]),s - def write_pecs_adf15(self,fil_name='', pec_inds = 0, num = 8, pecs_split=False): + def write_pecs_adf15(self,fil_name='', pec_inds = 0, num = 8, pecs_split=False, wave_lims=None): """ This function calls the write_adf15 function. @@ -2364,7 +2491,8 @@ def write_pecs_adf15(self,fil_name='', pec_inds = 0, num = 8, pecs_split=False): self.data['atomic']['charge_state'], self.data['user']['dens_grid'], self.data['user']['temp_grid'], self.data['atomic']['metas'], self.data['atomic']['ion_pot'], - user = self.data['user'], atomic = self.data['atomic'], num = num) + user = self.data['user'], atomic = self.data['atomic'], num = num, + wave_lims=wave_lims) else:#write un-split PECs @@ -2374,11 +2502,13 @@ def write_pecs_adf15(self,fil_name='', pec_inds = 0, num = 8, pecs_split=False): len(self.data['processed']['wave_air']),dtype='int') write_adf15(fil_name, pec_inds, self.data['processed']['wave_air']*10, - self.data['processed']['pecs'], self.data['atomic']['element'], + self.data['processed']['pecs'], self.data['processed']['pec_levels'], + self.data['atomic']['element'], self.data['atomic']['charge_state'], self.data['user']['dens_grid'], self.data['user']['temp_grid'], self.data['atomic']['metas'], self.data['atomic']['ion_pot'], - user = self.data['user'], atomic = self.data['atomic'], num = num) + user = self.data['user'], atomic = self.data['atomic'], num = num, + wave_lims=wave_lims) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py new file mode 100644 index 0000000..a2d4d9f --- /dev/null +++ b/colradpy/convolve_EEDF.py @@ -0,0 +1,514 @@ +''' + +Given cross-section data, this function calculates rate coefficients +given an electron energy distribution (EEDF) + +cjperks, Jan 18th, 2024 + +''' + +# Modules +import numpy as np +import scipy.constants as cnt +from scipy.special import kn +from scipy.interpolate import interp1d + +############################################################ +# +# Main +# +############################################################ + +def convolve_EEDF( + EEDF = None, # [1/eV] + Te = None, # [eV], dim(ntemp,) + XS = None, # [cm2], dim(nE,ntrans) + engyXS = None, # [eV], dim (nE,ntrans) + m = 0, # 0 or 1 + dE = 0, # [eV] + DC_flag = False, + limit = None, # (optional), high-energy asymptotic behavior + w_upr = None, # (optional) Upper level total angular momentum + w_lwr = None, # (optional) Upper level total angular momentum + ion_L = None, # (optional) ionized state orbital angular momentum + verbose = 1, + use_rel = True, # flag to use relativistic corrections + react = None, # reaction type in ['ce', 'rr', 'ci'] + ): + ''' + INPUTS: + EEDF -- [1/eV], + options: + 1) 'Maxwellian' --> pre-defined maximal entropy distribution function + if Te <10 keV = Maxwell-Boltzmann + if Te >10 keV = Maxwell-Juttner (relativistic) + + 2) 'Gaussian' --> electron beam, i.e. an EBIT + Therefore, Te is interpreted as the central beam voltage, [eV] + NOTE: currently hardcoded a beam width of 5eV. + + 3) NOT IMPLEMENTED YET (non-Maxwellians) + + Te -- [eV], dim(ntemp,), characteristic temperature of EEDF + if using non-Maxwellian then whatever characteristic + dimension that isn't the energy axis + + XS -- [cm2], dim(nE,ntrans), cross-section data + nE -- number of energy grid points + ntrans -- number of transitions considered + NOTE: Assumed to be given on a sufficient grid so + that setting any extrapolated values to zero + causes negligible error + + engyXS -- [eV], dim(nE,ntrans), cross-section energy axis + NOTE: Assumed to be not extremely fine, but fine enough + so that a log-log interpolation causes negligible error + + m -- flag on definition of cross-section energy axis + options: + 1) if m=0, incident electron energy + 2) if m=1, scattered electron energy + + dE -- (optional), [eV], dim(ntrans,), + if m=1 need the bound-bound transition energy difference + + DC_flag -- (optional), if XS is dielectronic capture strength + Assumes cross-section is a delta function in energy + NOTE: units of XS are then [eV *cm2] + + ''' + + # Check + if np.ndim(XS) == 1 and not DC_flag: + XS = XS[:,None] + engyXS = engyXS[:,None] + + # Prepares EEDF + if EEDF == 'Maxwellian': + if react == 'rr': + lwr = 1e-4 + ngrid = int(1e4) + else: + lwr = 1e-3 + ngrid = int(1e3) + EEDF, engyEEDF = _get_Max( + Te=Te, + use_rel=use_rel, + lwr=lwr, + ngrid=ngrid, + ) # dim(ngrid, ntemp), ([1/eV], [eV]) + + elif EEDF == 'Gaussian': + EEDF, engyEEDF = _get_Gauss( + Eb = Te, # [eV] + #dEb = 5, # [eV] + #ngrid = ngrid, + ) # dim(ngrid, ntemp), ([1/eV], [eV]) + + else: + print('NON-MAXWELLIAN ELECTRONS NOT IMPLEMENTED YET') + sys.exit(1) + + # If XS is dielectronic capture strength + if DC_flag: + return _calc_DC( + EEDF=EEDF, + engyEEDF=engyEEDF, + XS=XS, + engyXS=engyXS, + use_rel=use_rel, + ) # [cm3/s], dim(ntemp,ntrans) + + # Performs numerical integration + else: + return _calc_ratec( + EEDF=EEDF, + engyEEDF=engyEEDF, + XS=XS, + engyXS=engyXS, + m=m, + dE=dE, + limit = limit, + w_upr = w_upr, + w_lwr = w_lwr, + ion_L = ion_L, + verbose = verbose, + use_rel=use_rel, + react=react, + ) # [cm3/s], dim(ntemp,ntrans) + + +############################################################ +# +# Utilities +# +############################################################ + +# Calculaes dielectronic capture rate coefficient +def _calc_DC( + EEDF = None, # [1/eV], dim(ngrid,ntemp) + engyEEDF = None, # [eV], dim(ngrid, ntemp) + XS = None, # [eV*cm2], dim(ntrans,) + engyXS = None, # [eV], dim(ntrans,) + use_rel=None, + ): + + # Init output + ratec = np.zeros((EEDF.shape[1], len(XS))) # [cm3/s], dim(ntemp,ntrans) + + # Loop over transitions + for nn in np.arange(len(XS)): + # Incident electron velocity + vel = _get_vel( + E_inc = engyXS[nn], + use_rel = use_rel, + ) # [cm/s], dim(ngrid,) + + # Loop over temperatures + for tt in np.arange(EEDF.shape[1]): + # Interpolates EEDF + EEDF_tmp = 10**interp1d( + np.log10(engyEEDF[:,tt]), + np.log10(EEDF[:,tt]), + bounds_error = False, + fill_value = (-1e5,-1e5) + )(np.log10(engyXS[nn])) # [1/eV] + + # Calculates rate coefficient, [cm3/s] + ratec[tt,nn] = XS[nn] *vel *EEDF_tmp + + # Output, [cm3/s], dim(ntemp,ntrans) + return ratec + +# Calculate rate coefficient integral +def _calc_ratec( + EEDF = None, # [1/eV], dim(ngrid,ntemp) + engyEEDF = None, # [eV], dim(ngrid, ntemp) + XS = None, # [cm2], dim(nE, ntrans) + engyXS = None, # [eV], dim(nE, ntrans), incident electron energy + m = None, + dE = None, # [eV], dim(ntrans,) + limit = None, # dim(ntrans,2) if ce or dim(ntrans,4) if rr, ci + w_upr = None, + w_lwr = None, + ion_L = None, + verbose = None, + use_rel = None, + react = None, + ): + + # Init output + ratec = np.zeros((EEDF.shape[1], XS.shape[1])) # [cm3/s], dim(ntemp,ntrans) + if verbose == 1: + XS_out = np.zeros( + (EEDF.shape[0], EEDF.shape[1], XS.shape[1]) + ) # dim(ngrid, ntemp, ntrans) + engy_out = np.zeros( + (EEDF.shape[0], EEDF.shape[1]) + ) # dim(ngrid, ntemp) + + # Loop over temperatures + for tt in np.arange(EEDF.shape[1]): + # Incident electron velocity + vel = _get_vel( + E_inc = engyEEDF[:,tt], + use_rel = use_rel, + ) # [cm/s], dim(ngrid,) + + # Loop over transitions + for nn in np.arange(XS.shape[1]): + # Interpolates cross-section onto finer grid of EEDF + if m == 0: + engy_tmp = engyXS[:,nn] + elif m == 1: + engy_tmp = engyXS[:,nn] + dE[nn] + + XS_tmp = 10**interp1d( + np.log10(engy_tmp), + np.log10(XS[:,nn]), + bounds_error=False, + fill_value = (-1e5,-1e5) + )(np.log10(engyEEDF[:,tt])) # dim(ngrid,), [cm2] + + # Fill missing data right at threshold + ## NOTE: collisional ioniz quickly drops to 0 near E_inc = dE + ## NOTE: collisional excit is pretty nonlinear but flat enough this should be fine + ## NOTE: radiative recomb blows up near E_inc = dE, should fill with _get_limit + # E_inc is the incident photon energy (photoionization) with min(E_inc) = dE + if react != 'ci': + indE = np.where( + (engyEEDF[:,tt] >= dE[nn]) + & (engyEEDF[:,tt] <= engy_tmp[0]) + )[0] + XS_tmp[indE] = XS[0,nn] + + # Fill values with high-energy asymptotic behavior if available (low-energy for radiative recombination) + if limit is not None: + XS_tmp = _get_limit( + XS_tmp = XS_tmp, + engy_tmp = engy_tmp, + engyEEDF = engyEEDF[:,tt], + dE = dE[nn], + limit = limit[nn,:], + w_upr = w_upr[nn], + w_lwr = w_lwr[nn], + ion_L = ion_L, + react = react, + ) + + # Preforms integration + ratec[tt,nn] = np.trapz( + XS_tmp *vel * EEDF[:,tt], + engyEEDF[:,tt] + ) + + if verbose == 1: + XS_out[:,tt,nn] = XS_tmp + engy_out[:,tt] = engyEEDF[:,tt] + + # Output, [cm3/s], dim(ntemp,ntrans) + if verbose == 1: + return ratec, XS_out, engy_out + else: + return ratec + +############################################################ +# +# Calculation +# +############################################################ + +# Calculate high-energy asymptotic behavior of cross-section +## NOTE: So far, FAC-specific +def _get_limit( + XS_tmp = None, + engy_tmp = None, # [eV], incident electron energy always + engyEEDF = None, # [eV], incident electron energy always + dE = None, # [eV], threshold energy + limit = None, + w_upr = None, + w_lwr = None, + ion_L = None, + react = None, + ): + + # Useful constants + fine_struct2 = cnt.physical_constants['fine-structure constant'][0]**2 # [] + eV2Hartree = 1/cnt.physical_constants['Hartree energy in eV'][0] # [Hartree/eV] + a02 = cnt.physical_constants['Bohr radius'][0]**2 *1e4 # [cm^2] + + # FAC cross-section data isn't really trustworthy about 10x transition energy + tol = 10 + if engy_tmp[-1] >= tol*dE: + indE = np.where(engyEEDF > tol*dE)[0] + else: + indE = np.where(engyEEDF> engy_tmp[-1])[0] + + # If unnecessary + if (len(indE) == 0 and react != 'rr'): + return XS_tmp + + if react == 'ce': + # Electron kinetic momentum squared with fine structure correction + k02 = 2 * engyEEDF[indE] *eV2Hartree *( + 1+ + 0.5 *fine_struct2 *engyEEDF[indE] *eV2Hartree + ) # [atomic units] + + # Bethe asymptotic form of collision strength for (first term) optically + # allowed and (second term) forbidden transitions + omega = ( + limit[0] + * np.log(engyEEDF[indE]/dE) + + limit[1] + ) # [] + + # Cross-section + XS_tmp[indE] = ( + np.pi *omega + / k02 + / (1 +2*w_lwr) + ) * a02 # [cm2] + + elif react == 'rr': + # Note with radiative recombination, you want the low-energy behavior + indE = np.where(engyEEDF < engy_tmp[0])[0] + + # Incident photon energy, [eV] + E_gamma = (engyEEDF[indE] + dE) + + # Photon-electron energy, [eV] + E_e = engyEEDF[indE] + + # Formula for bound-free oscillator strength, [1/Hartree] + xx = (E_e + limit[3]) /limit[3] + yy = (1 +limit[2]) /(np.sqrt(xx) + limit[2]) + + dgf_dE = ( + E_gamma /(E_e +limit[3]) + * limit[0] + * xx**(-3.5 -ion_L +0.5 *limit[1]) + * yy**limit[1] + ) + + eps = E_e *eV2Hartree # [Hartree] + omega = E_gamma *eV2Hartree # [Hartree] + + # Function for photo-ionization cross-section, [atomic units] + XS_PI = ( + 2 *np.pi *np.sqrt(fine_struct2) + /(1 +2*w_lwr) + *(1 +fine_struct2 *eps) + /(1 +0.5 *fine_struct2 *eps) + *dgf_dE + ) + + # Function for radiative recombinarion cross-section, [cm2] + XS_tmp[indE] = ( + fine_struct2/2 + * (1 +2*w_lwr)/(1 +2*w_upr) + * omega**2 + / eps + / (1 +0.5*fine_struct2 *eps) + * XS_PI + ) *a02 + + elif react == 'ci': + # Electron kinetic momentum squared with fine structure correction + k02 = 2 * engyEEDF[indE] *eV2Hartree *( + 1+ + 0.5 * fine_struct2 * engyEEDF[indE] *eV2Hartree + ) # [atomic units] + + # Formula for collision strength + xx = engyEEDF[indE]/dE + yy = 1 - 1/xx + + omega = ( + limit[0] *np.log(xx) + + limit[1] *yy**2 + + limit[2] *1/xx *yy + + limit[3] *1/xx**2 *yy + ) + + # Cross-section + XS_tmp[indE] = ( + omega + / k02 + / (1 +2*w_lwr) + ) * a02 # [cm2] + + # Error check + XS_tmp[XS_tmp<0] = 0 + + # Output + return XS_tmp + +# Calculate Gaussian energy distribution function +def _get_Gauss( + Eb = None, # [eV], dim(nEb,), central beam voltage + dEb = 5, # [eV], (scalar), beam width + ngrid = int(1e3), # energy grid resolution + nstd = 5, # number of standard deviations to make energy grid + ): + + # Calc bounds of energy grid + _dEb = dEb*nstd + + # Energy mesh + engy = np.linspace( + Eb - _dEb, + Eb + _dEb, + ngrid, axis = 0 + ) # dim(ngrid, nEb), [eV] + + # Calculates normalized Gaussian + EEDF = ( + 1/(dEb * np.sqrt(2*np.pi)) + * np.exp( + -(engy-Eb[None,:])**2 + /(2*dEb**2) + ) + ) # dim(ngrid, nEb), [1/eV] + + # Output + return EEDF, engy + +# Calculate Maxwellian energy distribution function +def _get_Max( + Te=None, # [eV], dim(ntemp,) + # Energy grid settings + ngrid = int(1e3), # energy grid resolution + lwr = 1e-3, # energy grid limits wrt multiple of Te + upr = 5e1, # !!!!!! Make sure these make sense for cross-sections + use_rel = None, # Use relativistic corrections + ): + ''' + NOTE: EEDF energy axis defined as incident electron energy + ''' + + # Useful constants + eV2electron = cnt.e /(cnt.m_e*cnt.c**2) # [1/eV] + + # Init output + EEDF = np.zeros((ngrid, len(Te))) # dim(ngrid, ntemp), [1/eV] + engy = np.zeros((ngrid, len(Te))) # dim(ngrid, ntemp), [eV] + + # Loop over temperature + for tt in np.arange(len(Te)): + engy[:,tt] = np.logspace( + np.log10(Te[tt]*lwr), + np.log10(Te[tt]*upr), + ngrid + ) # [eV] + + # relativistic form + if Te[tt] >= 10e3 and use_rel: + # Factors + theta = Te[tt] * eV2electron + gamma = 1 + engy[:,tt] * eV2electron + beta = np.sqrt(1-1/gamma**2) + + EEDF[:,tt] = ( + gamma**2 + * beta + /theta + /kn(2, 1/theta) + *np.exp(-gamma/theta) + ) * eV2electron # dim(ngrid, ntemp); [1/eV] + + # low-energy form + else: + EEDF[:,tt] = ( + 2*np.sqrt(engy[:,tt]/np.pi) + * Te[tt]**(-3/2) + * np.exp(-engy[:,tt]/Te[tt]) + ) # dim(ngrid, ntemp); [1/eV] + + # Output + return EEDF, engy + +# Incident electron velocity +def _get_vel( + E_inc = None, # [eV] + use_rel = None + ): + + # Useful constants + eV2electron = cnt.e /(cnt.m_e*cnt.c**2) # [1/eV] + + # Relativistic form + if use_rel: + return ( + cnt.c *100 *np.sqrt( + 1 - + 1/(E_inc *eV2electron +1)**2 + ) + ) # [cm/s], dim(ngrid,) + # Classical form + else: + return ( + cnt.c *100 *np.sqrt( + 2*E_inc *eV2electron + ) + ) # [cm/s], dim(ngrid,) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py new file mode 100644 index 0000000..1e68bbd --- /dev/null +++ b/colradpy/read_FAC.py @@ -0,0 +1,800 @@ +''' + +read_FAC.py is a subfunction that reads atomic data created by +the Flexible Atomic Code (FAC) to be used in collisional- +radiative modeling in ColRadPy + +NOTE: Assumes you have a local version of FAC + +cjperks +Dec 18, 2023 + +TO DO: + 1) Add non-Maxwellian convolution + 2) Missing energy level data (zpla, zpla1) + 3) Charge exchange data + +''' + +# Module +from pfac import rfac, crm +import sys, os +import numpy as np +import copy +import colradpy.read_FAC_utils as utils +from colradpy.convolve_EEDF import convolve_EEDF + +############################################################ +# +# Main +# +############################################################ + +def read_FAC( + # Ion species of interest because not stored in FAC data files.... + ele = None, # Species name + nele = None, # Number of electrons + Zele = None, # Nuclear charge + # File Management + fil = None, # Common path to FAC files, excluding physics extentions + # Physics controls + EEDF = None, # if None -> assumes Maxwell-averages rates from pfac.fac.MaxwellRate + reacts = None, # if None -> looks for all file suffixes + Te = None, # if not using MaxwellRate files, [eV], dim(ntemp,) + verbose = 0, + ): + + ######## -------- Determines which files to search for -------- ######## + + # Electron energy distribution settings + if EEDF == 'Maxwellian_mr': + # Use Maxwell-averaged rates from pfac.fac.MaxwellRate + use_mr = True + elif EEDF in ['Maxwellian', 'Gaussian']: + use_mr = False + else: + print('REQUESTED ELECTRON ENERGY DISTRIBUTION NOT IMPLEMENTED YET') + sys.exit(1) + + # FAC data files to search for + if reacts == 'incl_all': + reacts = [ + 'en', # Energy levels + 'tr', # Einstein coefficients + 'ce', # Collisional excitation + 'rr', # Radiative recombination + 'ai', # Autoionization/dielectronic recombination + 'ci', # Collision ionization + ] + + ######## -------- Reads data -------- ######## + + # Initialize output + FAC = {} + FAC['rates'] = {} + + # Energy levels + if 'en' in reacts: + FAC = _en( + FAC=FAC, + fil=fil, + ele=ele, + nele=nele, + Zele=Zele, + ) + # Error check + else: + print('NEED TO INCLUDE ENERGY LEVEL DATA IN MODELING!!!') + + # Gets rate coefficient data + # Includes: collisional excit, radiative recomb, collisional ioniz + if use_mr: + # Use Maxwell-averaged rates from pfac.fac.MaxwellRate + FAC = _get_mr( + FAC=FAC, + fil=fil, + verbose=verbose, + reacts=reacts, + ) + else: + # Use cross-section data files + FAC = _get_xs( + FAC=FAC, + fil=fil, + EEDF=EEDF, + Te=Te, + verbose=verbose, + reacts=reacts, + ) + + # Einstein coefficients + if 'tr' in reacts: + FAC = _tr( + FAC=FAC, + fil=fil, + nele=nele, + Zele=Zele, + ) + # Error check + else: + print('NEED TO INCLUDE EINSTEIN COEFFICIENT DATA IN MODELING!!!') + + # Autoionization/dielectronic recombination + if 'ai' in reacts: + FAC = _ai( + nele=nele, + FAC=FAC, + fil=fil, + EEDF=EEDF, + Te=Te, + verbose=verbose, + ) + + # Charge exchange + if 'cx' in reacts: + print('CHARGE EXCHANGE IS NOT IMPLEMENTED YET!!!') + # If empty + else: + FAC['rates']['cx'] = {} + FAC['rates']['cx']['cx_transitions'] = np.asarray([]) + FAC['rates']['cx']['cx_excit'] = np.asarray([]) + + # Infinite energy points + # I assume you made your FAC data setting the appropriate energies + FAC['rates']['inf_engy'] = np.array([], dtype='float64') + + ######## -------- Formats Output -------- ######## + + out = {} + out['atomic'] = copy.deepcopy(FAC['atomic']) + out['input_file'] = copy.deepcopy(FAC) + out['rates'] = copy.deepcopy(FAC['rates']) + + # Output + return out + +############################################################ +# +# Transition Data Files +# +############################################################ + +# Reads energy level file +def _en( + FAC=None, + fil=None, + ele=None, + nele=None, + Zele=None, + ): + + # Useful constants + eV2invcm = 8065.73 # [1/cm/eV] + + # Initializes output + FAC['atomic'] = {} + + # Reads FAC energy levelfile + en = rfac.read_en( + fil+'a.en' + ) + + # Charge states + FAC['atomic']['element'] = ele # Species name + FAC['atomic']['charge_state'] = Zele-nele # Electronic charge + FAC['atomic']['iz0'] = Zele # Nuclear charge + FAC['atomic']['iz1'] = Zele-nele +1 # Spectroscopic charge + + # Initializes data arrays + ion_pot = [] # [1/cm], dim(nion,), ionization potentials + ion_term = [] # dim(nion,), ionization states + ion_pot_lvl = [] # dim(nion,), ionization state index + ion_L = [] # dim(nion,), ionization state orbital angular momentum + + config = [] # dim(ntran,), state configuration + L = [] # dim(ntran,), state L quantum number + S = [] # dim(ntran,), state S quantum number + w = [] # dim(ntran,), state J quantum number + energy = [] # dim(ntran,), state energy level wrt ground + #zpla = [] # dim(ntran,nion), !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + #zpla1 = [] # dim(ntran,nion), !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + # Convert FAC indices to form assumed by ColRadPy + FAC['lvl_indices'] = {} + FAC['lvl_indices']['FAC'] = np.array([], dtype='int') + FAC['lvl_indices']['ColRadPy'] = np.array([], dtype='int') + ind_e = 0 + ind_i = 0 + + # Loop over blocks + for blk in np.arange(len(en[1])): + # If block of excited states + if en[1][blk]['NELE'] == nele: + # Loop over states + for st in np.arange(len(en[1][blk]['ILEV'])): + config.append( + en[1][blk]['sname'][st].decode("utf-8") + ) + + L.append( + en[1][blk]['VNL'][st]%100 + ) + + w.append( + en[1][blk]['2J'][st]/2 + ) + + S.append( + 2*abs(w[-1]-L[-1]) +1 + ) + + energy.append( + en[1][blk]['ENERGY'][st] * eV2invcm + ) + + # Index conversion + FAC['lvl_indices']['FAC'] = np.append( + FAC['lvl_indices']['FAC'], + en[1][blk]['ILEV'][st] + ) + ind_e += 1 + FAC['lvl_indices']['ColRadPy'] = np.append( + FAC['lvl_indices']['ColRadPy'], + ind_e + ) + + + # If block of ionized states + else: + # Loop over states + for st in np.arange(len(en[1][blk]['ILEV'])): + ion_term.append( + en[1][blk]['sname'][st].decode("utf-8") + ) + + ion_pot.append( + en[1][blk]['ENERGY'][st] * eV2invcm + ) + + ion_L.append( + en[1][blk]['VNL'][st]%100 + ) + + # Index conversion + FAC['lvl_indices']['FAC'] = np.append( + FAC['lvl_indices']['FAC'], + en[1][blk]['ILEV'][st] + ) + ind_i += 1 + FAC['lvl_indices']['ColRadPy'] = np.append( + FAC['lvl_indices']['ColRadPy'], + ind_i + ) + + ion_pot_lvl.append( + ind_i + ) + + + # Stores data + FAC['atomic']['ion_pot'] = np.asarray(ion_pot) + FAC['atomic']['ion_term'] = np.asarray(ion_term) + FAC['atomic']['config'] = np.asarray(config) + FAC['atomic']['L'] = np.asarray(L) + FAC['atomic']['S'] = np.asarray(S) + FAC['atomic']['w'] = np.asarray(w) + FAC['atomic']['energy'] = np.asarray(energy) + #FAC['atomic']['zpla'] = np.asarray(zpla) + #FAC['atomic']['zpla1'] = np.asarray(zpla1) + FAC['atomic']['zpla'] = -1*np.ones((len(config), len(ion_pot))) + FAC['atomic']['zpla1'] = -1*np.ones((len(config), len(ion_pot))) + FAC['atomic']['ion_pot_lvl'] = np.asarray(ion_pot_lvl) + FAC['atomic']['ion_L'] = np.asarray(ion_L) + + # Output + return FAC + +# Reads Einstein coefficients +def _tr( + FAC=None, + fil=None, + nele=None, + Zele=None, + ): + + # Init output dictionary + upr = [] # dim(ntrans,) + lwr = [] # dim(ntrans,) + a_val = [] # [1/s], dim(ntrans,) + upr_ColRadPy = [] + lwr_ColRadPy = [] + + # Reads transition rate data file + tr = rfac.read_tr(fil+'a.tr') + + # Loop over blocks + for blk in np.arange(len(tr[1])): + # Loop over transitions + for tran in np.arange(len(tr[1][blk]['lower_index'])): + upr.append( + tr[1][blk]['upper_index'][tran] + ) + + lwr.append( + tr[1][blk]['lower_index'][tran] + ) + + a_val.append( + tr[1][blk]['rate'][tran] + ) + + # Converts FAC indices to ColRadPy notation + ind_upr = np.where(FAC['lvl_indices']['FAC'] == upr[-1])[0][0] + ind_lwr = np.where(FAC['lvl_indices']['FAC'] == lwr[-1])[0][0] + upr_ColRadPy.append( + FAC['lvl_indices']['ColRadPy'][ind_upr] + ) + lwr_ColRadPy.append( + FAC['lvl_indices']['ColRadPy'][ind_lwr] + ) + + # Transition array from TRTable() + trans_ColRadPy = np.vstack((upr_ColRadPy, lwr_ColRadPy)).T # dim(ntrans,2) + + ## NOTE: ColRadPy assumes the index of 'a_val' corresponds to 'col_transitions' + # But, not all transitions have collisional exication or Einstein coefficients + # Therefore need to fill + col_trans = FAC['rates']['excit']['col_transitions'].copy() + col_excit = FAC['rates']['excit']['col_excit'].copy() + all_trans = np.unique( + np.concatenate(( + trans_ColRadPy, + col_trans, + ), axis=0), + axis=0) + all_a_val = np.ones(all_trans.shape[0])*1e-30 # dim(ntrans,) + all_excit = np.ones((all_trans.shape[0], col_excit.shape[1]))*1e-30 # dim(ntrans, ntemp) + + for tt in np.arange(all_trans.shape[0]): + # Finds the indices for this transition + inda = np.where( + np.all( + trans_ColRadPy == all_trans[tt,:], + axis = 1 + ) + )[0] + indc = np.where( + np.all( + col_trans == all_trans[tt,:], + axis = 1 + ) + )[0] + + # Error check + if len(inda) != 1 and len(indc) != 1: + print('TR error') + print(all_trans[tt,:]) + print('xxxx') + + # Fills arrays + if len(inda) == 1: + all_a_val[tt] = a_val[inda[0]] + + if len(indc) == 1: + all_excit[tt,:] = col_excit[indc[0],:] + + # Includes two-photon emission + # Only works for H-like and He-like + if nele <= 2: + a_2photon = crm.TwoPhoton(Zele, nele-1) # [1/s] + + # 2s(0.5) for H-like + if nele == 1: + ind = np.where(FAC['atomic']['config'] == '2s1')[0][0] + # 1s 2s (J=0) for He-like + elif nele == 2: + ind = np.where( + (FAC['atomic']['config'] == '1s1.2s1') + & (FAC['atomic']['w'] == 0) + )[0][0] + + ind2 = np.where( + np.all( + all_trans == np.asarray([ + FAC['lvl_indices']['ColRadPy'][ind], 1 + ]), + axis = 1 + ) + )[0][0] + + all_a_val[ind2] += a_2photon + + # Output + FAC['rates']['a_val'] = all_a_val # [1/s], dim(ntrans,) + FAC['rates']['excit']['col_transitions'] = all_trans # dim(ntrans,2) + FAC['rates']['excit']['col_excit'] = all_excit # [upsilon], dim(ntrans, ntemp) + + # Output + return FAC + +############################################################ +# +# Data Loading & Prep +# +############################################################ + +# Reads dielectronic recombination strength data files +def _ai( + nele=None, + EEDF=None, + Te=None, # [eV], dim(ntemp,) + FAC=None, + fil=None, + verbose=None, + ): + + # Error check if H-like + if nele == 1: + return FAC + + # Init output data + bnd = [] + fre = [] + bnd_ColRadPy = [] + fre_ColRadPy = [] + dE = [] # [eV] + DC = [] # [eV *cm2] + AI = [] # [1/s] + + # Loads FAC AI data file + try: + ai = rfac.read_ai(fil+'a.ai') + except: + print('No AI/DR rates found. Skipping') + return FAC + + # Loop over blocks + for blk in np.arange(len(ai[1])): + # Loop over transitions + for trn in np.arange(len(ai[1][blk]['Delta E'])): + # Stores data + dE.append( + ai[1][blk]['Delta E'][trn] + ) + + DC.append( + ai[1][blk]['DC strength'][trn]*1e-20 + ) + + AI.append( + ai[1][blk]['AI rate'][trn] + ) + + bnd.append( + ai[1][blk]['bound_index'][trn] + ) + fre.append( + ai[1][blk]['free_index'][trn] + ) + + # Converts FAC indices to ColRadPy notation + ind_bnd = np.where(FAC['lvl_indices']['FAC'] == bnd[-1])[0][0] + ind_fre = np.where(FAC['lvl_indices']['FAC'] == fre[-1])[0][0] + bnd_ColRadPy.append( + FAC['lvl_indices']['ColRadPy'][ind_bnd] + ) + fre_ColRadPy.append( + FAC['lvl_indices']['ColRadPy'][ind_fre] + ) + + # Calculates rate coefficient assuming delta function energy-dependent cross-section + ratec = convolve_EEDF( + EEDF=EEDF, + Te=Te, + XS=DC, + engyXS=dE, + DC_flag=True, + use_rel=True, + ).T # [cm3/s],dim(ntrans,ntemp) + + # Transition array from AITable() + trans_ColRadPy = np.vstack((fre_ColRadPy, bnd_ColRadPy)).T # dim(ntrans,2) + + # Combining together autoionization and collisional ionization + ion_trans = FAC['rates']['ioniz']['ion_transitions'].copy() + ion_col = FAC['rates']['ioniz']['ion_excit'].copy() + all_ion_trans = np.unique( + np.concatenate(( + np.flip(trans_ColRadPy,axis=1), + ion_trans, + ), axis=0), + axis=0) # (bound, free) indexing + all_ion_colls = np.ones((all_ion_trans.shape[0], ion_col.shape[1]))*1e-30 # dim(ntrans, ntemp) + all_ion_rates = np.ones(( + len(FAC['atomic']['energy']), len(FAC['atomic']['ion_pot']) + ))*1e-30 # dim(nlvls, nions) + + for tt in np.arange(all_ion_trans.shape[0]): + # Finds the indices for this transition + ind_ai = np.where( + np.all( + np.flip(trans_ColRadPy,axis=1) == all_ion_trans[tt,:], + axis = 1 + ) + )[0] + ind_ci = np.where( + np.all( + ion_trans == all_ion_trans[tt,:], + axis = 1 + ) + )[0] + + # Error check + if len(ind_ai) != 1 and len(ind_ci) != 1: + print('AI Error') + print(all_trans[tt,:]) + print('xxxx') + + # Fills arrays + if len(ind_ci) == 1: + all_ion_colls[tt,:] += ion_col[ind_ci[0],:] + + if len(ind_ai) == 1: + all_ion_rates[trans_ColRadPy[ind_ai,1]-1, trans_ColRadPy[ind_ai,0]-1] += AI[ind_ai[0]] + + # Combining together rad recomb and DC recomb + rad_trans = FAC['rates']['recomb']['recomb_transitions'].copy() + rad_recomb = FAC['rates']['recomb']['recomb_excit'].copy() + all_trans = np.unique( + np.concatenate(( + trans_ColRadPy, + rad_trans, + ), axis=0), + axis=0) + all_recomb = np.ones((all_trans.shape[0], rad_recomb.shape[1]))*1e-30 # dim(ntrans, ntemp) + + for tt in np.arange(all_trans.shape[0]): + # Finds the indices for this transition + ind_dc = np.where( + np.all( + trans_ColRadPy == all_trans[tt,:], + axis = 1 + ) + )[0] + ind_rr = np.where( + np.all( + rad_trans == all_trans[tt,:], + axis = 1 + ) + )[0] + + # Error check + if len(ind_rr) != 1 and len(ind_dc) != 1: + print('DC Error') + print(all_trans[tt,:]) + print('xxxx') + + # Fills arrays + if len(ind_rr) == 1: + all_recomb[tt,:] += rad_recomb[ind_rr[0],:] + + if len(ind_dc) == 1: + all_recomb[tt,:] += ratec[ind_dc[0],:] + + # Output + FAC['rates']['recomb']['recomb_transitions'] = all_trans # dim(ntrans,2) + FAC['rates']['recomb']['recomb_excit'] = all_recomb # [cm3/s], dim(ntrans, ntemp) + + FAC['rates']['ioniz']['ion_transitions'] = all_ion_trans # dim(ntrans, 2) + FAC['rates']['ioniz']['ion_excit'] = all_ion_colls # [reduced], dim(ntrans, 2) + FAC['rates']['ioniz']['autoioniz'] = all_ion_rates # [1/s], dim(nlvls, nions) + + if verbose == 1: + FAC['rates']['recomb']['RR_trans'] = rad_trans # dim(ntrans,2) + FAC['rates']['recomb']['RR_ratec_cm3/s'] = rad_recomb # [cm3/s], dim(ntrans, ntemp) + FAC['rates']['recomb']['DC_trans'] = trans_ColRadPy # dim(ntrans,2) + FAC['rates']['recomb']['DC_ratec_cm3/s'] = ratec # [cm3/s], dim(ntrans,ntemp) + + FAC['rates']['ioniz']['ion_trans_unfil'] = ion_trans # dim(ntrans,2) + + return FAC + +# Reads cross-section data files +def _get_xs( + EEDF=None, + Te=None, # [eV], dim(ntemp,) + FAC=None, + fil=None, + verbose = None, + reacts = None, + ): + + # Error check + if 'ce' not in reacts: + print('NEED TO INCLUDE COLLISIONAL EXCITATION DATA IN MODELING!!!') + + # Fill blank + if 'rr' not in reacts: + FAC['rates']['recomb'] = {} + FAC['rates']['recomb']['recomb_transitions'] = np.asarray([]) + FAC['rates']['recomb']['recomb_excit'] = np.asarray([]) + if 'ci' not in reacts: + FAC['rates']['ioniz'] = {} + FAC['rates']['ioniz']['ion_transitions'] = np.asarray([]) + FAC['rates']['ioniz']['ion_excit'] = np.asarray([]) + + # Useful constants + eV2K = 11604.5 + + # Loop over reactions + for react in reacts: + # Skip + if react not in ['ce', 'rr', 'ci']: + continue + + # Loads date from ascii file + XSdata = utils._read_ascii( + fil = fil, + react = react, + ) + + # Saves temperature grid data + if 'temp_grid' not in FAC.keys(): + # NOTE: ColRadPy assumes Te is in Kelvin within the data file (per ADF04 standard) + FAC['temp_grid'] = Te*eV2K # [K], dim(nt,) + + # Performs EEDF convolution to cross-sections + ( + data, + st_upr, st_lwr, + ratec, XS, engy + ) = utils._conv_ascii2colradpy( + FAC = FAC, + XSdata = XSdata, + EEDF = EEDF, + Te = Te, + react = react, + verbose = verbose, + ) + + # Dictionary labeling + if react == 'ce': + lbl = 'excit' + lbl_trans = 'col_transitions' + lbl_excit = 'col_excit' + elif react == 'rr': + lbl = 'recomb' + lbl_trans = 'recomb_transitions' + lbl_excit = 'recomb_excit' + elif react == 'ci': + lbl = 'ioniz' + lbl_trans = 'ion_transitions' + lbl_excit = 'ion_excit' + + # Formats output + FAC['rates'][lbl] = {} + + # Stores transitions arrays + if react in ['ce', 'rr']: + # If excit -- (upr,lwr) states in ColRadPy indices + # If recomb -- Z+1 state -> Z state + FAC['rates'][lbl][lbl_trans] = np.vstack( + (np.asarray(st_upr), np.asarray(st_lwr)) + ).T # dim(ntrans,2) + else: + # If ioniz -- Z state -> Z+1 state + FAC['rates'][lbl][lbl_trans] = np.vstack( + (np.asarray(st_lwr), np.asarray(st_upr)) + ).T # dim(ntrans,2) + + # Stores rate coefficient data + # If excit -- [upsilon], if recomb -- [cm3/s], if ioniz -- [reduced] + FAC['rates'][lbl][lbl_excit] = np.asarray(data) # dim(ntrans, nt) + + if verbose == 1: + FAC['rates'][lbl]['XS_cm2'] = np.asarray(XS) # dim(ntrans, nE, ntemp), [cm2] + FAC['rates'][lbl]['engy_eV'] = np.asarray(engy) # dim(ntrans, nE, ntemp), [eV] + if react == 'ce': + # Stores transition array before padding to match Einstein coefficient + FAC['rates'][lbl]['col_trans_unfill'] = FAC['rates'][lbl][lbl_trans].copy() + if react != 'rr': + # Stores rate coefficients in SI units + FAC['rates'][lbl]['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, nt), [cm3/s] + + # Output + return FAC + +# Reads Maxwell-averaged data files +def _get_mr( + FAC = None, + fil = None, + verbose = None, + reacts = None, + ): + + # Error check + if 'ce' not in reacts: + print('NEED TO INCLUDE COLLISIONAL EXCITATION DATA IN MODELING!!!') + + # Fill blank + if 'rr' not in reacts: + FAC['rates']['recomb'] = {} + FAC['rates']['recomb']['recomb_transitions'] = np.asarray([]) + FAC['rates']['recomb']['recomb_excit'] = np.asarray([]) + if 'ci' not in reacts: + FAC['rates']['ioniz'] = {} + FAC['rates']['ioniz']['ion_transitions'] = np.asarray([]) + FAC['rates']['ioniz']['ion_excit'] = np.asarray([]) + + # Useful constants + eV2K = 11604.5 + + # Loop over reactions + for react in reacts: + # Skip + if react not in ['ce', 'rr', 'ci']: + continue + + # Reads data file + mr = utils._read_mr( + fil = fil, + react = react, + ) + + # Saves temperature grid data + if 'temp_grid' not in FAC.keys(): + # NOTE: ColRadPy assumes Te is in Kelvin within the data file (per ADF04 standard) + FAC['temp_grid'] = np.asarray(mr['Te_eV'])*eV2K # [K], dim(nt,) + + # Converts data file to ColRadPy form + data, st_upr, st_lwr, ratec = utils._conv_mr2colradpy( + FAC = FAC, + mr = mr, + react = react, + verbose = verbose, + ) + + # Dictionary labeling + if react == 'ce': + lbl = 'excit' + lbl_trans = 'col_transitions' + lbl_excit = 'col_excit' + elif react == 'rr': + lbl = 'recomb' + lbl_trans = 'recomb_transitions' + lbl_excit = 'recomb_excit' + elif react == 'ci': + lbl = 'ioniz' + lbl_trans = 'ion_transitions' + lbl_excit = 'ion_excit' + + # Formats output + FAC['rates'][lbl] = {} + + # Stores transitions arrays + if react in ['ce', 'rr']: + # If excit -- (upr,lwr) states in ColRadPy indices + # If recomb -- Z+1 state -> Z state + FAC['rates'][lbl][lbl_trans] = np.vstack( + (np.asarray(st_upr), np.asarray(st_lwr)) + ).T # dim(ntrans,2) + else: + # If ioniz -- Z state -> Z+1 state + FAC['rates'][lbl][lbl_trans] = np.vstack( + (np.asarray(st_lwr), np.asarray(st_upr)) + ).T # dim(ntrans,2) + + # Stores rate coefficient data + # If excit -- [upsilon], if recomb -- [cm3/s], if ioniz -- [reduced] + FAC['rates'][lbl][lbl_excit] = np.asarray(data) # dim(ntrans, nt) + + if verbose == 1: + if react == 'ce': + # Stores transition array before padding to match Einstein coefficient + FAC['rates'][lbl]['col_trans_unfill'] = FAC['rates'][lbl][lbl_trans].copy() + if react != 'rr': + # Stores rate coefficients in SI units + FAC['rates'][lbl]['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, nt), [cm3/s] + + # Output + return FAC \ No newline at end of file diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py new file mode 100644 index 0000000..bd236d6 --- /dev/null +++ b/colradpy/read_FAC_utils.py @@ -0,0 +1,453 @@ +''' + +read_FAC_utils.py is a script that houses general-purpose +functions for reading FAC data + +cjperks +Jan 24, 2024 + + +''' + +# Module +from pfac import rfac +import os +import numpy as np +import copy +from colradpy.convolve_EEDF import convolve_EEDF +import scipy.constants as cnt + +############################################################ +# +# Converting Data to ColRadPy format +# +############################################################ + +# Converts data read from _read_ascii to form wanted by ColRadPy +# NOTE: Performs EEDF convolution from cross-section to rate coeffs +def _conv_ascii2colradpy( + FAC = None, + XSdata = None, + EEDF = None, + Te = None, + react = None, + verbose = None, + ): + + # Init output + st_lwr = [] + st_upr = [] + data = [] + ratec = [] + XS = [] + engy = [] + + # Loop over lower states + for lwr in XSdata.keys(): + # Converts indices + ind_lwr = np.where(FAC['lvl_indices']['FAC'] == lwr)[0][0] + + # Loop over upper states + for upr in XSdata[lwr].keys(): + # Converts indices + ind_upr = np.where(FAC['lvl_indices']['FAC'] == upr)[0][0] + st_lwr.append( + FAC['lvl_indices']['ColRadPy'][ind_lwr] + ) + st_upr.append( + FAC['lvl_indices']['ColRadPy'][ind_upr] + ) + + # Ionized state oribtal angular momentum for rr cross-section extrapolation + if react == 'rr': + ion_L = FAC['atomic']['ion_L'][st_upr[-1] -1] + else: + ion_L = 0 + + # Calculates Rate coefficient data, [cm3/s], dim(ntrans, ntemp) + if verbose == 1: + ( + tmp_data, + tmp_XS, + tmp_engy + )= convolve_EEDF( + EEDF = EEDF, + Te=Te, + XS = XSdata[lwr][upr]['XS'], + engyXS = XSdata[lwr][upr]['engy'], + m = 0, + dE = np.asarray([XSdata[lwr][upr]['dE']]), + limit = XSdata[lwr][upr]['limit'][None,:], + w_upr = np.asarray([XSdata[lwr][upr]['w_upr']]), + w_lwr = np.asarray([XSdata[lwr][upr]['w_lwr']]), + ion_L = ion_L, + verbose=verbose, + use_rel = True, + react = react, + ) + data.append(tmp_data[:,0]) + XS.append(tmp_XS[:,:,0]) + engy.append(tmp_engy) + + else: + data.append( + convolve_EEDF( + EEDF = EEDF, + Te=Te, + XS = XSdata[lwr][upr]['XS'], + engyXS = XSdata[lwr][upr]['engy'], + m = 0, + dE = np.asarray([XSdata[lwr][upr]['dE']]), + limit = XSdata[lwr][upr]['limit'][None,:], + w_upr = np.asarray([XSdata[lwr][upr]['w_upr']]), + w_lwr = np.asarray([XSdata[lwr][upr]['w_lwr']]), + ion_L = ion_L, + verbose=verbose, + use_rel = True, + react = react, + )[:,0] + ) + + if verbose == 1: + ratec.append(data[-1].copy()) + + if react == 'ce': + # Convert to Upsilon form + data[-1] = _conv_rate2upsilon( + data = data[-1], + Te_eV = Te, + ind_lwr = st_lwr[-1] -1, + ind_upr = st_upr[-1] -1, + FAC = FAC, + ) + + elif react == 'ci': + # Converts ionization data to adf04 reduced form + data_tmp = _conv_rate2reduct( + data = data[-1], + Te_eV = Te, + ind_st = st_lwr[-1] -1, + ind_ion = st_upr[-1] -1, + FAC = FAC, + ) + data[-1] = [x if not np.isnan(x) else 0.0 for x in data_tmp] # fixes overflow error + + # Output + return data, st_upr, st_lwr, ratec, XS, engy + +# Converts data read from _read_mr to form wanted by ColRadPy +def _conv_mr2colradpy( + FAC = None, + mr = None, + react = None, + verbose = None, + ): + + # Init output + st_lwr = [] + st_upr = [] + data = [] + ratec = [] + + if react == 'ce': + lbl = 'coll_excit' + elif react == 'rr': + lbl = 'rad_recomb' + elif react == 'ci': + lbl = 'coll_ion' + + # Loop over states with charge Z + for lwr in mr.keys(): + # Skips temp grid + if lwr == 'Te_eV': + continue + + # Converts indices + ind_lwr = np.where(FAC['lvl_indices']['FAC'] == lwr)[0][0] + + # Loop over states with charge Z+1 + for upr in mr[lwr][lbl].keys(): + # Converts indices + ind_upr = np.where(FAC['lvl_indices']['FAC'] == upr)[0][0] + st_lwr.append( + FAC['lvl_indices']['ColRadPy'][ind_lwr] + ) + st_upr.append( + FAC['lvl_indices']['ColRadPy'][ind_upr] + ) + + data.append( + mr[lwr][lbl][upr] + ) # [cm3/s] + + if verbose == 1: + ratec.append(data[-1].copy()) + + if react == 'ce': + # Convert to Upsilon form + data[-1] = _conv_rate2upsilon( + data = data[-1], + Te_eV = mr['Te_eV'], + ind_lwr = st_lwr[-1] -1, + ind_upr = st_upr[-1] -1, + FAC = FAC, + ) + + elif react == 'ci': + # Converts ionization data to adf04 reduced form + data_tmp = _conv_rate2reduct( + data = data[-1], + Te_eV = mr['Te_eV'], + ind_st = st_lwr[-1] -1, + ind_ion = st_upr[-1] -1, + FAC = FAC, + ) + data[-1] = [x if not np.isnan(x) else 0.0 for x in data_tmp] # fixes overflow error + + # Output + return data, st_upr, st_lwr, ratec + +############################################################ +# +# Reformating Rate Coefficient Data +# +############################################################ + + +# Converts rate coefficeint to adf04 reduced form +def _conv_rate2reduct( + data = None, # [cm3/s], dim(ntemp,), list + Te_eV = None, # [eV], dim(ntemp,) + ind_st = None, + ind_ion = None, + FAC = None, + ): + + # Useful constants + eV2invcm = 8065.73 # [1/cm/eV] + + # Init output + out = [] + + for tt in np.arange(len(Te_eV)): + # Calculates reducted rate + out.append( + data[tt] * np.exp( + (FAC['atomic']['ion_pot'][ind_ion] - FAC['atomic']['energy'][ind_st]) + / eV2invcm + / Te_eV[tt] + ) + ) + + # Output + return out + +# Converts rate coefficient to adf04 Upsilon form +def _conv_rate2upsilon( + data = None, # [cm3/s], dim(ntemp,), array + Te_eV = None, # [eV], dim(ntemp,) + ind_lwr = None, + ind_upr = None, + FAC = None, + ): + ''' + Ref -- A. Burgess and J.A. Tully, Astronomy and Astrophysics, Vol.254, NO. FEB(I), P. 436, 1992 + Eqn 20 + + ''' + + # Useful constants + eV2invcm = 8065.73 # [1/cm/eV] + factor = ( + 2*np.sqrt(np.pi) + *cnt.physical_constants['Bohr radius'][0]*1e2 # [cm] + *cnt.hbar/cnt.m_e *1e4 # [cm^2/s] + ) # [cm^3/s], 2.1716e-8 + + tmp = data * ( + np.sqrt(np.asarray(Te_eV)/13.6058) + /factor + *(1 + 2*FAC['atomic']['w'][ind_lwr]) + * np.exp( + abs( + FAC['atomic']['energy'][ind_upr] + -FAC['atomic']['energy'][ind_lwr] + ) + /(np.asarray(Te_eV)*eV2invcm) + ) + ) + + tmp[np.isnan(tmp)] = 0 + return tmp + +############################################################ +# +# Reading Data Files +# +############################################################ + +# Reads cross-section data files +def _read_ascii( + fil = None, + react = None, + ): + + # Reads data file + if react == 'ce': + data_fil = rfac.read_ce(fil+'a.ce') + lwr_lbl = 'lower_index' + upr_lbl = 'upper_index' + data_lbl = 'crosssection' + elif react == 'rr': + data_fil = rfac.read_rr(fil+'a.rr') + lwr_lbl = 'bound_index' + upr_lbl = 'free_index' + data_lbl = 'RR crosssection' + elif react == 'ci': + data_fil = rfac.read_ci(fil+'a.ci') + lwr_lbl = 'bound_index' + upr_lbl = 'free_index' + data_lbl = 'crosssection' + + # Initializes output dictionary + out = {} + + # Loop over blocks + for blk in np.arange(len(data_fil[1])): + # Loop over transitions + for trn in np.arange(len(data_fil[1][blk]['Delta E'])): + # Transition indices + lwr = data_fil[1][blk][lwr_lbl][trn] + upr = data_fil[1][blk][upr_lbl][trn] + + # If new lower state + if lwr not in out.keys(): + out[lwr] = {} + out[lwr][upr] = {} + + # Stores cross-section data, [cm2] + out[lwr][upr]['XS'] = data_fil[1][blk][data_lbl][trn,:]*1e-20 + + # Stores transition energy, [eV] + out[lwr][upr]['dE'] = data_fil[1][blk]['Delta E'][trn] + if react == 'ce': + out[lwr][upr]['dE'] *= 10**int(np.log10(data_fil[1][blk]['TE0'])) + + if react == 'ce': + # Stores energy grid in terms of incident electron energy, [eV] + if data_fil[1][blk]['ETYPE'] == 0: + out[lwr][upr]['engy'] = data_fil[1][blk]['EGRID'] + elif data_fil[1][blk]['ETYPE'] == 1: + out[lwr][upr]['engy'] = data_fil[1][blk]['EGRID'] + data_fil[1][blk]['TE0'] + + # Stores total angular momentum + out[lwr][upr]['w_upr'] = data_fil[1][blk]['upper_2J'][trn]/2 + out[lwr][upr]['w_lwr'] = data_fil[1][blk]['lower_2J'][trn]/2 + + # Stores parameters for high-energy behavior + out[lwr][upr]['limit'] = np.asarray([ + data_fil[1][blk]['bethe'][trn], + data_fil[1][blk]['born'][trn,0] + ]) + + elif react == 'rr': + # Stores energy grid in terms of photo-electron energy, [eV] + if data_fil[1][blk]['ETYPE'] == 0: + out[lwr][upr]['engy'] = data_fil[1][blk]['EGRID'] - data_fil[1][blk]['Delta E'][trn] + elif data_fil[1][blk]['ETYPE'] == 1: + out[lwr][upr]['engy'] = data_fil[1][blk]['EGRID'] + + # Stores total angular momentum + out[lwr][upr]['w_upr'] = data_fil[1][blk]['free_2J'][trn]/2 + out[lwr][upr]['w_lwr'] = data_fil[1][blk]['bound_2J'][trn]/2 + + # Stores parameters for high-energy behavior, dim(4,) + out[lwr][upr]['limit'] = data_fil[1][blk]['parameters'][trn,:] + + elif react == 'ci': + # Stores energy grid in terms of incident electron energy, [eV] + if data_fil[1][blk]['ETYPE'] == 0: + out[lwr][upr]['engy'] = data_fil[1][blk]['EGRID'] + elif data_fil[1][blk]['ETYPE'] == 1: + out[lwr][upr]['engy'] = data_fil[1][blk]['EGRID'] + data_fil[1][blk]['Delta E'][trn] + + # Stores total angular momentum + out[lwr][upr]['w_upr'] = data_fil[1][blk]['free_2J'][trn]/2 + out[lwr][upr]['w_lwr'] = data_fil[1][blk]['bound_2J'][trn]/2 + + # Stores parameters for high-energy behavior, dim(4,) + out[lwr][upr]['limit'] = data_fil[1][blk]['parameters'][trn,:] + + # Output + return out + +# Read Maxwellian-averaged data files +def _read_mr( + fil = None, + react = None, + ): + + # Reads data file + f = open( + fil+data+'.mr', + 'r' + ) + + # Initializes output dictionary + out = {} + + # Data orginaization labels + if react == 'ce': + label = 'coll_excit' + elif react == 'rr': + label = 'rad_recomb' + elif react == 'ci': + label = 'coll_ion' + + # Loop over lines + for line in f: + # Skip line breaks + if line == '\n': + continue + + # If reading a header + if line.split(' ')[0] == '#': + # Lower level + lwr = int(line.split('\t')[0][1:]) + # Upper level + upr = int(line.split('\t')[2]) + + # Number of temperature points + ntemp = int(line.split('\t')[5]) + indt = 0 + + # If temperature mesh not yet included + if 'Te_eV' not in out.keys(): + out['Te_eV'] = [] + + # If new lower state + if lwr not in out.keys(): + out[lwr] = {} + out[lwr][label] = {} + + out[lwr][label][upr] = [] + + # If reading data + else: + line = line.replace('\t', ' ') + # If need to add temperature mesh + if len(out['Te_eV']) < ntemp: + out['Te_eV'].append( + float(line.split(' ')[1]) + ) + + # Adds rate coefficient data, [cm3/s] + out[lwr][label][upr].append( + float(line.replace(' ', ' ').split(' ')[-1])*1e-10 + ) + + # Increases index to next temp point + indt += 1 + + # Output + return out diff --git a/colradpy/write_adf15.py b/colradpy/write_adf15.py index dcbc270..e8b7783 100644 --- a/colradpy/write_adf15.py +++ b/colradpy/write_adf15.py @@ -1,11 +1,11 @@ import numpy as np import itertools - +from datetime import date # group subset of an array # this is so there are 8 entries per line # in the adf15 file -# see stackexchange +# see stackexchangeG def grouper(n, iterable): it = iter(iterable) while True: @@ -15,8 +15,9 @@ def grouper(n, iterable): yield chunk -def write_adf15(fil_name,pec_inds, wave, pecs, element_sym, - charge_state, dens, temp, metas, ion_pot, user='', atomic='', num=8): +def write_adf15(fil_name,pec_inds, wave, pecs, pec_lvls, element_sym, + charge_state, dens, temp, metas, ion_pot, user='', atomic='', num=8, + wave_lims=None): """Write all or a subset of PECs to an adf15 formatted file. standard adf15 units wave in [A] @@ -68,16 +69,40 @@ def write_adf15(fil_name,pec_inds, wave, pecs, element_sym, f = open(fil_name,'w+')#open the file - pecs = pecs[pec_inds]#only take PECs requested - wave = wave[pec_inds]#only take wavelengths requested + pecs[pecs<=1e-40] = 1e-40 # sets effective zero + + # Filter PECs by strength + tol = 1e-17 + pec_inds_fil = np.array([], dtype='int') + for pp in pec_inds: + if np.mean(pecs[pp, :,:,:]) >= tol: + pec_inds_fil = np.append(pec_inds_fil, pp) + + # Further filters PECs by wavelength + if wave_lims is not None: + pec_inds_fil2 = np.where( + (wave >= wave_lims[0]) + & (wave <= wave_lims[1]) + )[0] + + # Takes intersection + pec_inds_fil = list(set(pec_inds_fil) & set(pec_inds_fil2)) + + pecs = pecs[pec_inds_fil]#only take PECs requested, dim(ntrans,nmeta,ntemp,ndens) + wave = wave[pec_inds_fil]#only take wavelengths requested + pec_lvls = pec_lvls[pec_inds_fil,:]#only take PEC levels requested #first line in the adf15 file format # nothing really important here as its all redunant information later in file - first_line = ' ' + str(len(wave)) + ' /' + element_sym + \ + first_line = ' ' + str(len(pec_inds_fil)*pecs.shape[1]) + ' /' + element_sym + \ ' ' + str(charge_state) + ' PHOTON EMISSIVITY COEFFICIENTS/\n' f.write(first_line) + # Array to store transition info + trans = [] + isel = 0 + for j in range(0,np.shape(pecs)[1]):#loop over metastables for i in range(0, len(wave)):#loop over PECs #letting the user know if the PEC is excitation or recombination @@ -88,11 +113,38 @@ def write_adf15(fil_name,pec_inds, wave, pecs, element_sym, else: typ = 'CHEXC' + # Increases isel index + isel += 1 + + # Transition infor + upr = pec_lvls[i,0] + lwr = pec_lvls[i,1] + trans.append( + 'C'+str(isel).rjust(6, ' ') + + "{:1.5F}".format(wave[i]).rjust(15, ' ') + + ( + str(upr) + +'('+str(atomic['S'][upr])+')' + +str(atomic['L'][upr]) + +'('+str(atomic['w'][upr])+')' + ).rjust(17, ' ') + +'-' + + ( + str(lwr) + +'('+str(atomic['S'][lwr])+')' + +str(atomic['L'][lwr]) + +'('+str(atomic['w'][lwr])+')' + ).rjust(15, ' ') + + typ.rjust(8, ' ') + + '\n' + ) + #header information for the PEC, important things here are the wavelength # number of temps and dens and metastable level - pec_header = ' ' + "{0:.1f}".format(wave[i]*10) + ' A ' + str(len(dens)) +\ + pec_header = ' ' + "{0:.5f}".format(wave[i]) + ' A ' + str(len(dens)) +\ ' ' + str(len(temp)) + ' /FILMEM = bottom /TYPE = ' +typ + ' /INDM = ' +\ - str(j) + ' ISEL = ' + str(i+1)+'\n' + 'T'+ ' /ISEL = ' + str(isel)+'\n' + #str(j) + ' /ISEL = ' + str(isel)+'\n' f.write(pec_header) @@ -121,31 +173,58 @@ def write_adf15(fil_name,pec_inds, wave, pecs, element_sym, f.write('C -----------\n') f.write('C\n') f.write('C ***MADE WITH ColRadPy***\n') - f.write('C NUCLEAR CHARGE = 4\n') - f.write('C ION CHARGE +1 = 3\n') + f.write('C NUCLEAR CHARGE = '+str(element_sym)+'\n') + f.write('C ION CHARGE +1 = '+str(charge_state+1)+'\n') f.write('C\n') - f.write('C SPECIFIC ION FILE :\n') + f.write('C SPECIFIC ION FILE : '+user['file_loc']+'\n') f.write('C\n') + f.write('C TABULATION : photon emissivity coefft (te,ne)\n') + f.write('C UNITS : ph. emis coef[cm^3 s^-1]; te [eV]; ne [cm^-3]\n') f.write('C\n') if(type(user) == dict): for i in user.keys(): f.write('C '+ i + ' ' +str(user[i]) + "\n") if(type(atomic) == dict): - f.write('C ' + 'lvl' +' '+ 'config' + - ' ' + 'S'+ - ' ' + 'L'+ - ' ' + 'w'+ - ' ' + 'eng\n') + f.write('C\n') + f.write('C' + 'lvl'.rjust(5, ' ') +'config'.rjust(20,' ') + + 'S'.rjust(5,' ')+ + 'L'.rjust(5,' ')+ + 'w'.rjust(5,' ')+ + 'energy [cm^-1]\n'.rjust(18, ' ')) f.write('C-----------------------------------------------------------------------\n') f.write('C\n') for ii in range(0,len(atomic['config'])): - f.write('C ' + str(ii) +' '+ str(atomic['config'][ii])+\ - ' ' + str(atomic['S'][ii]) + \ - ' ' + str(atomic['L'][ii]) + \ - ' ' + str(atomic['w'][ii]) + \ - ' ' + str(atomic['energy'][ii])+'\n') + f.write('C' + str(ii).rjust(5, ' ') + str(atomic['config'][ii]).rjust(20, ' ')+\ + str(atomic['S'][ii]).rjust(5, ' ') + \ + str(atomic['L'][ii]).rjust(5, ' ') + \ + str(atomic['w'][ii]).rjust(5, ' ') + \ + "{:0.2F}".format(atomic['energy'][ii]).rjust(18, ' ')+'\n') + + # Transitions header + f.write('C\n') + f.write('C-----------------------------------------------------------------------\n') + f.write('C\n') + f.write( + "C ISEL WVLEN [A] TRANSITION lvl(S)L(w) TYPE\n" + ) + f.write( + "C ----- ------------ -------------------------------- -----\n" + ) + for tt in np.arange(len(trans)): + f.write(trans[tt]) + + # Disclaimer + today = date.today().strftime("%b-%d-%Y") + f.write('C\n') + f.write('C-----------------------------------------------------------------------\n') + f.write('C\n') + f.write('C\tPRODUCER:\t Conor Perks, cjperks@psfc.mit.edu\n') + f.write('C\tDATE: \t'+ today) + f.write('\n') + f.write('C\n') + f.write('C-----------------------------------------------------------------------\n') f.close()