From d7dba63e8d8a751d04d87c2ce95ec1e97590db7f Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 18 Dec 2023 14:12:29 -0500 Subject: [PATCH 01/46] Init function to read FAC data --- colradpy/__init__.py | 1 + colradpy/read_FAC.py | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+) create mode 100644 colradpy/read_FAC.py diff --git a/colradpy/__init__.py b/colradpy/__init__.py index ae0e1e3..6ac38b6 100644 --- a/colradpy/__init__.py +++ b/colradpy/__init__.py @@ -7,3 +7,4 @@ from .burgess_tully_rates import * from .split_multiplet import * from .nist_read_txt import * +from .read_FAC import * \ No newline at end of file diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py new file mode 100644 index 0000000..0291e11 --- /dev/null +++ b/colradpy/read_FAC.py @@ -0,0 +1,35 @@ +''' + +read_FAC.py is a subfunction that reads atomic data created by +the Flexibly 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 + +''' + +# Module +from pfac import rfac +import os +import numpy as np + + +######################################### +# +# Main +# +######################################### + +def read_FAC( + # File Management + fil = None, + # Physics controls + EEDF = None, # if None -> assumes Maxwell-averages rates from pfac.fac.MaxwellRate + physics = None, # if None -> looks for all + ): \ No newline at end of file From bcca9086b8e566d4a5dd5a608169f357de03cfe4 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 18 Dec 2023 14:30:06 -0500 Subject: [PATCH 02/46] Started formating for FAC file management --- colradpy/read_FAC.py | 39 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 36 insertions(+), 3 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 0291e11..99ca7e1 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -28,8 +28,41 @@ def read_FAC( # File Management - fil = None, + fil = None, # Common path to FAC files, excluding physics extentions # Physics controls EEDF = None, # if None -> assumes Maxwell-averages rates from pfac.fac.MaxwellRate - physics = None, # if None -> looks for all - ): \ No newline at end of file + physics = None, # if None -> looks for all file suffixes + ): + + #### ---- Determines which files to search for ---- #### + + # Electron energy distribution settings + if EEDF is None: + # Use Maxwell-averaged rates from pfac.fac.MaxwellRate + use_mr = True + else: + print('NON-MAXWELLIAN ELECTRONS NOT IMPLEMENTED YET') + sys.exit(1) + + # FAC data file suffixes to search for + if physics is None: + # If already Maxwell-averaged rate coefficients + if use_mr: + physics = [ + 'a.en', # ASCII-format Energy levels + 'a.tr', # ASCII-format Einstein coefficients + 'ce.mr', # Collisional excitation + 'rr.mr', # Radiative recombination + #'ai.mr' # Autoionization/dielectronic recombination + 'ci.mr', # Collision ionization + ] + # If using cross-section files + else: + physics = [ + 'a.en', # ASCII-format Energy levels + 'a.tr', # ASCII-format Einstein coefficients + 'a.ce', # ASCII-format Collisional excitation + 'a.rr', # ASCII-format Radiative recombination + #'a.ai' # ASCII-format Autoionization/dielectronic recombination + 'a.ci', # ASCII-format Collision ionization + ] From daf1a912e8216a055ebc09dae779e765d4cb4269 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 18 Dec 2023 14:47:18 -0500 Subject: [PATCH 03/46] Init file loading management --- colradpy/read_FAC.py | 116 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 112 insertions(+), 4 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 99ca7e1..6cc2d12 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -20,11 +20,11 @@ import numpy as np -######################################### +############################################################ # -# Main +# Main # -######################################### +############################################################ def read_FAC( # File Management @@ -34,7 +34,7 @@ def read_FAC( physics = None, # if None -> looks for all file suffixes ): - #### ---- Determines which files to search for ---- #### + ######## -------- Determines which files to search for -------- ######## # Electron energy distribution settings if EEDF is None: @@ -66,3 +66,111 @@ def read_FAC( #'a.ai' # ASCII-format Autoionization/dielectronic recombination 'a.ci', # ASCII-format Collision ionization ] + + ######## -------- Reads data -------- ######## + + # Initialize output + FAC = {} + + # Energy levels + if 'a.en' in physics: + FAC = _en( + FAC=FAC, + fil=fil, + ) + # Error check + else: + print('NEED TO INCLUDE ENERGY LEVEL DATA IN MODELING!!!') + sys.exit(1) + + # Einstein coefficients + if 'a.tr' in physics: + FAC = _tr( + FAC=FAC, + fil=fil, + ) + # Error check + else: + print('NEED TO INCLUDE EINSTEIN COEFFICIENT DATA IN MODELING!!!') + sys.exit(1) + + # Collisional excitation + if 'a.ce' in physics: + FAC = _ce( + FAC=FAC, + fil=fil, + EEDF=EEDF, + ) + elif 'ce.mr' in physics: + FAC = _ce_mr( + FAC=FAC, + fil=fil + ) + # Error check + else: + print('NEED TO INCLUDE COLLISIONAL EXCITATION DATA IN MODELING!!!') + sys.exit(1) + + # Radiative recombination + if 'a.rr' in physics: + FAC = _rr( + FAC=FAC, + fil=fil, + EEDF=EEDF, + ) + elif 'rr.mr' in physics: + FAC = _rr_mr( + FAC=FAC, + fil=fil + ) + + # Autoionization/dielectronic recombination + if 'a.ai' in physics: + FAC = _ai( + FAC=FAC, + fil=fil, + EEDF=EEDF, + ) + elif 'ai.mr' in physics: + FAC = _ai_mr( + FAC=FAC, + fil=fil + ) + + # Collisional ionization + if 'a.ci' in physics: + FAC = _ai( + FAC=FAC, + fil=fil, + EEDF=EEDF, + ) + elif 'ci.mr' in physics: + FAC = _ci_mr( + FAC=FAC, + fil=fil + ) + + ######## -------- Formats Output -------- ######## + + # Output + return FAC + +############################################################ +# +# Utilities +# +############################################################ + +# Reads energy level file +def _en( + FAC=None, + fil=None, + ): + + # Initializes output + FAC['atomic'] = {} + + # Reads FAC energy levelfile + en = rfac.read_en( + fil+'a.en' + ) From 826e4e56e6e47e07e0db4f6af185764d8c219135 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 18 Dec 2023 15:14:51 -0500 Subject: [PATCH 04/46] Init function to read energy level FAC files --- colradpy/read_FAC.py | 90 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 89 insertions(+), 1 deletion(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 6cc2d12..358cb55 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -11,6 +11,7 @@ TO DO: 1) Add non-Maxwellian convolution + 2) Missing energy level data ''' @@ -27,6 +28,10 @@ ############################################################ 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 @@ -77,6 +82,9 @@ def read_FAC( FAC = _en( FAC=FAC, fil=fil, + ele=ele, + nele=nele, + Zele=Zele, ) # Error check else: @@ -157,7 +165,7 @@ def read_FAC( ############################################################ # -# Utilities +# File Reading # ############################################################ @@ -165,8 +173,14 @@ def read_FAC( def _en( FAC=None, fil=None, + ele=None, + nele=None, + Zele=None, ): + # Useful constants + ev2invcm = 8065.73 # [/1cm/eV] + # Initializes output FAC['atomic'] = {} @@ -174,3 +188,77 @@ def _en( en = rfac.read_en( fil+'a.en' ) + + # Charge states + FAC['atomic']['element'] = ele # Species name + FAC['atomic']['charge_state'] = nele # Electronic charge + FAC['atomic']['iz0'] = Zele # Nuclear charge + FAC['atomic']['iz1'] = 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 + + 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), !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + # Loop over blocks + for blk in en[1].keys(): + # 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][-1] + ) + + w.append( + en[1][blk]['2J'][st]/2 + ) + + energy.append( + en[1][blk]['ENERGY'][st] * eV2invcm + ) + + + # 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_pot_lvl.append( + en[1][blk]['ILEV'][st] + ) + + # 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']['ion_pot_lvl'] = np.asarray(ion_pot_lvl) + + # Output + return FAC + From 457100c13aa200564bd4bacf97f958191a2d7383 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 18 Dec 2023 15:37:02 -0500 Subject: [PATCH 05/46] Init function to read Einstein coefficient FAC files --- colradpy/read_FAC.py | 39 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 358cb55..bb8e8b0 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -76,6 +76,7 @@ def read_FAC( # Initialize output FAC = {} + FAC['rates'] = {} # Energy levels if 'a.en' in physics: @@ -244,7 +245,7 @@ def _en( ) ion_pot_lvl.append( - en[1][blk]['ILEV'][st] + en[1][blk]['ILEV'][st] +1 ) # Stores data @@ -262,3 +263,39 @@ def _en( # Output return FAC +# Reads Einstein coefficients +def _tr( + FAC=None, + fil=None, + ): + + # Init output dictionary + upr = [] # dim(ntrans,) + lwr = [] # dim(ntrans,) + a_val = [] # [1/s], dim(ntrans,) + + # 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] +1 + ) + + lwr.append( + tr[1][blk]['lower_index'][tran] +1 + ) + + a_val.append( + tr[1][blk]['rate'][tran] + ) + + # Formats output + FAC['rates']['a_val'] = np.asarray(a_val) # [1/s], dim(ntrans,) + trans = np.vstack((upr,lwr)).T # dim(ntrans,2) -> form for coll. excit transitions + + # Output + return FAC, trans \ No newline at end of file From 101ee79d0efe1d0cd61be9a593aa739c89a1d126 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 18 Dec 2023 15:56:52 -0500 Subject: [PATCH 06/46] Init function to read collisional excitation FAC files --- colradpy/read_FAC.py | 152 ++++++++++++++++++++++++++++++++++++++----- 1 file changed, 136 insertions(+), 16 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index bb8e8b0..7d5bfef 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -54,8 +54,8 @@ def read_FAC( # If already Maxwell-averaged rate coefficients if use_mr: physics = [ - 'a.en', # ASCII-format Energy levels - 'a.tr', # ASCII-format Einstein coefficients + 'en', # ASCII-format Energy levels + 'tr', # ASCII-format Einstein coefficients 'ce.mr', # Collisional excitation 'rr.mr', # Radiative recombination #'ai.mr' # Autoionization/dielectronic recombination @@ -64,12 +64,12 @@ def read_FAC( # If using cross-section files else: physics = [ - 'a.en', # ASCII-format Energy levels - 'a.tr', # ASCII-format Einstein coefficients - 'a.ce', # ASCII-format Collisional excitation - 'a.rr', # ASCII-format Radiative recombination - #'a.ai' # ASCII-format Autoionization/dielectronic recombination - 'a.ci', # ASCII-format Collision ionization + 'en', # ASCII-format Energy levels + 'tr', # ASCII-format Einstein coefficients + 'ce', # ASCII-format Collisional excitation + 'rr', # ASCII-format Radiative recombination + #'ai' # ASCII-format Autoionization/dielectronic recombination + 'ci', # ASCII-format Collision ionization ] ######## -------- Reads data -------- ######## @@ -79,7 +79,7 @@ def read_FAC( FAC['rates'] = {} # Energy levels - if 'a.en' in physics: + if 'en' in physics: FAC = _en( FAC=FAC, fil=fil, @@ -93,8 +93,8 @@ def read_FAC( sys.exit(1) # Einstein coefficients - if 'a.tr' in physics: - FAC = _tr( + if 'tr' in physics: + FAC, trans = _tr( FAC=FAC, fil=fil, ) @@ -104,16 +104,18 @@ def read_FAC( sys.exit(1) # Collisional excitation - if 'a.ce' in physics: + if 'ce' in physics: FAC = _ce( FAC=FAC, fil=fil, + trans=trans, EEDF=EEDF, ) elif 'ce.mr' in physics: FAC = _ce_mr( FAC=FAC, - fil=fil + fil=fil, + trans=trans, ) # Error check else: @@ -121,7 +123,7 @@ def read_FAC( sys.exit(1) # Radiative recombination - if 'a.rr' in physics: + if 'rr' in physics: FAC = _rr( FAC=FAC, fil=fil, @@ -134,7 +136,7 @@ def read_FAC( ) # Autoionization/dielectronic recombination - if 'a.ai' in physics: + if 'ai' in physics: FAC = _ai( FAC=FAC, fil=fil, @@ -298,4 +300,122 @@ def _tr( trans = np.vstack((upr,lwr)).T # dim(ntrans,2) -> form for coll. excit transitions # Output - return FAC, trans \ No newline at end of file + return FAC, trans + +# Reads Maxwell-averaged collisional excitation data files +def _ce_mr( + FAC=None, + fil=None, + trans=None, + ): + + # Reads data file + mr = _rad_mr( + fil=fil, + data='ce' + ) + + # Saves temperature grid data + FAC['temp_grid'] = np.asarray(mr['Te_eV']) # [eV], dim(nt,) + + # Initializes rate data + data = np.zeros((trans.shape[0], len(mr['Te_eV']))) # dim(ntrans,nt) + + # Loop over transitions + for tt in np.arange(trans.shape[0]): + # Upper and lower level indices + upr = int(trans[tt,0] -1) + lwr = int(trans[tt,1] - 1) + + # Saves transition rate coefficients, [cm3/s] + data[tt,:] = np.asarray( + data[lwr]['coll_excit'][upr] + ) + + # Formats output + FAC['rates']['excit'] = {} + FAC['rates']['excit']['col_transitions'] = trans # dim(ntrans,2), (upr,lwr) states + FAC['rates']['excit']['col_excit'] = data # dim(ntrans,nt), [cm3/s] + + # Output + return FAC + + + +############################################################ +# +# Utilities +# +############################################################ + + +# Read Maxwellian-averaged data files +def _read_mr( + fil = None, + data = None, + ): + + # Reads data file + f = open( + fil+data+'.mr', + 'r' + ) + + # Initializes output dictionary + out = {} + + # Data orginaization labels + if data == 'ce': + label = 'coll_excit' + elif data == 'rr': + label = 'rad_recomb' + elif data == '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(' ')[0]) + ) + + # Adds rate coefficient data, [cm3/s] + out[lwr][label][upr].append( + float(line.split(' ')[-1])*1e-10 + ) + + # Increases index to next temp point + indt += 1 + + # Output + return out From 6f608317b753ee921eaa133cdb6730bf48dfb14a Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 18 Dec 2023 16:09:47 -0500 Subject: [PATCH 07/46] Init function to read radiative recombination FAC files --- colradpy/read_FAC.py | 53 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 51 insertions(+), 2 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 7d5bfef..d1fd780 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -310,7 +310,7 @@ def _ce_mr( ): # Reads data file - mr = _rad_mr( + mr = _read_mr( fil=fil, data='ce' ) @@ -325,7 +325,7 @@ def _ce_mr( for tt in np.arange(trans.shape[0]): # Upper and lower level indices upr = int(trans[tt,0] -1) - lwr = int(trans[tt,1] - 1) + lwr = int(trans[tt,1] -1) # Saves transition rate coefficients, [cm3/s] data[tt,:] = np.asarray( @@ -340,6 +340,55 @@ def _ce_mr( # Output return FAC +# Reads Maxwell-averaged Radiative recombination data +def _rr_mr( + FAC=None, + fil=None, + ): + + # Reads data file + mr = _read_mr( + fil=fil, + data='rr' + ) + + # Init output + state = [] + ion = [] + data = [] + + # Loop over states with charge Z + for st in mr.keys(): + # Skips temp grid + if st == 'Te_eV': + continue + + state.append( + st+1 + ) + + # Loop over states with charge Z+1 + for ii in mr[st]['rad_recomb'].keys(): + ion.append( + ii+1 + ) + + data.append( + mr[st]['rad_recomb'][ii] + ) # [cm3/s] + + # Formats output + FAC['rates']['recomb'] = {} + FAC['rates']['recomb']['recomb_transitions'] = np.vstack( + (np.asarray(ion), np.asarray(st)) + ).T # dim(ntrans,2), Z+1 state -> Z state + FAC['rates']['recomb']['recomb_excit'] = np.asarray( + data + ) # dim(ntrans, nt), [cm3/s] + + # Ouput + return FAC + ############################################################ From 361b68ea43bd00ee6c17c55de45bcc8a03b20bfc Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 18 Dec 2023 16:21:39 -0500 Subject: [PATCH 08/46] Init function to read collisional ionization FAC files --- colradpy/read_FAC.py | 68 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 67 insertions(+), 1 deletion(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index d1fd780..9964c08 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -134,6 +134,11 @@ def read_FAC( FAC=FAC, fil=fil ) + # If empty + else: + FAC['rates']['recomb'] = {} + FAC['rates']['recomb']['recomb_transitions'] = np.asarray([]) + FAC['rates']['recomb']['recomb_excit'] = np.asarray([]) # Autoionization/dielectronic recombination if 'ai' in physics: @@ -149,7 +154,7 @@ def read_FAC( ) # Collisional ionization - if 'a.ci' in physics: + if 'ci' in physics: FAC = _ai( FAC=FAC, fil=fil, @@ -160,6 +165,20 @@ def read_FAC( FAC=FAC, fil=fil ) + # If empty + else: + FAC['rates']['ioniz'] = {} + FAC['rates']['ioniz']['ion_transitions'] = np.asarray([]) + FAC['rates']['ioniz']['ion_excit'] = np.asarray([]) + + # Charge exchange + if 'cx' in physics: + 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([]) ######## -------- Formats Output -------- ######## @@ -389,7 +408,54 @@ def _rr_mr( # Ouput return FAC +# Reads Maxwell-averaged collisional ionization data +def _ci_mr( + FAC=None, + fil=None, + ): + + # Reads data file + mr = _read_mr( + fil=fil, + data='ci' + ) + + # Init output + state = [] + ion = [] + data = [] + + # Loop over states with charge Z + for st in mr.keys(): + # Skips temp grid + if st == 'Te_eV': + continue + + state.append( + st+1 + ) + + # Loop over states with charge Z+1 + for ii in mr[st]['coll_ion'].keys(): + ion.append( + ii+1 + ) + + data.append( + mr[st]['coll_ion'][ii] + ) # [cm3/s] + + # Formats output + FAC['rates']['ioniz'] = {} + FAC['rates']['ioniz']['ion_transitions'] = np.vstack( + (np.asarray(st), np.asarray(ion)) + ).T # dim(ntrans,2), Z state -> Z+1 state + FAC['rates']['ioniz']['ion_excit'] = np.asarray( + data + ) # dim(ntrans, nt), [cm3/s] + # Ouput + return FAC ############################################################ # From 66a58008bc10806092f92b70c14d13e1380f9647 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 18 Dec 2023 16:28:42 -0500 Subject: [PATCH 09/46] Init function to format FAC output for ColRadPy --- colradpy/read_FAC.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 9964c08..dc762dd 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -1,7 +1,7 @@ ''' read_FAC.py is a subfunction that reads atomic data created by -the Flexibly Atomic Code (FAC) to be used in collisional- +the Flexible Atomic Code (FAC) to be used in collisional- radiative modeling in ColRadPy NOTE: Assumes you have a local version of FAC @@ -12,6 +12,7 @@ TO DO: 1) Add non-Maxwellian convolution 2) Missing energy level data + 3) Dielectronic recombination data ''' @@ -19,7 +20,7 @@ from pfac import rfac import os import numpy as np - +import copy ############################################################ # @@ -182,8 +183,13 @@ def read_FAC( ######## -------- Formats Output -------- ######## + out = {} + out['atomic'] = copy.deepcop(FAC['atomic']) + out['input_file'] = copy.deepcopy(FAC) + out['rates'] = copy.deepcopy(FAC['rates']) + # Output - return FAC + return out ############################################################ # @@ -463,7 +469,6 @@ def _ci_mr( # ############################################################ - # Read Maxwellian-averaged data files def _read_mr( fil = None, From 283cad6b666553eae66689d029bea1cbd1d3b8d8 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 18 Dec 2023 17:06:09 -0500 Subject: [PATCH 10/46] Init using FAC reading functions in populating CRM --- colradpy/colradpy_class.py | 31 ++++++++++++++++++++++++++++--- colradpy/read_FAC.py | 8 ++++---- 2 files changed, 32 insertions(+), 7 deletions(-) diff --git a/colradpy/colradpy_class.py b/colradpy/colradpy_class.py index ff47afe..a78d11e 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,21 @@ 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.updata_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'], + physics=self.data['user']['atomic_physics'], + ) + ) self.data['user']['file_loc'] = str(fil)#make this a string changed for hdfdict else: if(self.data): diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index dc762dd..1ee4d60 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -43,7 +43,7 @@ def read_FAC( ######## -------- Determines which files to search for -------- ######## # Electron energy distribution settings - if EEDF is None: + if EEDF == 'Maxwellian': # Use Maxwell-averaged rates from pfac.fac.MaxwellRate use_mr = True else: @@ -51,12 +51,12 @@ def read_FAC( sys.exit(1) # FAC data file suffixes to search for - if physics is None: + if physics == 'incl_all': # If already Maxwell-averaged rate coefficients if use_mr: physics = [ - 'en', # ASCII-format Energy levels - 'tr', # ASCII-format Einstein coefficients + 'en', # ASCII-format Energy levels + 'tr', # ASCII-format Einstein coefficients 'ce.mr', # Collisional excitation 'rr.mr', # Radiative recombination #'ai.mr' # Autoionization/dielectronic recombination From 044fa633e8b32ef4bd22e4e020c9b595ab521101 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 8 Jan 2024 15:50:27 -0500 Subject: [PATCH 11/46] [#18] debug reading FAC data --- colradpy/colradpy_class.py | 5 ++-- colradpy/read_FAC.py | 54 +++++++++++++++++++++++--------------- 2 files changed, 36 insertions(+), 23 deletions(-) diff --git a/colradpy/colradpy_class.py b/colradpy/colradpy_class.py index a78d11e..c38d278 100644 --- a/colradpy/colradpy_class.py +++ b/colradpy/colradpy_class.py @@ -342,7 +342,7 @@ def populate_data(self,fil): 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.updata_dict( + self.data = self.update_dict( self.data, read_FAC( ele=self.data['user']['ele'], nele=self.data['user']['nele'], @@ -2399,7 +2399,8 @@ 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'], diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 1ee4d60..c17aa15 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -11,8 +11,10 @@ TO DO: 1) Add non-Maxwellian convolution - 2) Missing energy level data - 3) Dielectronic recombination data + 2) Missing energy level data (zpla, zpla1) + 3) Dielectronic recombination/autoionization data + 4) Charge exchange data + 5) 2-photon emission data ''' @@ -184,7 +186,7 @@ def read_FAC( ######## -------- Formats Output -------- ######## out = {} - out['atomic'] = copy.deepcop(FAC['atomic']) + out['atomic'] = copy.deepcopy(FAC['atomic']) out['input_file'] = copy.deepcopy(FAC) out['rates'] = copy.deepcopy(FAC['rates']) @@ -207,7 +209,7 @@ def _en( ): # Useful constants - ev2invcm = 8065.73 # [/1cm/eV] + eV2invcm = 8065.73 # [/1cm/eV] # Initializes output FAC['atomic'] = {} @@ -218,10 +220,10 @@ def _en( ) # Charge states - FAC['atomic']['element'] = ele # Species name - FAC['atomic']['charge_state'] = nele # Electronic charge - FAC['atomic']['iz0'] = Zele # Nuclear charge - FAC['atomic']['iz1'] = nele +1 # Spectroscopic charge + 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 @@ -230,14 +232,14 @@ def _en( config = [] # dim(ntran,), state configuration L = [] # dim(ntran,), state L quantum number - S = [] # dim(ntran,), state S 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), !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + #zpla = [] # dim(ntran,nion), !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + #zpla1 = [] # dim(ntran,nion), !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # Loop over blocks - for blk in en[1].keys(): + for blk in np.arange(len(en[1])): # If block of excited states if en[1][blk]['NELE'] == nele: # Loop over states @@ -247,13 +249,17 @@ def _en( ) L.append( - en[1][blk]['VNL'][st][-1] + 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 ) @@ -283,8 +289,10 @@ def _en( 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'] = 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) # Output @@ -353,9 +361,13 @@ def _ce_mr( lwr = int(trans[tt,1] -1) # Saves transition rate coefficients, [cm3/s] - data[tt,:] = np.asarray( - data[lwr]['coll_excit'][upr] - ) + try: + data[tt,:] = np.asarray( + mr[lwr]['coll_excit'][upr] + ) + # Not all possible transitions included in file + except: + blah = 0 # Formats output FAC['rates']['excit'] = {} @@ -405,7 +417,7 @@ def _rr_mr( # Formats output FAC['rates']['recomb'] = {} FAC['rates']['recomb']['recomb_transitions'] = np.vstack( - (np.asarray(ion), np.asarray(st)) + (np.asarray(ion), np.asarray(state)) ).T # dim(ntrans,2), Z+1 state -> Z state FAC['rates']['recomb']['recomb_excit'] = np.asarray( data @@ -454,7 +466,7 @@ def _ci_mr( # Formats output FAC['rates']['ioniz'] = {} FAC['rates']['ioniz']['ion_transitions'] = np.vstack( - (np.asarray(st), np.asarray(ion)) + (np.asarray(state), np.asarray(ion)) ).T # dim(ntrans,2), Z state -> Z+1 state FAC['rates']['ioniz']['ion_excit'] = np.asarray( data @@ -531,7 +543,7 @@ def _read_mr( # Adds rate coefficient data, [cm3/s] out[lwr][label][upr].append( - float(line.split(' ')[-1])*1e-10 + float(line.replace(' ', ' ').split(' ')[-1])*1e-10 ) # Increases index to next temp point From 090aa5505b7c723372d6969e3edf1cfeae9b4d3c Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Tue, 9 Jan 2024 10:42:52 -0500 Subject: [PATCH 12/46] [#18] debugged level indexing --- colradpy/read_FAC.py | 91 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 70 insertions(+), 21 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index c17aa15..844f4c3 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -97,7 +97,7 @@ def read_FAC( # Einstein coefficients if 'tr' in physics: - FAC, trans = _tr( + FAC, trans_FAC = _tr( FAC=FAC, fil=fil, ) @@ -111,14 +111,14 @@ def read_FAC( FAC = _ce( FAC=FAC, fil=fil, - trans=trans, + trans_FAC=trans_FAC, EEDF=EEDF, ) elif 'ce.mr' in physics: FAC = _ce_mr( FAC=FAC, fil=fil, - trans=trans, + trans_FAC=trans_FAC, ) # Error check else: @@ -135,7 +135,7 @@ def read_FAC( elif 'rr.mr' in physics: FAC = _rr_mr( FAC=FAC, - fil=fil + fil=fil, ) # If empty else: @@ -153,7 +153,7 @@ def read_FAC( elif 'ai.mr' in physics: FAC = _ai_mr( FAC=FAC, - fil=fil + fil=fil, ) # Collisional ionization @@ -166,7 +166,7 @@ def read_FAC( elif 'ci.mr' in physics: FAC = _ci_mr( FAC=FAC, - fil=fil + fil=fil, ) # If empty else: @@ -183,6 +183,10 @@ def read_FAC( 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 = {} @@ -238,6 +242,13 @@ def _en( #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 @@ -264,6 +275,17 @@ def _en( 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: @@ -277,10 +299,22 @@ def _en( en[1][blk]['ENERGY'][st] * eV2invcm ) + # 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( - en[1][blk]['ILEV'][st] +1 + ind_i ) + # Stores data FAC['atomic']['ion_pot'] = np.asarray(ion_pot) FAC['atomic']['ion_term'] = np.asarray(ion_term) @@ -317,11 +351,11 @@ def _tr( # Loop over transitions for tran in np.arange(len(tr[1][blk]['lower_index'])): upr.append( - tr[1][blk]['upper_index'][tran] +1 + tr[1][blk]['upper_index'][tran] ) lwr.append( - tr[1][blk]['lower_index'][tran] +1 + tr[1][blk]['lower_index'][tran] ) a_val.append( @@ -330,16 +364,16 @@ def _tr( # Formats output FAC['rates']['a_val'] = np.asarray(a_val) # [1/s], dim(ntrans,) - trans = np.vstack((upr,lwr)).T # dim(ntrans,2) -> form for coll. excit transitions + trans_FAC = np.vstack((upr,lwr)).T # dim(ntrans,2) -> form for coll. excit transitions in FAC indices # Output - return FAC, trans + return FAC, trans_FAC # Reads Maxwell-averaged collisional excitation data files def _ce_mr( FAC=None, fil=None, - trans=None, + trans_FAC=None, ): # Reads data file @@ -352,13 +386,20 @@ def _ce_mr( FAC['temp_grid'] = np.asarray(mr['Te_eV']) # [eV], dim(nt,) # Initializes rate data - data = np.zeros((trans.shape[0], len(mr['Te_eV']))) # dim(ntrans,nt) + data = np.zeros((trans_FAC.shape[0], len(mr['Te_eV']))) # dim(ntrans,nt) + trans_ColRadPy = np.zeros(trans_FAC.shape, dtype='int') # Loop over transitions - for tt in np.arange(trans.shape[0]): + for tt in np.arange(trans_FAC.shape[0]): # Upper and lower level indices - upr = int(trans[tt,0] -1) - lwr = int(trans[tt,1] -1) + upr = int(trans_FAC[tt,0]) + lwr = int(trans_FAC[tt,1]) + + # Converts indices + ind_upr = np.where(FAC['lvl_indices']['FAC'] == upr) + ind_lwr = np.where(FAC['lvl_indices']['FAC'] == lwr) + trans_ColRadPy[tt,0] = FAC['lvl_indices']['ColRadPy'][ind_upr] + trans_ColRadPy[tt,1] = FAC['lvl_indices']['ColRadPy'][ind_lwr] # Saves transition rate coefficients, [cm3/s] try: @@ -371,7 +412,7 @@ def _ce_mr( # Formats output FAC['rates']['excit'] = {} - FAC['rates']['excit']['col_transitions'] = trans # dim(ntrans,2), (upr,lwr) states + FAC['rates']['excit']['col_transitions'] = trans_ColRadPy # dim(ntrans,2), (upr,lwr) states in ColRadPy indices FAC['rates']['excit']['col_excit'] = data # dim(ntrans,nt), [cm3/s] # Output @@ -400,14 +441,18 @@ def _rr_mr( if st == 'Te_eV': continue + # Converts indices + ind_st = np.where(FAC['lvl_indices']['FAC'] == st)[0][0] state.append( - st+1 + FAC['lvl_indices']['ColRadPy'][ind_st] ) # Loop over states with charge Z+1 for ii in mr[st]['rad_recomb'].keys(): + # Converts indices + ind_ion = np.where(FAC['lvl_indices']['FAC'] == ii)[0][0] ion.append( - ii+1 + FAC['lvl_indices']['ColRadPy'][ind_ion] ) data.append( @@ -449,14 +494,18 @@ def _ci_mr( if st == 'Te_eV': continue + # Converts indices + ind_st = np.where(FAC['lvl_indices']['FAC'] == st)[0][0] state.append( - st+1 + FAC['lvl_indices']['ColRadPy'][ind_st] ) # Loop over states with charge Z+1 for ii in mr[st]['coll_ion'].keys(): + # Converts indices + ind_ion = np.where(FAC['lvl_indices']['FAC'] == ii)[0][0] ion.append( - ii+1 + FAC['lvl_indices']['ColRadPy'][ind_ion] ) data.append( From 4503dc9fb1a173618ebab3746494457b0c776730 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Tue, 9 Jan 2024 15:01:40 -0500 Subject: [PATCH 13/46] [#18] debugged expected input file Te units --- colradpy/read_FAC.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 844f4c3..762f246 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -213,7 +213,7 @@ def _en( ): # Useful constants - eV2invcm = 8065.73 # [/1cm/eV] + eV2invcm = 8065.73 # [1/cm/eV] # Initializes output FAC['atomic'] = {} @@ -383,7 +383,9 @@ def _ce_mr( ) # Saves temperature grid data - FAC['temp_grid'] = np.asarray(mr['Te_eV']) # [eV], dim(nt,) + # NOTE: ColRadPy assumes Te is in Kelvin within the data file (per ADF04 standard) + eV2K = 11604.5 + FAC['temp_grid'] = np.asarray(mr['Te_eV'])*eV2K # [K], dim(nt,) # Initializes rate data data = np.zeros((trans_FAC.shape[0], len(mr['Te_eV']))) # dim(ntrans,nt) From c1f0b90484e4c5f1b082e526f1316983db9fa26e Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 18 Jan 2024 14:09:08 -0500 Subject: [PATCH 14/46] [#18] debugged having to convert excit rates into adf04 upsilon form --- colradpy/read_FAC.py | 40 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 38 insertions(+), 2 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 762f246..ce759fb 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -376,6 +376,9 @@ def _ce_mr( trans_FAC=None, ): + # Useful constants + eV2K = 11604.5 + # Reads data file mr = _read_mr( fil=fil, @@ -384,7 +387,6 @@ def _ce_mr( # Saves temperature grid data # NOTE: ColRadPy assumes Te is in Kelvin within the data file (per ADF04 standard) - eV2K = 11604.5 FAC['temp_grid'] = np.asarray(mr['Te_eV'])*eV2K # [K], dim(nt,) # Initializes rate data @@ -412,6 +414,15 @@ def _ce_mr( except: blah = 0 + # Converts data to expected reduced form (Upsilons per ADF04 standard) + data[tt,:] = _conv_rate2upsilon( + data = data[tt,:], + Te_eV = mr['Te_eV'], + ind_upr = trans_ColRadPy[tt,0] -1, + ind_lwr = trans_ColRadPy[tt,1] -1, + FAC = FAC, + ) + # Formats output FAC['rates']['excit'] = {} FAC['rates']['excit']['col_transitions'] = trans_ColRadPy # dim(ntrans,2), (upr,lwr) states in ColRadPy indices @@ -532,6 +543,31 @@ def _ci_mr( # ############################################################ +# Converts rate coefficient to adf04 Upsilon form +def _conv_rate2upsilon( + data = None, # [cm3/s], dim(ntemp,) + Te_eV = None, # [eV], dim(ntemp,) + ind_upr = None, + ind_lwr = None, + FAC = None, + ): + + # Useful constants + eV2invcm = 8065.73 # [1/cm/eV] + + return data * ( + np.sqrt(np.asarray(Te_eV)/13.6058) + /2.1716e-8 + *(1 + 2*FAC['atomic']['w'][ind_upr]) + * np.exp( + abs( + FAC['atomic']['energy'][ind_upr] + -FAC['atomic']['energy'][ind_lwr] + ) + /(np.asarray(Te_eV)*eV2invcm) + ) + ) + # Read Maxwellian-averaged data files def _read_mr( fil = None, @@ -589,7 +625,7 @@ def _read_mr( # If need to add temperature mesh if len(out['Te_eV']) < ntemp: out['Te_eV'].append( - float(line.split(' ')[0]) + float(line.split(' ')[1]) ) # Adds rate coefficient data, [cm3/s] From 961aed88a201ea16e81bb30dbb1a9346010b764b Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 18 Jan 2024 14:29:34 -0500 Subject: [PATCH 15/46] [#18] init function to convolve cross-sections with EEDF --- colradpy/convolve_EEDF.py | 69 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 colradpy/convolve_EEDF.py diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py new file mode 100644 index 0000000..93af188 --- /dev/null +++ b/colradpy/convolve_EEDF.py @@ -0,0 +1,69 @@ +''' + +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,) + engyXS = None, # [eV], dim (nE,) + m = 0, # 0 or 1 + ): + ''' + INPUTS: + EEDF -- [1/eV], + options: + 1) 'Maxwellian' --> pre-defined distribution function + if Te <10 keV = Maxwell-Boltzmann + if Te >10 keV = Maxwell-Juttner (relativistic) + + 2) 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,), cross-section data + + engyXS -- [eV], dim(nE,), cross-section energy axis + + m -- flag on definition of cross-section energy axis + options: + 1) if m=0, incident electron energy + 2) if m=1, scattered electron energy + + ''' + + # Prepares EEDF + if EEDF == 'Maxwellian': + EEDF = _get_Max( + Te=Te, + ) + else: + print('NON-MAXWELLIAN ELECTRONS NOT IMPLEMENTED YET') + sys.exit(1) + + +############################################################ +# +# Utilities +# +############################################################ \ No newline at end of file From 91fb794a168c0a3cc459383fa3872e0072729250 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 18 Jan 2024 15:14:47 -0500 Subject: [PATCH 16/46] [#18] function to calc Maxwellian dist --- colradpy/convolve_EEDF.py | 87 +++++++++++++++++++++++++++++++++++---- 1 file changed, 79 insertions(+), 8 deletions(-) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index 93af188..f6fb87b 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -23,15 +23,16 @@ def convolve_EEDF( EEDF = None, # [1/eV] Te = None, # [eV], dim(ntemp,) - XS = None, # [cm2], dim(nE,) - engyXS = None, # [eV], dim (nE,) + XS = None, # [cm2], dim(nE,ntrans) + engyXS = None, # [eV], dim (nE,ntrans) m = 0, # 0 or 1 + dE = 0, # [eV] ): ''' INPUTS: EEDF -- [1/eV], options: - 1) 'Maxwellian' --> pre-defined distribution function + 1) 'Maxwellian' --> pre-defined maximal entropy distribution function if Te <10 keV = Maxwell-Boltzmann if Te >10 keV = Maxwell-Juttner (relativistic) @@ -41,29 +42,99 @@ def convolve_EEDF( if using non-Maxwellian then whatever characteristic dimension that isn't the energy axis - XS -- [cm2], dim(nE,), cross-section data + XS -- [cm2], dim(nE,ntrans), cross-section data + nE -- number of energy grid points + ntrnas -- number of transitions considered - engyXS -- [eV], dim(nE,), cross-section energy axis + engyXS -- [eV], dim(nE,ntrans), cross-section energy axis 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 + ''' + # Useful constants + eV2electron = cnt.e /(cnt.m_e*cnt.c**2) # [1/eV] + # Prepares EEDF if EEDF == 'Maxwellian': - EEDF = _get_Max( + EEDF, engyEEDF = _get_Max( Te=Te, - ) + ) # dim(ngrid, ntemp) else: print('NON-MAXWELLIAN ELECTRONS NOT IMPLEMENTED YET') sys.exit(1) + # Loop over temperatures + for tt in np.arange(len(Te)): + # Incident electron velocity, relativistic form + vel = cnt.c * np.sqrt( + 1 - + 1/(engyEEDF[:,tt] *eV2electron +1)**2 + ) # [m/s], dim(ngrid,) + ############################################################ # # Utilities # -############################################################ \ No newline at end of file +############################################################ + +# Calculate Maxwellian energy distribution function +def _get_Max( + Te=None, # [eV], dim(ntemp,) + # Energy grid settings + ngrid = int(1e4), # energy grid resolution + lwr = 1e-3, # energy grid limits wrt multiple of Te + upr = 5e1, # !!!!!! Make sure these make sense for lines of interest + ): + ''' + 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] + + # low-energy form + if Te[tt] <= 10e3: + EEDF[:,tt] = ( + 2*np.sqrt(Einc/np.pi) + * temps[tt]**(-3/2) + * np.exp(-Einc/temps[tt]) + ) # dim(ngrid, ntemp); [1/eV] + + # relativistic form + else: + # 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] + + # Output + return EEDF, engy + From a27895bf456b894e594912868ae18ee8d12550d6 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 18 Jan 2024 15:29:17 -0500 Subject: [PATCH 17/46] [#18] function to preform rate coefficient integration --- colradpy/convolve_EEDF.py | 41 +++++++++++++++++++++++++++++++++++---- 1 file changed, 37 insertions(+), 4 deletions(-) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index f6fb87b..8449cca 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -47,6 +47,7 @@ def convolve_EEDF( ntrnas -- number of transitions considered engyXS -- [eV], dim(nE,ntrans), cross-section energy axis + NOTE: Assumed to be not extremely fine m -- flag on definition of cross-section energy axis options: @@ -61,22 +62,54 @@ def convolve_EEDF( # Useful constants eV2electron = cnt.e /(cnt.m_e*cnt.c**2) # [1/eV] + # Check + if ndim(XS) == 1: + XS = XS[:,None] + engyXS = engyXS[:,None] + # Prepares EEDF if EEDF == 'Maxwellian': EEDF, engyEEDF = _get_Max( Te=Te, - ) # dim(ngrid, ntemp) + ) # dim(ngrid, ntemp), ([1/eV], [eV]) else: print('NON-MAXWELLIAN ELECTRONS NOT IMPLEMENTED YET') sys.exit(1) + # Init output + ratec = np.zeros((len(Te))) # [cm3/s], dim(ntemp,ntrans) + # Loop over temperatures for tt in np.arange(len(Te)): # Incident electron velocity, relativistic form - vel = cnt.c * np.sqrt( + vel = cnt.c *100 *np.sqrt( 1 - 1/(engyEEDF[:,tt] *eV2electron +1)**2 - ) # [m/s], dim(ngrid,) + ) # [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 = (0,0) + )(np.log10(engyEEDF[:,tt])) # dim(ngrid,), [cm2] + + # Preforms integration + ratec[tt,nn] = np.trapz( + XS_tmp *vel * EEDF[:,tt], + engyEEDF[:,tt] + ) + + # Output, [cm3/s], dim(ntemp,ntrans) + return ratec ############################################################ @@ -91,7 +124,7 @@ def _get_Max( # Energy grid settings ngrid = int(1e4), # energy grid resolution lwr = 1e-3, # energy grid limits wrt multiple of Te - upr = 5e1, # !!!!!! Make sure these make sense for lines of interest + upr = 5e1, # !!!!!! Make sure these make sense for cross-sections ): ''' NOTE: EEDF energy axis defined as incident electron energy From ec67969fb8fc53364e4d2352146cfbb28bf33654 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 18 Jan 2024 15:55:20 -0500 Subject: [PATCH 18/46] [#18] init function to convolve EEDF for Dielectronic capture strengths --- colradpy/convolve_EEDF.py | 103 +++++++++++++++++++++++++++++++++----- 1 file changed, 90 insertions(+), 13 deletions(-) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index 8449cca..cc436a5 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -27,6 +27,7 @@ def convolve_EEDF( engyXS = None, # [eV], dim (nE,ntrans) m = 0, # 0 or 1 dE = 0, # [eV] + DC_flag = False, ): ''' INPUTS: @@ -57,13 +58,13 @@ def convolve_EEDF( 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 - # Useful constants - eV2electron = cnt.e /(cnt.m_e*cnt.c**2) # [1/eV] + ''' # Check - if ndim(XS) == 1: + if ndim(XS) == 1 and not DC_flag: XS = XS[:,None] engyXS = engyXS[:,None] @@ -76,11 +77,94 @@ def convolve_EEDF( 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 + ) # [cm3/s], dim(ntemp,ntrans) + + # Performs numerical integration + else: + return _calc_ratec( + EEDF=EEDF, + engyEEDF=engyEEDF, + XS=XS, + engyXS=engyXS + ) # [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, # [cm2], dim(ntrans,) + engyXS = None, # [eV], dim(ntrans,) + m = None, + dE = None, # [eV], dim(ntrans,) + ): + + # Useful constants + eV2electron = cnt.e /(cnt.m_e*cnt.c**2) # [1/eV] + + # 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)): + if m == 0: + engy_tmp = engyXS[nn] + elif m == 1: + engy_tmp = engyXS[nn] + dE[nn] + + # Incident electron velocity, relativistic form + vel = cnt.c *100 *np.sqrt( + 1 - + 1/(engy_tmp *eV2electron +1)**2 + ) # [cm/s] + + # 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 = (0,0) + )(np.log10(engy_tmp)) # [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) + m = None, + dE = None, # [eV], dim(ntrans,) + ): + + # Useful constants + eV2electron = cnt.e /(cnt.m_e*cnt.c**2) # [1/eV] + # Init output - ratec = np.zeros((len(Te))) # [cm3/s], dim(ntemp,ntrans) + ratec = np.zeros((EEDF.shape[1], XS.shape[1])) # [cm3/s], dim(ntemp,ntrans) # Loop over temperatures - for tt in np.arange(len(Te)): + for tt in np.arange(EEDF.shape[1]): # Incident electron velocity, relativistic form vel = cnt.c *100 *np.sqrt( 1 - @@ -111,13 +195,6 @@ def convolve_EEDF( # Output, [cm3/s], dim(ntemp,ntrans) return ratec - -############################################################ -# -# Utilities -# -############################################################ - # Calculate Maxwellian energy distribution function def _get_Max( Te=None, # [eV], dim(ntemp,) From 6f2c6be98697176d12f46a0fe0c6d81aae8f734a Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 18 Jan 2024 16:36:34 -0500 Subject: [PATCH 19/46] [#18] added 2-photon decay to A-values for H-like, He-like systems --- colradpy/read_FAC.py | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index ce759fb..deb3ee6 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -19,7 +19,7 @@ ''' # Module -from pfac import rfac +from pfac import rfac, crm import os import numpy as np import copy @@ -100,6 +100,8 @@ def read_FAC( FAC, trans_FAC = _tr( FAC=FAC, fil=fil, + nele=nele, + Zele=Zele, ) # Error check else: @@ -336,6 +338,8 @@ def _en( def _tr( FAC=None, fil=None, + nele=None, + Zele=None, ): # Init output dictionary @@ -366,6 +370,23 @@ def _tr( FAC['rates']['a_val'] = np.asarray(a_val) # [1/s], dim(ntrans,) trans_FAC = np.vstack((upr,lwr)).T # dim(ntrans,2) -> form for coll. excit transitions in FAC indices + # 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] + + FAC['rates']['a_val'][ind] += a_2photon + # Output return FAC, trans_FAC From 799fa6f9b6d89e638d9599a10bb4f668ee4474e8 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Fri, 19 Jan 2024 12:09:04 -0500 Subject: [PATCH 20/46] [#18] debugged convolve EEDF --- colradpy/__init__.py | 3 ++- colradpy/convolve_EEDF.py | 33 +++++++++++++++++++++------------ 2 files changed, 23 insertions(+), 13 deletions(-) diff --git a/colradpy/__init__.py b/colradpy/__init__.py index 6ac38b6..1ebebba 100644 --- a/colradpy/__init__.py +++ b/colradpy/__init__.py @@ -7,4 +7,5 @@ from .burgess_tully_rates import * from .split_multiplet import * from .nist_read_txt import * -from .read_FAC import * \ No newline at end of file +from .read_FAC import * +from .convolve_EEDF import * \ No newline at end of file diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index cc436a5..adf0ea5 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -45,10 +45,14 @@ def convolve_EEDF( XS -- [cm2], dim(nE,ntrans), cross-section data nE -- number of energy grid points - ntrnas -- number of transitions considered + 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 + 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: @@ -60,11 +64,12 @@ def convolve_EEDF( 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 ndim(XS) == 1 and not DC_flag: + if np.ndim(XS) == 1 and not DC_flag: XS = XS[:,None] engyXS = engyXS[:,None] @@ -83,7 +88,9 @@ def convolve_EEDF( EEDF=EEDF, engyEEDF=engyEEDF, XS=XS, - engyXS=engyXS + engyXS=engyXS, + m=m, + dE=dE, ) # [cm3/s], dim(ntemp,ntrans) # Performs numerical integration @@ -92,7 +99,9 @@ def convolve_EEDF( EEDF=EEDF, engyEEDF=engyEEDF, XS=XS, - engyXS=engyXS + engyXS=engyXS, + m=m, + dE=dE, ) # [cm3/s], dim(ntemp,ntrans) @@ -106,7 +115,7 @@ def convolve_EEDF( def _calc_DC( EEDF = None, # [1/eV], dim(ngrid,ntemp) engyEEDF = None, # [eV], dim(ngrid, ntemp) - XS = None, # [cm2], dim(ntrans,) + XS = None, # [eV*cm2], dim(ntrans,) engyXS = None, # [eV], dim(ntrans,) m = None, dE = None, # [eV], dim(ntrans,) @@ -138,7 +147,7 @@ def _calc_DC( np.log10(engyEEDF[:,tt]), np.log10(EEDF[:,tt]), bounds_error = False, - fill_value = (0,0) + fill_value = (-1e5,-1e5) )(np.log10(engy_tmp)) # [1/eV] # Calculates rate coefficient, [cm3/s] @@ -183,7 +192,7 @@ def _calc_ratec( np.log10(engy_tmp), np.log10(XS[:,nn]), bounds_error=False, - fill_value = (0,0) + fill_value = (-1e5,-1e5) )(np.log10(engyEEDF[:,tt])) # dim(ngrid,), [cm2] # Preforms integration @@ -199,7 +208,7 @@ def _calc_ratec( def _get_Max( Te=None, # [eV], dim(ntemp,) # Energy grid settings - ngrid = int(1e4), # energy grid resolution + ngrid = int(5e2), # energy grid resolution lwr = 1e-3, # energy grid limits wrt multiple of Te upr = 5e1, # !!!!!! Make sure these make sense for cross-sections ): @@ -225,9 +234,9 @@ def _get_Max( # low-energy form if Te[tt] <= 10e3: EEDF[:,tt] = ( - 2*np.sqrt(Einc/np.pi) - * temps[tt]**(-3/2) - * np.exp(-Einc/temps[tt]) + 2*np.sqrt(engy[:,tt]/np.pi) + * Te[tt]**(-3/2) + * np.exp(-engy[:,tt]/Te[tt]) ) # dim(ngrid, ntemp); [1/eV] # relativistic form From bfcdfbc1341552bde00f71d98854c0d611f11637 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Fri, 19 Jan 2024 13:26:55 -0500 Subject: [PATCH 21/46] [#18] corrected form of ionization rates to adf04 reduced form --- colradpy/read_FAC.py | 43 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 40 insertions(+), 3 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index deb3ee6..447fa5a 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -447,7 +447,7 @@ def _ce_mr( # Formats output FAC['rates']['excit'] = {} FAC['rates']['excit']['col_transitions'] = trans_ColRadPy # dim(ntrans,2), (upr,lwr) states in ColRadPy indices - FAC['rates']['excit']['col_excit'] = data # dim(ntrans,nt), [cm3/s] + FAC['rates']['excit']['col_excit'] = data # dim(ntrans,nt), [upsilon] # Output return FAC @@ -546,6 +546,15 @@ def _ci_mr( mr[st]['coll_ion'][ii] ) # [cm3/s] + # Converts ionization data to adf04 reduced form + data[-1] = _conv_rate2reduct( + data = data[-1], + Te_eV = mr['Te_eV'], + ind_st = state[-1] -1, + ind_ion = ion[-1] -1, + FAC = FAC, + ) + # Formats output FAC['rates']['ioniz'] = {} FAC['rates']['ioniz']['ion_transitions'] = np.vstack( @@ -553,7 +562,7 @@ def _ci_mr( ).T # dim(ntrans,2), Z state -> Z+1 state FAC['rates']['ioniz']['ion_excit'] = np.asarray( data - ) # dim(ntrans, nt), [cm3/s] + ) # dim(ntrans, nt), [reduced rate] # Ouput return FAC @@ -564,9 +573,37 @@ def _ci_mr( # ############################################################ +# 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,) + data = None, # [cm3/s], dim(ntemp,), array Te_eV = None, # [eV], dim(ntemp,) ind_upr = None, ind_lwr = None, From d3d4a17bb3563edf161a581a4ae2b9b3c2db561d Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Wed, 24 Jan 2024 13:37:58 -0500 Subject: [PATCH 22/46] [#18] fixed indexing bug when loading ce data causing it to skip some transitions --- colradpy/read_FAC.py | 388 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 329 insertions(+), 59 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 447fa5a..0c4be48 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -23,6 +23,7 @@ import os import numpy as np import copy +from colradpy.convolve_EEDF import convolve_EEDF as cE ############################################################ # @@ -40,14 +41,17 @@ def read_FAC( # Physics controls EEDF = None, # if None -> assumes Maxwell-averages rates from pfac.fac.MaxwellRate physics = None, # if None -> looks for all file suffixes + Te = None, # if not using MaxwellRate files ): ######## -------- Determines which files to search for -------- ######## # Electron energy distribution settings - if EEDF == 'Maxwellian': + if EEDF == 'Maxwellian_mr': # Use Maxwell-averaged rates from pfac.fac.MaxwellRate use_mr = True + elif EEDF == 'Maxwellian': + use_mr = False else: print('NON-MAXWELLIAN ELECTRONS NOT IMPLEMENTED YET') sys.exit(1) @@ -61,7 +65,7 @@ def read_FAC( 'tr', # ASCII-format Einstein coefficients 'ce.mr', # Collisional excitation 'rr.mr', # Radiative recombination - #'ai.mr' # Autoionization/dielectronic recombination + #'ai', # Autoionization/dielectronic recombination 'ci.mr', # Collision ionization ] # If using cross-section files @@ -70,9 +74,9 @@ def read_FAC( 'en', # ASCII-format Energy levels 'tr', # ASCII-format Einstein coefficients 'ce', # ASCII-format Collisional excitation - 'rr', # ASCII-format Radiative recombination - #'ai' # ASCII-format Autoionization/dielectronic recombination - 'ci', # ASCII-format Collision ionization + #'rr', # ASCII-format Radiative recombination + #'ai', # ASCII-format Autoionization/dielectronic recombination + #'ci', # ASCII-format Collision ionization ] ######## -------- Reads data -------- ######## @@ -95,38 +99,37 @@ def read_FAC( print('NEED TO INCLUDE ENERGY LEVEL DATA IN MODELING!!!') sys.exit(1) - # Einstein coefficients - if 'tr' in physics: - FAC, trans_FAC = _tr( - FAC=FAC, - fil=fil, - nele=nele, - Zele=Zele, - ) - # Error check - else: - print('NEED TO INCLUDE EINSTEIN COEFFICIENT DATA IN MODELING!!!') - sys.exit(1) - # Collisional excitation if 'ce' in physics: FAC = _ce( FAC=FAC, fil=fil, - trans_FAC=trans_FAC, EEDF=EEDF, + Te=Te, ) elif 'ce.mr' in physics: FAC = _ce_mr( FAC=FAC, fil=fil, - trans_FAC=trans_FAC, ) # Error check else: print('NEED TO INCLUDE COLLISIONAL EXCITATION DATA IN MODELING!!!') sys.exit(1) + # Einstein coefficients + if 'tr' in physics: + FAC = _tr( + FAC=FAC, + fil=fil, + nele=nele, + Zele=Zele, + ) + # Error check + else: + print('NEED TO INCLUDE EINSTEIN COEFFICIENT DATA IN MODELING!!!') + sys.exit(1) + # Radiative recombination if 'rr' in physics: FAC = _rr( @@ -346,6 +349,8 @@ def _tr( 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') @@ -366,9 +371,18 @@ def _tr( tr[1][blk]['rate'][tran] ) - # Formats output - FAC['rates']['a_val'] = np.asarray(a_val) # [1/s], dim(ntrans,) - trans_FAC = np.vstack((upr,lwr)).T # dim(ntrans,2) -> form for coll. excit transitions in FAC indices + # 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) # Includes two-photon emission # Only works for H-like and He-like @@ -385,16 +399,203 @@ def _tr( & (FAC['atomic']['w'] == 0) )[0][0] - FAC['rates']['a_val'][ind] += a_2photon + ind2 = np.where( + np.all( + trans_ColRadPy == np.asarray([ + FAC['lvl_indices']['ColRadPy'][ind], 1 + ]), + axis = 1 + ) + )[0][0] + + a_val[ind2] += a_2photon + + ## 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(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],:] + + # 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, trans_FAC + return FAC + +# Reads collisional excitation cross-section data files +def _ce( + EEDF=None, + Te=None, # [eV], dim(ntemp,) + FAC=None, + fil=None, + trans_FAC=None, + vebose = 1, + ): + + # Useful constants + eV2K = 11604.5 + + # Saves temperature grid data + # NOTE: ColRadPy assumes Te is in Kelvin within the data file (per ADF04 standard) + FAC['temp_grid'] = Te*eV2K # [K], dim(nt,) + + # Initializes rate data + data = np.zeros((trans_FAC.shape[0], len(Te))) # dim(ntrans,nt) + trans_ColRadPy = np.zeros(trans_FAC.shape, dtype='int') + XS = None + + # Reads FAC data file + ce = rfac.read_ce(fil+'a.ce') + + # Loop over blocks + for blk in np.arange(len(ce[1])): + # Loop over transitions in block + for trn in np.arange(len(ce[1][blk]['Delta E'])): + # Init data storage + if XS is None: + XS = np.zeros((len(ce[1][blk]['EGRID']),trans_FAC.shape[0])) # dim(nE, ntrans) + engyXS = np.zeros((len(ce[1][blk]['EGRID']),trans_FAC.shape[0])) # dim(nE, ntrans) + dE = np.zeros(trans_FAC.shape[0]) + Bethe = np.zeros((trans_FAC.shape[0], 2)) + w_upr = np.zeros(trans_FAC.shape[0]) + + + # Upper and lower level indices + upr = ce[1][blk]['upper_index'][trn] + lwr = ce[1][blk]['lower_index'][trn] + + # Converts indices + ind_upr = np.where(FAC['lvl_indices']['FAC'] == upr) + ind_lwr = np.where(FAC['lvl_indices']['FAC'] == lwr) + + try: + ind_tt = np.where( + (trans_FAC[:,0] == upr) + & (trans_FAC[:,1] == lwr) + )[0][0] + except: + print(upr) + print(lwr) + print('xxxxx') + continue + + #import pdb + #pdb.set_trace() + + trans_ColRadPy[ind_tt,0] = FAC['lvl_indices']['ColRadPy'][ind_upr] + trans_ColRadPy[ind_tt,1] = FAC['lvl_indices']['ColRadPy'][ind_lwr] + + # Stores data + XS[:,ind_tt] = ce[1][blk]['crosssection'][trn,:]*1e-20 + engyXS[:,ind_tt] = ce[1][blk]['EGRID'] + dE[ind_tt] = ce[1][blk]['TE0'] + Bethe[ind_tt,0] = ce[1][blk]['bethe'][trn] + Bethe[ind_tt,1] = ce[1][blk]['born'][trn,0] + w_upr[ind_tt] = FAC['atomic']['w'][ind_upr] + + # Calculates Rate coefficient data, [cm3/s], dim(ntrans, ntemp) + if verbose == 1: + ( + data, + )= cE( + EEDF = EEDF, + Te=Te, + XS = XS, + engyXS = engyXS, + m = ce[1][blk]['ETYPE'], + dE = dE, + Bethe = Bethe, + #Bethe = None, + w_upr = w_upr, + verbose=verbose, + ) + data = data.T + + else: + data = cE( + EEDF = EEDF, + Te=Te, + XS = XS, + engyXS = engyXS, + m = ce[1][blk]['ETYPE'], + dE = dE, + Bethe = Bethe, + #Bethe = None, + w_upr = w_upr, + verbose=verbose, + ).T + + #import pdb + #pdb.set_trace() + if verbose == 1: + ratec = data.copy() + + + # Converts data to expected reduced form (Upsilons per ADF04 standard) + for tt in np.arange(trans_FAC.shape[0]): + data[tt,:] = _conv_rate2upsilon( + data = data[tt,:], + Te_eV = Te, + ind_upr = trans_ColRadPy[tt,0] -1, + ind_lwr = trans_ColRadPy[tt,1] -1, + FAC = FAC, + ) + + + # Formats output + FAC['rates']['excit'] = {} + FAC['rates']['excit']['col_transitions'] = trans_ColRadPy # dim(ntrans,2), (upr,lwr) states in ColRadPy indices + FAC['rates']['excit']['col_excit'] = data # dim(ntrans,nt), [upsilon] + if verbose == 1: + FAC['rates']['excit']['col_trans_unfill'] = FAC['rates']['excit']['col_transitions'].copy() + FAC['rates']['excit']['ratec_cm3/s'] = ratec # dim(ntrans,nt), [cm3/s] + FAC['rates']['excit']['Bethe'] = Bethe # dim(ntrans, 2) + + # Output + return FAC # Reads Maxwell-averaged collisional excitation data files def _ce_mr( FAC=None, fil=None, - trans_FAC=None, + verbose = 1, ): # Useful constants @@ -410,44 +611,55 @@ def _ce_mr( # 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,) - # Initializes rate data - data = np.zeros((trans_FAC.shape[0], len(mr['Te_eV']))) # dim(ntrans,nt) - trans_ColRadPy = np.zeros(trans_FAC.shape, dtype='int') + # Init output + st_lwr = [] + st_upr = [] + ratec = [] - # Loop over transitions - for tt in np.arange(trans_FAC.shape[0]): - # Upper and lower level indices - upr = int(trans_FAC[tt,0]) - lwr = int(trans_FAC[tt,1]) + # Loop over lower states + for lwr in mr.keys(): + # Skips temp grid + if lwr == 'Te_eV': + continue # Converts indices - ind_upr = np.where(FAC['lvl_indices']['FAC'] == upr) - ind_lwr = np.where(FAC['lvl_indices']['FAC'] == lwr) - trans_ColRadPy[tt,0] = FAC['lvl_indices']['ColRadPy'][ind_upr] - trans_ColRadPy[tt,1] = FAC['lvl_indices']['ColRadPy'][ind_lwr] - - # Saves transition rate coefficients, [cm3/s] - try: - data[tt,:] = np.asarray( - mr[lwr]['coll_excit'][upr] + ind_lwr = np.where(FAC['lvl_indices']['FAC'] == lwr)[0][0] + + # Loop over upper states + for upr in mr[lwr]['coll_excit'].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] ) - # Not all possible transitions included in file - except: - blah = 0 - # Converts data to expected reduced form (Upsilons per ADF04 standard) + ratec.append( + mr[lwr]['coll_excit'][upr] + ) # [cm3/s] + + # Convert to Upsilon form + data = np.zeros((len(ratec), len(mr['Te_eV']))) # dim(ntrans, nt) + for tt in np.arange(len(ratec)): data[tt,:] = _conv_rate2upsilon( - data = data[tt,:], + data = ratec[tt], Te_eV = mr['Te_eV'], - ind_upr = trans_ColRadPy[tt,0] -1, - ind_lwr = trans_ColRadPy[tt,1] -1, + ind_upr = st_upr[tt] -1, + ind_lwr = st_lwr[tt] -1, FAC = FAC, ) - + # Formats output FAC['rates']['excit'] = {} - FAC['rates']['excit']['col_transitions'] = trans_ColRadPy # dim(ntrans,2), (upr,lwr) states in ColRadPy indices - FAC['rates']['excit']['col_excit'] = data # dim(ntrans,nt), [upsilon] + FAC['rates']['excit']['col_transitions'] = np.vstack( + (np.asarray(st_upr), np.asarray(st_lwr)) + ).T # dim(ntrans,2), (upr,lwr) states in ColRadPy indices + FAC['rates']['excit']['col_excit'] = data # dim(ntrans, nt), [cm3/s] + if verbose == 1: + FAC['rates']['excit']['col_trans_unfill'] = FAC['rates']['excit']['col_transitions'].copy() + FAC['rates']['excit']['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, nt), [cm3/s] # Output return FAC @@ -477,14 +689,14 @@ def _rr_mr( # Converts indices ind_st = np.where(FAC['lvl_indices']['FAC'] == st)[0][0] - state.append( - FAC['lvl_indices']['ColRadPy'][ind_st] - ) # Loop over states with charge Z+1 for ii in mr[st]['rad_recomb'].keys(): # Converts indices ind_ion = np.where(FAC['lvl_indices']['FAC'] == ii)[0][0] + state.append( + FAC['lvl_indices']['ColRadPy'][ind_st] + ) ion.append( FAC['lvl_indices']['ColRadPy'][ind_ion] ) @@ -530,14 +742,14 @@ def _ci_mr( # Converts indices ind_st = np.where(FAC['lvl_indices']['FAC'] == st)[0][0] - state.append( - FAC['lvl_indices']['ColRadPy'][ind_st] - ) # Loop over states with charge Z+1 for ii in mr[st]['coll_ion'].keys(): # Converts indices ind_ion = np.where(FAC['lvl_indices']['FAC'] == ii)[0][0] + state.append( + FAC['lvl_indices']['ColRadPy'][ind_st] + ) ion.append( FAC['lvl_indices']['ColRadPy'][ind_ion] ) @@ -626,6 +838,64 @@ def _conv_rate2upsilon( ) ) +# Reads cross-section data files +def _read_ascii( + fil = None, + data = None, + ): + + # Reads data file + if data == 'ce': + data_fil = rfac.read_ce(fil+'a.ce') + lwr_lbl = 'lower_index' + upr_lbl = 'upper_index' + data_lbl = 'crosssection' + elif data == 'rr': + data_fil = rfac.read_rr(fil+'a.ce') + lwr_lbl = 'bound_index' + upr_lbl = 'free_index' + data_lbl = 'RR crosssection' + elif data == 'ci': + data_fil = rfac.read_rr(fil+'a.ce') + 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 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 parameters for high-energy behavior + if data == 'ce': + out[lwr][upr]['limit'] = np.asarray([ + data_fil[1][blk]['bethe'][trn], + data_fil[1][blk]['born'][trn,0] + ]) + + # Output + return out + # Read Maxwellian-averaged data files def _read_mr( fil = None, From a77898d5b7cbd219c52cf986b5cb24683b41dbd8 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Wed, 24 Jan 2024 14:48:08 -0500 Subject: [PATCH 23/46] [#18] gneralized redundant code in handling mr file data --- colradpy/read_FAC.py | 395 ++++++++++++++++++++----------------------- 1 file changed, 185 insertions(+), 210 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 0c4be48..adcb878 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -23,7 +23,7 @@ import os import numpy as np import copy -from colradpy.convolve_EEDF import convolve_EEDF as cE +from colradpy.convolve_EEDF import convolve_EEDF ############################################################ # @@ -41,7 +41,8 @@ def read_FAC( # Physics controls EEDF = None, # if None -> assumes Maxwell-averages rates from pfac.fac.MaxwellRate physics = None, # if None -> looks for all file suffixes - Te = None, # if not using MaxwellRate files + Te = None, # if not using MaxwellRate files, [eV], dim(ntemp,) + verbose = 1, ): ######## -------- Determines which files to search for -------- ######## @@ -106,11 +107,13 @@ def read_FAC( fil=fil, EEDF=EEDF, Te=Te, + verbose=verbose, ) elif 'ce.mr' in physics: FAC = _ce_mr( FAC=FAC, fil=fil, + verbose=verbose, ) # Error check else: @@ -172,6 +175,7 @@ def read_FAC( FAC = _ci_mr( FAC=FAC, fil=fil, + verbose=verbose, ) # If empty else: @@ -466,7 +470,7 @@ def _ce( FAC=None, fil=None, trans_FAC=None, - vebose = 1, + vebose = None, ): # Useful constants @@ -476,126 +480,105 @@ def _ce( # NOTE: ColRadPy assumes Te is in Kelvin within the data file (per ADF04 standard) FAC['temp_grid'] = Te*eV2K # [K], dim(nt,) - # Initializes rate data - data = np.zeros((trans_FAC.shape[0], len(Te))) # dim(ntrans,nt) - trans_ColRadPy = np.zeros(trans_FAC.shape, dtype='int') - XS = None - - # Reads FAC data file - ce = rfac.read_ce(fil+'a.ce') - - # Loop over blocks - for blk in np.arange(len(ce[1])): - # Loop over transitions in block - for trn in np.arange(len(ce[1][blk]['Delta E'])): - # Init data storage - if XS is None: - XS = np.zeros((len(ce[1][blk]['EGRID']),trans_FAC.shape[0])) # dim(nE, ntrans) - engyXS = np.zeros((len(ce[1][blk]['EGRID']),trans_FAC.shape[0])) # dim(nE, ntrans) - dE = np.zeros(trans_FAC.shape[0]) - Bethe = np.zeros((trans_FAC.shape[0], 2)) - w_upr = np.zeros(trans_FAC.shape[0]) - - - # Upper and lower level indices - upr = ce[1][blk]['upper_index'][trn] - lwr = ce[1][blk]['lower_index'][trn] + # Loads date from ascii file + ce = _read_ascii( + fil = fil, + data = 'ce', + ) - # Converts indices - ind_upr = np.where(FAC['lvl_indices']['FAC'] == upr) - ind_lwr = np.where(FAC['lvl_indices']['FAC'] == lwr) - - try: - ind_tt = np.where( - (trans_FAC[:,0] == upr) - & (trans_FAC[:,1] == lwr) - )[0][0] - except: - print(upr) - print(lwr) - print('xxxxx') - continue - - #import pdb - #pdb.set_trace() - - trans_ColRadPy[ind_tt,0] = FAC['lvl_indices']['ColRadPy'][ind_upr] - trans_ColRadPy[ind_tt,1] = FAC['lvl_indices']['ColRadPy'][ind_lwr] - - # Stores data - XS[:,ind_tt] = ce[1][blk]['crosssection'][trn,:]*1e-20 - engyXS[:,ind_tt] = ce[1][blk]['EGRID'] - dE[ind_tt] = ce[1][blk]['TE0'] - Bethe[ind_tt,0] = ce[1][blk]['bethe'][trn] - Bethe[ind_tt,1] = ce[1][blk]['born'][trn,0] - w_upr[ind_tt] = FAC['atomic']['w'][ind_upr] - - # Calculates Rate coefficient data, [cm3/s], dim(ntrans, ntemp) + # Init output + st_lwr = [] + st_upr = [] + ratec = [] if verbose == 1: - ( - data, - )= cE( - EEDF = EEDF, - Te=Te, - XS = XS, - engyXS = engyXS, - m = ce[1][blk]['ETYPE'], - dE = dE, - Bethe = Bethe, - #Bethe = None, - w_upr = w_upr, - verbose=verbose, - ) - data = data.T - - else: - data = cE( - EEDF = EEDF, - Te=Te, - XS = XS, - engyXS = engyXS, - m = ce[1][blk]['ETYPE'], - dE = dE, - Bethe = Bethe, - #Bethe = None, - w_upr = w_upr, - verbose=verbose, - ).T + XS = [] + engy = [] - #import pdb - #pdb.set_trace() - if verbose == 1: - ratec = data.copy() + # Loop over lower states + for lwr in ce.keys(): + # Converts indices + ind_lwr = np.where(FAC['lvl_indices']['FAC'] == lwr)[0][0] + + # Loop over upper states + for upr in ce[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] + ) + # Calculates Rate coefficient data, [cm3/s], dim(ntrans, ntemp) + if verbose == 1: + ( + tmp_ratec, + tmp_XS, + tmp_engy + )= convolve_EEDF( + EEDF = EEDF, + Te=Te, + XS = ce[lwr][upr]['XS'], + engyXS = ce[lwr][upr]['engy'], + m = 0, + dE = np.asarray([ce[lwr][upr]['dE']]), + Bethe = ce[lwr][upr]['limit'][None,:], + w_upr = np.asarray([ce[lwr][upr]['w_upr']]), + verbose=verbose, + use_rel = True, + ) + ratec.append(tmp_ratec[:,0]) + XS.append(tmp_XS[:,0]) + engy.append(tmp_engy[:,0]) + + else: + ratec.append( + convolve_EEDF( + EEDF = EEDF, + Te=Te, + XS = ce[lwr][upr]['XS'], + engyXS = ce[lwr][upr]['engy'], + m = 0, + dE = np.asarray([ce[lwr][upr]['dE']]), + Bethe = ce[lwr][upr]['limit'][None,:], + w_upr = np.asarray([ce[lwr][upr]['w_upr']]), + verbose=verbose, + use_rel = True, + )[:,0] + ) - # Converts data to expected reduced form (Upsilons per ADF04 standard) - for tt in np.arange(trans_FAC.shape[0]): + # Convert to Upsilon form + data = np.zeros((len(ratec), len(Te))) # dim(ntrans, ntemp) + for tt in np.arange(len(ratec)): data[tt,:] = _conv_rate2upsilon( - data = data[tt,:], + data = ratec[tt], Te_eV = Te, - ind_upr = trans_ColRadPy[tt,0] -1, - ind_lwr = trans_ColRadPy[tt,1] -1, + ind_upr = st_upr[tt] -1, + ind_lwr = st_lwr[tt] -1, FAC = FAC, ) - # Formats output FAC['rates']['excit'] = {} - FAC['rates']['excit']['col_transitions'] = trans_ColRadPy # dim(ntrans,2), (upr,lwr) states in ColRadPy indices - FAC['rates']['excit']['col_excit'] = data # dim(ntrans,nt), [upsilon] + FAC['rates']['excit']['col_transitions'] = np.vstack( + (np.asarray(st_upr), np.asarray(st_lwr)) + ).T # dim(ntrans,2), (upr,lwr) states in ColRadPy indices + FAC['rates']['excit']['col_excit'] = data # dim(ntrans, ntemp), [upsilon] if verbose == 1: FAC['rates']['excit']['col_trans_unfill'] = FAC['rates']['excit']['col_transitions'].copy() - FAC['rates']['excit']['ratec_cm3/s'] = ratec # dim(ntrans,nt), [cm3/s] - FAC['rates']['excit']['Bethe'] = Bethe # dim(ntrans, 2) + FAC['rates']['excit']['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, ntemp), [cm3/s] + FAC['rates']['excit']['XS_cm2'] = np.asarray(XS) # dim(ntrans, nE), [cm2] + FAC['rates']['excit']['engy_eV'] = np.asarray(engy) # dim(ntrans, nE), [eV] - # Output + # Output return FAC # Reads Maxwell-averaged collisional excitation data files def _ce_mr( FAC=None, fil=None, - verbose = 1, + verbose = None, ): # Useful constants @@ -611,52 +594,20 @@ def _ce_mr( # 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,) - # Init output - st_lwr = [] - st_upr = [] - ratec = [] - - # Loop over lower states - 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 upper states - for upr in mr[lwr]['coll_excit'].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] - ) - - ratec.append( - mr[lwr]['coll_excit'][upr] - ) # [cm3/s] - - # Convert to Upsilon form - data = np.zeros((len(ratec), len(mr['Te_eV']))) # dim(ntrans, nt) - for tt in np.arange(len(ratec)): - data[tt,:] = _conv_rate2upsilon( - data = ratec[tt], - Te_eV = mr['Te_eV'], - ind_upr = st_upr[tt] -1, - ind_lwr = st_lwr[tt] -1, - FAC = FAC, - ) + # Converts data file to ColRadPy form + data, st_upr, st_lwr, ratec = _conv_mr2colradpy( + FAC = FAC, + mr = mr, + react = 'ce', + verbose=verbose, + ) # Formats output FAC['rates']['excit'] = {} FAC['rates']['excit']['col_transitions'] = np.vstack( (np.asarray(st_upr), np.asarray(st_lwr)) ).T # dim(ntrans,2), (upr,lwr) states in ColRadPy indices - FAC['rates']['excit']['col_excit'] = data # dim(ntrans, nt), [cm3/s] + FAC['rates']['excit']['col_excit'] = np.asarray(data) # dim(ntrans, nt), [cm3/s] if verbose == 1: FAC['rates']['excit']['col_trans_unfill'] = FAC['rates']['excit']['col_transitions'].copy() FAC['rates']['excit']['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, nt), [cm3/s] @@ -676,43 +627,19 @@ def _rr_mr( data='rr' ) - # Init output - state = [] - ion = [] - data = [] - - # Loop over states with charge Z - for st in mr.keys(): - # Skips temp grid - if st == 'Te_eV': - continue - - # Converts indices - ind_st = np.where(FAC['lvl_indices']['FAC'] == st)[0][0] - - # Loop over states with charge Z+1 - for ii in mr[st]['rad_recomb'].keys(): - # Converts indices - ind_ion = np.where(FAC['lvl_indices']['FAC'] == ii)[0][0] - state.append( - FAC['lvl_indices']['ColRadPy'][ind_st] - ) - ion.append( - FAC['lvl_indices']['ColRadPy'][ind_ion] - ) - - data.append( - mr[st]['rad_recomb'][ii] - ) # [cm3/s] + # Converts data file to ColRadPy form + data, ion, state, ratec = _conv_mr2colradpy( + FAC = FAC, + mr = mr, + react = 'rr', + ) # Formats output FAC['rates']['recomb'] = {} FAC['rates']['recomb']['recomb_transitions'] = np.vstack( (np.asarray(ion), np.asarray(state)) ).T # dim(ntrans,2), Z+1 state -> Z state - FAC['rates']['recomb']['recomb_excit'] = np.asarray( - data - ) # dim(ntrans, nt), [cm3/s] + FAC['rates']['recomb']['recomb_excit'] = np.asarray(data) # dim(ntrans, nt), [cm3/s] # Ouput return FAC @@ -721,6 +648,7 @@ def _rr_mr( def _ci_mr( FAC=None, fil=None, + verbose=None, ): # Reads data file @@ -729,61 +657,102 @@ def _ci_mr( data='ci' ) + # Converts data file to ColRadPy form + data, ion, state, ratec = _conv_mr2colradpy( + FAC = FAC, + mr = mr, + react = 'ci', + verbose=verbose, + ) + + # Formats output + FAC['rates']['ioniz'] = {} + FAC['rates']['ioniz']['ion_transitions'] = np.vstack( + (np.asarray(state), np.asarray(ion)) + ).T # dim(ntrans,2), Z state -> Z+1 state + FAC['rates']['ioniz']['ion_excit'] = np.asarray(data) # dim(ntrans, nt), [reduced rate] + if verbose == 1: + FAC['rates']['ioniz']['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, nt), [cm3/s] + + # Ouput + return FAC + +############################################################ +# +# Utilities +# +############################################################ + +# Converts data read from _read_mr to form wanted by ColRadPy +def _conv_mr2colradpy( + FAC = None, + mr = None, + react = None, + verbose = None, + ): + # Init output - state = [] - ion = [] + 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 st in mr.keys(): + for lwr in mr.keys(): # Skips temp grid - if st == 'Te_eV': + if lwr == 'Te_eV': continue # Converts indices - ind_st = np.where(FAC['lvl_indices']['FAC'] == st)[0][0] + ind_lwr = np.where(FAC['lvl_indices']['FAC'] == lwr)[0][0] # Loop over states with charge Z+1 - for ii in mr[st]['coll_ion'].keys(): + for upr in mr[lwr][lbl].keys(): # Converts indices - ind_ion = np.where(FAC['lvl_indices']['FAC'] == ii)[0][0] - state.append( - FAC['lvl_indices']['ColRadPy'][ind_st] + ind_upr = np.where(FAC['lvl_indices']['FAC'] == upr)[0][0] + st_lwr.append( + FAC['lvl_indices']['ColRadPy'][ind_lwr] ) - ion.append( - FAC['lvl_indices']['ColRadPy'][ind_ion] + st_upr.append( + FAC['lvl_indices']['ColRadPy'][ind_upr] ) data.append( - mr[st]['coll_ion'][ii] + mr[lwr][lbl][upr] ) # [cm3/s] - # Converts ionization data to adf04 reduced form - data[-1] = _conv_rate2reduct( - data = data[-1], - Te_eV = mr['Te_eV'], - ind_st = state[-1] -1, - ind_ion = ion[-1] -1, - FAC = FAC, - ) - - # Formats output - FAC['rates']['ioniz'] = {} - FAC['rates']['ioniz']['ion_transitions'] = np.vstack( - (np.asarray(state), np.asarray(ion)) - ).T # dim(ntrans,2), Z state -> Z+1 state - FAC['rates']['ioniz']['ion_excit'] = np.asarray( - data - ) # dim(ntrans, nt), [reduced rate] + 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, + ) - # Ouput - return FAC + elif react == 'ci': + # Converts ionization data to adf04 reduced form + data[-1] = _conv_rate2reduct( + data = data[-1], + Te_eV = mr['Te_eV'], + ind_st = st_lwr[-1] -1, + ind_ion = st_upr[-1] -1, + FAC = FAC, + ) -############################################################ -# -# Utilities -# -############################################################ + # Output + return data, st_upr, st_lwr, ratec # Converts rate coefficeint to adf04 reduced form def _conv_rate2reduct( @@ -817,8 +786,8 @@ def _conv_rate2reduct( def _conv_rate2upsilon( data = None, # [cm3/s], dim(ntemp,), array Te_eV = None, # [eV], dim(ntemp,) - ind_upr = None, ind_lwr = None, + ind_upr = None, FAC = None, ): @@ -886,6 +855,12 @@ def _read_ascii( elif data_fil[1][blk]['ETYPE'] == 1: out[lwr][upr]['engy'] = data_fil[1][blk]['EGRID'] + data_fil[1][blk]['TE0'] + # Stores transition energy, [eV] + out[lwr][upr]['dE'] = data_fil[1][blk]['Delta E'][trn] + + # Stores total angular momentum + out[lwr][upr]['w_upr'] = data_fil[1][blk]['upper_2J'][trn]/2 + # Stores parameters for high-energy behavior if data == 'ce': out[lwr][upr]['limit'] = np.asarray([ From c2a6d0dc41bbbbdf54e682cded8b856113d8d2c6 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Wed, 24 Jan 2024 14:59:03 -0500 Subject: [PATCH 24/46] [#18] to shorten the script and helpful readability, I broke off a bunch of general purpose functions into a utilities script --- colradpy/read_FAC.py | 297 +++---------------------------------- colradpy/read_FAC_utils.py | 295 ++++++++++++++++++++++++++++++++++++ 2 files changed, 317 insertions(+), 275 deletions(-) create mode 100644 colradpy/read_FAC_utils.py diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index adcb878..3f66a51 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -14,7 +14,6 @@ 2) Missing energy level data (zpla, zpla1) 3) Dielectronic recombination/autoionization data 4) Charge exchange data - 5) 2-photon emission data ''' @@ -24,6 +23,7 @@ import numpy as np import copy from colradpy.convolve_EEDF import convolve_EEDF +import colradpy.read_FAC_utils as utils ############################################################ # @@ -208,7 +208,7 @@ def read_FAC( ############################################################ # -# File Reading +# Transition Data Files # ############################################################ @@ -463,6 +463,12 @@ def _tr( # Output return FAC +############################################################ +# +# Cross-section Data Files +# +############################################################ + # Reads collisional excitation cross-section data files def _ce( EEDF=None, @@ -481,7 +487,7 @@ def _ce( FAC['temp_grid'] = Te*eV2K # [K], dim(nt,) # Loads date from ascii file - ce = _read_ascii( + ce = utils._read_ascii( fil = fil, data = 'ce', ) @@ -551,7 +557,7 @@ def _ce( # Convert to Upsilon form data = np.zeros((len(ratec), len(Te))) # dim(ntrans, ntemp) for tt in np.arange(len(ratec)): - data[tt,:] = _conv_rate2upsilon( + data[tt,:] = utils._conv_rate2upsilon( data = ratec[tt], Te_eV = Te, ind_upr = st_upr[tt] -1, @@ -574,6 +580,12 @@ def _ce( # Output return FAC +############################################################ +# +# fac.MaxwellRate Data Files +# +############################################################ + # Reads Maxwell-averaged collisional excitation data files def _ce_mr( FAC=None, @@ -585,7 +597,7 @@ def _ce_mr( eV2K = 11604.5 # Reads data file - mr = _read_mr( + mr = utils._read_mr( fil=fil, data='ce' ) @@ -595,7 +607,7 @@ def _ce_mr( FAC['temp_grid'] = np.asarray(mr['Te_eV'])*eV2K # [K], dim(nt,) # Converts data file to ColRadPy form - data, st_upr, st_lwr, ratec = _conv_mr2colradpy( + data, st_upr, st_lwr, ratec = utils._conv_mr2colradpy( FAC = FAC, mr = mr, react = 'ce', @@ -622,13 +634,13 @@ def _rr_mr( ): # Reads data file - mr = _read_mr( + mr = utils._read_mr( fil=fil, data='rr' ) # Converts data file to ColRadPy form - data, ion, state, ratec = _conv_mr2colradpy( + data, ion, state, ratec = utils._conv_mr2colradpy( FAC = FAC, mr = mr, react = 'rr', @@ -652,13 +664,13 @@ def _ci_mr( ): # Reads data file - mr = _read_mr( + mr = utils._read_mr( fil=fil, data='ci' ) # Converts data file to ColRadPy form - data, ion, state, ratec = _conv_mr2colradpy( + data, ion, state, ratec = utils._conv_mr2colradpy( FAC = FAC, mr = mr, react = 'ci', @@ -676,268 +688,3 @@ def _ci_mr( # Ouput return FAC - -############################################################ -# -# Utilities -# -############################################################ - -# 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[-1] = _conv_rate2reduct( - data = data[-1], - Te_eV = mr['Te_eV'], - ind_st = st_lwr[-1] -1, - ind_ion = st_upr[-1] -1, - FAC = FAC, - ) - - # Output - return data, st_upr, st_lwr, ratec - -# 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, - ): - - # Useful constants - eV2invcm = 8065.73 # [1/cm/eV] - - return data * ( - np.sqrt(np.asarray(Te_eV)/13.6058) - /2.1716e-8 - *(1 + 2*FAC['atomic']['w'][ind_upr]) - * np.exp( - abs( - FAC['atomic']['energy'][ind_upr] - -FAC['atomic']['energy'][ind_lwr] - ) - /(np.asarray(Te_eV)*eV2invcm) - ) - ) - -# Reads cross-section data files -def _read_ascii( - fil = None, - data = None, - ): - - # Reads data file - if data == 'ce': - data_fil = rfac.read_ce(fil+'a.ce') - lwr_lbl = 'lower_index' - upr_lbl = 'upper_index' - data_lbl = 'crosssection' - elif data == 'rr': - data_fil = rfac.read_rr(fil+'a.ce') - lwr_lbl = 'bound_index' - upr_lbl = 'free_index' - data_lbl = 'RR crosssection' - elif data == 'ci': - data_fil = rfac.read_rr(fil+'a.ce') - 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 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 transition energy, [eV] - out[lwr][upr]['dE'] = data_fil[1][blk]['Delta E'][trn] - - # Stores total angular momentum - out[lwr][upr]['w_upr'] = data_fil[1][blk]['upper_2J'][trn]/2 - - # Stores parameters for high-energy behavior - if data == 'ce': - out[lwr][upr]['limit'] = np.asarray([ - data_fil[1][blk]['bethe'][trn], - data_fil[1][blk]['born'][trn,0] - ]) - - # Output - return out - -# Read Maxwellian-averaged data files -def _read_mr( - fil = None, - data = None, - ): - - # Reads data file - f = open( - fil+data+'.mr', - 'r' - ) - - # Initializes output dictionary - out = {} - - # Data orginaization labels - if data == 'ce': - label = 'coll_excit' - elif data == 'rr': - label = 'rad_recomb' - elif data == '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/read_FAC_utils.py b/colradpy/read_FAC_utils.py new file mode 100644 index 0000000..3452883 --- /dev/null +++ b/colradpy/read_FAC_utils.py @@ -0,0 +1,295 @@ +''' + +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 + +############################################################ +# +# Converting Data to ColRadPy format +# +############################################################ + +# 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[-1] = _conv_rate2reduct( + data = data[-1], + Te_eV = mr['Te_eV'], + ind_st = st_lwr[-1] -1, + ind_ion = st_upr[-1] -1, + FAC = FAC, + ) + + # 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, + ): + + # Useful constants + eV2invcm = 8065.73 # [1/cm/eV] + + return data * ( + np.sqrt(np.asarray(Te_eV)/13.6058) + /2.1716e-8 + *(1 + 2*FAC['atomic']['w'][ind_upr]) + * np.exp( + abs( + FAC['atomic']['energy'][ind_upr] + -FAC['atomic']['energy'][ind_lwr] + ) + /(np.asarray(Te_eV)*eV2invcm) + ) + ) + +############################################################ +# +# Reading Data Files +# +############################################################ + +# Reads cross-section data files +def _read_ascii( + fil = None, + data = None, + ): + + # Reads data file + if data == 'ce': + data_fil = rfac.read_ce(fil+'a.ce') + lwr_lbl = 'lower_index' + upr_lbl = 'upper_index' + data_lbl = 'crosssection' + elif data == 'rr': + data_fil = rfac.read_rr(fil+'a.ce') + lwr_lbl = 'bound_index' + upr_lbl = 'free_index' + data_lbl = 'RR crosssection' + elif data == 'ci': + data_fil = rfac.read_rr(fil+'a.ce') + 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 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 transition energy, [eV] + out[lwr][upr]['dE'] = data_fil[1][blk]['Delta E'][trn] + + # Stores total angular momentum + out[lwr][upr]['w_upr'] = data_fil[1][blk]['upper_2J'][trn]/2 + + # Stores parameters for high-energy behavior + if data == 'ce': + out[lwr][upr]['limit'] = np.asarray([ + data_fil[1][blk]['bethe'][trn], + data_fil[1][blk]['born'][trn,0] + ]) + + # Output + return out + +# Read Maxwellian-averaged data files +def _read_mr( + fil = None, + data = None, + ): + + # Reads data file + f = open( + fil+data+'.mr', + 'r' + ) + + # Initializes output dictionary + out = {} + + # Data orginaization labels + if data == 'ce': + label = 'coll_excit' + elif data == 'rr': + label = 'rad_recomb' + elif data == '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 From 0d22557b4f3bb3bb5f588c2d73d9daeded67afb8 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Wed, 24 Jan 2024 15:19:49 -0500 Subject: [PATCH 25/46] [#18] made functions to use FAC excit cross-sections for convolution --- colradpy/read_FAC.py | 94 ++++++----------------------- colradpy/read_FAC_utils.py | 119 +++++++++++++++++++++++++++++++++---- 2 files changed, 126 insertions(+), 87 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 3f66a51..7b5cccf 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -22,7 +22,6 @@ import os import numpy as np import copy -from colradpy.convolve_EEDF import convolve_EEDF import colradpy.read_FAC_utils as utils ############################################################ @@ -487,82 +486,23 @@ def _ce( FAC['temp_grid'] = Te*eV2K # [K], dim(nt,) # Loads date from ascii file - ce = utils._read_ascii( + XSdata = utils._read_ascii( fil = fil, - data = 'ce', + react = 'ce', ) - # Init output - st_lwr = [] - st_upr = [] - ratec = [] - if verbose == 1: - XS = [] - engy = [] - - # Loop over lower states - for lwr in ce.keys(): - # Converts indices - ind_lwr = np.where(FAC['lvl_indices']['FAC'] == lwr)[0][0] - - # Loop over upper states - for upr in ce[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] - ) - - # Calculates Rate coefficient data, [cm3/s], dim(ntrans, ntemp) - if verbose == 1: - ( - tmp_ratec, - tmp_XS, - tmp_engy - )= convolve_EEDF( - EEDF = EEDF, - Te=Te, - XS = ce[lwr][upr]['XS'], - engyXS = ce[lwr][upr]['engy'], - m = 0, - dE = np.asarray([ce[lwr][upr]['dE']]), - Bethe = ce[lwr][upr]['limit'][None,:], - w_upr = np.asarray([ce[lwr][upr]['w_upr']]), - verbose=verbose, - use_rel = True, - ) - ratec.append(tmp_ratec[:,0]) - XS.append(tmp_XS[:,0]) - engy.append(tmp_engy[:,0]) - - else: - ratec.append( - convolve_EEDF( - EEDF = EEDF, - Te=Te, - XS = ce[lwr][upr]['XS'], - engyXS = ce[lwr][upr]['engy'], - m = 0, - dE = np.asarray([ce[lwr][upr]['dE']]), - Bethe = ce[lwr][upr]['limit'][None,:], - w_upr = np.asarray([ce[lwr][upr]['w_upr']]), - verbose=verbose, - use_rel = True, - )[:,0] - ) - - # Convert to Upsilon form - data = np.zeros((len(ratec), len(Te))) # dim(ntrans, ntemp) - for tt in np.arange(len(ratec)): - data[tt,:] = utils._conv_rate2upsilon( - data = ratec[tt], - Te_eV = Te, - ind_upr = st_upr[tt] -1, - ind_lwr = st_lwr[tt] -1, + # 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 = 'ce', + verbose = verbose, ) # Formats output @@ -570,7 +510,7 @@ def _ce( FAC['rates']['excit']['col_transitions'] = np.vstack( (np.asarray(st_upr), np.asarray(st_lwr)) ).T # dim(ntrans,2), (upr,lwr) states in ColRadPy indices - FAC['rates']['excit']['col_excit'] = data # dim(ntrans, ntemp), [upsilon] + FAC['rates']['excit']['col_excit'] = np.asarray(data) # dim(ntrans, ntemp), [upsilon] if verbose == 1: FAC['rates']['excit']['col_trans_unfill'] = FAC['rates']['excit']['col_transitions'].copy() FAC['rates']['excit']['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, ntemp), [cm3/s] @@ -599,7 +539,7 @@ def _ce_mr( # Reads data file mr = utils._read_mr( fil=fil, - data='ce' + react='ce' ) # Saves temperature grid data @@ -636,7 +576,7 @@ def _rr_mr( # Reads data file mr = utils._read_mr( fil=fil, - data='rr' + react='rr' ) # Converts data file to ColRadPy form @@ -666,7 +606,7 @@ def _ci_mr( # Reads data file mr = utils._read_mr( fil=fil, - data='ci' + react='ci' ) # Converts data file to ColRadPy form diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py index 3452883..56ff065 100644 --- a/colradpy/read_FAC_utils.py +++ b/colradpy/read_FAC_utils.py @@ -22,9 +22,108 @@ # ############################################################ +# 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] + ) + + # 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']]), + Bethe = XSdata[lwr][upr]['limit'][None,:], + w_upr = np.asarray([XSdata[lwr][upr]['w_upr']]), + verbose=verbose, + use_rel = True, + ) + data.append(tmp_data[:,0]) + XS.append(tmp_XS[:,0]) + engy.append(tmp_engy[:,0]) + + 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']]), + Bethe = XSdata[lwr][upr]['limit'][None,:], + w_upr = np.asarray([XSdata[lwr][upr]['w_upr']]), + verbose=verbose, + use_rel = True, + )[:,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[-1] = _conv_rate2reduct( + data = data[-1], + Te_eV = Te, + ind_st = st_lwr[-1] -1, + ind_ion = st_upr[-1] -1, + FAC = FAC, + ) + + # 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, + FAC = None, mr = None, react = None, verbose = None, @@ -162,21 +261,21 @@ def _conv_rate2upsilon( # Reads cross-section data files def _read_ascii( fil = None, - data = None, + react = None, ): # Reads data file - if data == 'ce': + if react == 'ce': data_fil = rfac.read_ce(fil+'a.ce') lwr_lbl = 'lower_index' upr_lbl = 'upper_index' data_lbl = 'crosssection' - elif data == 'rr': + elif react == 'rr': data_fil = rfac.read_rr(fil+'a.ce') lwr_lbl = 'bound_index' upr_lbl = 'free_index' data_lbl = 'RR crosssection' - elif data == 'ci': + elif react == 'ci': data_fil = rfac.read_rr(fil+'a.ce') lwr_lbl = 'bound_index' upr_lbl = 'free_index' @@ -214,7 +313,7 @@ def _read_ascii( out[lwr][upr]['w_upr'] = data_fil[1][blk]['upper_2J'][trn]/2 # Stores parameters for high-energy behavior - if data == 'ce': + if react == 'ce': out[lwr][upr]['limit'] = np.asarray([ data_fil[1][blk]['bethe'][trn], data_fil[1][blk]['born'][trn,0] @@ -226,7 +325,7 @@ def _read_ascii( # Read Maxwellian-averaged data files def _read_mr( fil = None, - data = None, + react = None, ): # Reads data file @@ -239,11 +338,11 @@ def _read_mr( out = {} # Data orginaization labels - if data == 'ce': + if react == 'ce': label = 'coll_excit' - elif data == 'rr': + elif react == 'rr': label = 'rad_recomb' - elif data == 'ci': + elif react == 'ci': label = 'coll_ion' # Loop over lines From c9248a50f47955feeaee2d64d66fdbabb96e25db Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Wed, 24 Jan 2024 15:25:28 -0500 Subject: [PATCH 26/46] [#18] added high-energy fill to the cross-section data by Bethe asymptotic limit for excit and ioniz --- colradpy/convolve_EEDF.py | 115 ++++++++++++++++++++++++++++++------- colradpy/read_FAC_utils.py | 2 + 2 files changed, 96 insertions(+), 21 deletions(-) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index adf0ea5..b48937f 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -28,6 +28,11 @@ def convolve_EEDF( m = 0, # 0 or 1 dE = 0, # [eV] DC_flag = False, + Bethe = None, # (optional), Bethe asymptotic behavior + w_upr = None, # (optional) Upper level statistical weight + verbose = 1, + use_rel = True, # flag to use relativistic corrections + react = None, # reaction type in ['ce', 'rr', 'ci'] ): ''' INPUTS: @@ -77,6 +82,7 @@ def convolve_EEDF( if EEDF == 'Maxwellian': EEDF, engyEEDF = _get_Max( Te=Te, + use_rel=use_rel, ) # dim(ngrid, ntemp), ([1/eV], [eV]) else: print('NON-MAXWELLIAN ELECTRONS NOT IMPLEMENTED YET') @@ -91,6 +97,7 @@ def convolve_EEDF( engyXS=engyXS, m=m, dE=dE, + use_rel=use_rel, ) # [cm3/s], dim(ntemp,ntrans) # Performs numerical integration @@ -102,6 +109,11 @@ def convolve_EEDF( engyXS=engyXS, m=m, dE=dE, + Bethe = Bethe, + w_upr = w_upr, + verbose = verbose, + use_rel=use_rel, + react=react, ) # [cm3/s], dim(ntemp,ntrans) @@ -119,6 +131,7 @@ def _calc_DC( engyXS = None, # [eV], dim(ntrans,) m = None, dE = None, # [eV], dim(ntrans,) + use_rel=None, ): # Useful constants @@ -134,11 +147,16 @@ def _calc_DC( elif m == 1: engy_tmp = engyXS[nn] + dE[nn] - # Incident electron velocity, relativistic form - vel = cnt.c *100 *np.sqrt( - 1 - - 1/(engy_tmp *eV2electron +1)**2 - ) # [cm/s] + # Incident electron velocity + if use_rel: # relativistic form + vel = cnt.c *100 *np.sqrt( + 1 - + 1/(engyEEDF[:,tt] *eV2electron +1)**2 + ) # [cm/s], dim(ngrid,) + else: # classical form + vel = cnt.c *100 *np.sqrt( + 2*engyEEDF[:,tt] *eV2electron + ) # [cm/s], dim(ngrid,) # Loop over temperatures for tt in np.arange(EEDF.shape[1]): @@ -164,6 +182,11 @@ def _calc_ratec( engyXS = None, # [eV], dim(nE, ntrans) m = None, dE = None, # [eV], dim(ntrans,) + Bethe = None, # dim(ntrans,2) + w_upr = None, + verbose = None, + use_rel = None, + react = None, ): # Useful constants @@ -171,14 +194,26 @@ def _calc_ratec( # 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, relativistic form - vel = cnt.c *100 *np.sqrt( - 1 - - 1/(engyEEDF[:,tt] *eV2electron +1)**2 - ) # [cm/s], dim(ngrid,) + # Incident electron velocity + if use_rel: # relativistic form + vel = cnt.c *100 *np.sqrt( + 1 - + 1/(engyEEDF[:,tt] *eV2electron +1)**2 + ) # [cm/s], dim(ngrid,) + else: # classical form + vel = cnt.c *100 *np.sqrt( + 2*engyEEDF[:,tt] *eV2electron + ) # [cm/s], dim(ngrid,) # Loop over transitions for nn in np.arange(XS.shape[1]): @@ -195,22 +230,60 @@ def _calc_ratec( fill_value = (-1e5,-1e5) )(np.log10(engyEEDF[:,tt])) # dim(ngrid,), [cm2] + # Fill values with Bethe asymptotic behavior if available + if Bethe is not None: + # FAC cross-section data isn't really trustworthy about 10x transition energy + tol = 10 + if engy_tmp[-1] >= tol*dE[nn]: + indE = np.where(engyEEDF[:,tt] > tol*dE[nn])[0] + else: + indE = np.where(engyEEDF[:,tt] > engy_tmp[-1])[0] + if len(indE) == 0: + continue + # Asymptotic form of collision strength for (first term) optically + # allowed and (second term) forbidden transitions + if react in ['ce', 'ci']: + omega = ( + Bethe[nn,0] + * np.log(engyEEDF[indE,tt]/dE[nn]) + + Bethe[nn,1] + ) + + XS_tmp[indE] = ( + omega + *np.pi *cnt.physical_constants['Bohr radius'][0]**2 + * 13.6 /engyEEDF[indE,tt] + /(1 +2*w_upr[nn]) + ) *1e4 + + # Account for the different normalization of bound & free states + if react == 'ci': + XS_tmp[indE] /= np.pi + # 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) - return ratec + if verbose == 1: + return ratec, XS_out, engy_out + else: + return ratec # Calculate Maxwellian energy distribution function def _get_Max( Te=None, # [eV], dim(ntemp,) # Energy grid settings - ngrid = int(5e2), # energy grid resolution + 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 @@ -231,16 +304,8 @@ def _get_Max( ngrid ) # [eV] - # low-energy form - if Te[tt] <= 10e3: - EEDF[:,tt] = ( - 2*np.sqrt(engy[:,tt]/np.pi) - * Te[tt]**(-3/2) - * np.exp(-engy[:,tt]/Te[tt]) - ) # dim(ngrid, ntemp); [1/eV] - # relativistic form - else: + if Te[tt] >= 10e3 and use_rel: # Factors theta = Te[tt] * eV2electron gamma = 1 + engy[:,tt] * eV2electron @@ -254,6 +319,14 @@ def _get_Max( *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 diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py index 56ff065..95ecdd3 100644 --- a/colradpy/read_FAC_utils.py +++ b/colradpy/read_FAC_utils.py @@ -74,6 +74,7 @@ def _conv_ascii2colradpy( w_upr = np.asarray([XSdata[lwr][upr]['w_upr']]), verbose=verbose, use_rel = True, + react = react, ) data.append(tmp_data[:,0]) XS.append(tmp_XS[:,0]) @@ -92,6 +93,7 @@ def _conv_ascii2colradpy( w_upr = np.asarray([XSdata[lwr][upr]['w_upr']]), verbose=verbose, use_rel = True, + react = react, )[:,0] ) From 30e31cb74df5f9fd0f095bd4d673aa7223a665e6 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Wed, 24 Jan 2024 15:35:47 -0500 Subject: [PATCH 27/46] [#18] added functions to convolve rad recomb and coll ioniz cross-sections --- colradpy/read_FAC.py | 94 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 92 insertions(+), 2 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 7b5cccf..75d76fb 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -138,6 +138,8 @@ def read_FAC( FAC=FAC, fil=fil, EEDF=EEDF, + Te=Te, + verbose=verbose, ) elif 'rr.mr' in physics: FAC = _rr_mr( @@ -156,6 +158,7 @@ def read_FAC( FAC=FAC, fil=fil, EEDF=EEDF, + Te=Te, ) elif 'ai.mr' in physics: FAC = _ai_mr( @@ -165,10 +168,12 @@ def read_FAC( # Collisional ionization if 'ci' in physics: - FAC = _ai( + FAC = _ci( FAC=FAC, fil=fil, EEDF=EEDF, + Te=Te, + verbose=verbose, ) elif 'ci.mr' in physics: FAC = _ci_mr( @@ -474,7 +479,6 @@ def _ce( Te=None, # [eV], dim(ntemp,) FAC=None, fil=None, - trans_FAC=None, vebose = None, ): @@ -520,6 +524,92 @@ def _ce( # Output return FAC +# Reads radiative recombination cross-section data files +def _rr( + EEDF=None, + Te=None, # [eV], dim(ntemp,) + FAC=None, + fil=None, + verbose=None, + ): + + # Loads date from ascii file + XSdata = utils._read_ascii( + fil = fil, + react = 'rr', + ) + + # Performs EEDF convolution to cross-sections + ( + data, + ion, state, + ratec, XS, engy + ) = utils._conv_ascii2colradpy( + FAC = FAC, + XSdata = XSdata, + EEDF = EEDF, + Te = Te, + react = 'rr', + verbose = verbose, + ) + + # Formats output + FAC['rates']['recomb'] = {} + FAC['rates']['recomb']['recomb_transitions'] = np.vstack( + (np.asarray(ion), np.asarray(state)) + ).T # dim(ntrans,2), Z+1 state -> Z state + FAC['rates']['recomb']['recomb_excit'] = np.asarray(data) # dim(ntrans, nt), [cm3/s] + if verbose == 1: + FAC['rates']['recomb']['XS_cm2'] = np.asarray(XS) # dim(ntrans, nE), [cm2] + FAC['rates']['recomb']['engy_eV'] = np.asarray(engy) # dim(ntrans, nE), [eV] + + # Output + return FAC + + +# Reads collisional ionization cross-section data files +def _ci( + EEDF=None, + Te=None, # [eV], dim(ntemp,) + FAC=None, + fil=None, + verbose=None, + ): + + # Loads date from ascii file + XSdata = utils._read_ascii( + fil = fil, + react = 'ci', + ) + + # Performs EEDF convolution to cross-sections + ( + data, + ion, state, + ratec, XS, engy + ) = utils._conv_ascii2colradpy( + FAC = FAC, + XSdata = XSdata, + EEDF = EEDF, + Te = Te, + react = 'ci', + verbose = verbose, + ) + + # Formats output + FAC['rates']['ioniz'] = {} + FAC['rates']['ioniz']['ion_transitions'] = np.vstack( + (np.asarray(state), np.asarray(ion)) + ).T # dim(ntrans,2), Z state -> Z+1 state + FAC['rates']['ioniz']['ion_excit'] = np.asarray(data) # dim(ntrans, nt), [reduced rate] + if verbose == 1: + FAC['rates']['ioniz']['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, nt), [cm3/s] + FAC['rates']['ioniz']['XS_cm2'] = np.asarray(XS) # dim(ntrans, nE), [cm2] + FAC['rates']['ioniz']['engy_eV'] = np.asarray(engy) # dim(ntrans, nE), [eV] + + # Output + return FAC + ############################################################ # # fac.MaxwellRate Data Files From 2a381765ecb60c196618085178afe58d29d6ce78 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Wed, 24 Jan 2024 16:21:10 -0500 Subject: [PATCH 28/46] [#18] init high-energy asymptotic cross-section for rad recomb and coll ioniz --- colradpy/convolve_EEDF.py | 106 +++++++++++++++++++++++++------------ colradpy/read_FAC.py | 4 +- colradpy/read_FAC_utils.py | 53 +++++++++++++++---- 3 files changed, 116 insertions(+), 47 deletions(-) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index b48937f..28ad4ee 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -28,8 +28,9 @@ def convolve_EEDF( m = 0, # 0 or 1 dE = 0, # [eV] DC_flag = False, - Bethe = None, # (optional), Bethe asymptotic behavior - w_upr = None, # (optional) Upper level statistical weight + 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 verbose = 1, use_rel = True, # flag to use relativistic corrections react = None, # reaction type in ['ce', 'rr', 'ci'] @@ -109,8 +110,9 @@ def convolve_EEDF( engyXS=engyXS, m=m, dE=dE, - Bethe = Bethe, + limit = limit, w_upr = w_upr, + w_lwr = w_lwr, verbose = verbose, use_rel=use_rel, react=react, @@ -179,11 +181,12 @@ 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) + engyXS = None, # [eV], dim(nE, ntrans), incident electron energy m = None, dE = None, # [eV], dim(ntrans,) - Bethe = None, # dim(ntrans,2) + limit = None, # dim(ntrans,2) if ce or dim(ntrans,4) if rr, ci w_upr = None, + w_lwr = None, verbose = None, use_rel = None, react = None, @@ -230,35 +233,17 @@ def _calc_ratec( fill_value = (-1e5,-1e5) )(np.log10(engyEEDF[:,tt])) # dim(ngrid,), [cm2] - # Fill values with Bethe asymptotic behavior if available - if Bethe is not None: - # FAC cross-section data isn't really trustworthy about 10x transition energy - tol = 10 - if engy_tmp[-1] >= tol*dE[nn]: - indE = np.where(engyEEDF[:,tt] > tol*dE[nn])[0] - else: - indE = np.where(engyEEDF[:,tt] > engy_tmp[-1])[0] - if len(indE) == 0: - continue - # Asymptotic form of collision strength for (first term) optically - # allowed and (second term) forbidden transitions - if react in ['ce', 'ci']: - omega = ( - Bethe[nn,0] - * np.log(engyEEDF[indE,tt]/dE[nn]) - + Bethe[nn,1] - ) - - XS_tmp[indE] = ( - omega - *np.pi *cnt.physical_constants['Bohr radius'][0]**2 - * 13.6 /engyEEDF[indE,tt] - /(1 +2*w_upr[nn]) - ) *1e4 - - # Account for the different normalization of bound & free states - if react == 'ci': - XS_tmp[indE] /= np.pi + # Fill values with high-energy asymptotic behavior if available + 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], + ) # Preforms integration ratec[tt,nn] = np.trapz( @@ -276,6 +261,59 @@ def _calc_ratec( else: return ratec +############################################################ +# +# Calculation +# +############################################################ + +# Calculate high-energy asymptotic behavior of cross-section +def _get_limit( + XS_tmp = None, + engy_tmp = None, + engyEEDF = None, + dE = None, + limit = None, + w_upr = None, + w_lwr = None, + ): + + # 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: + return XS_tmp + + if react == 'ce': + # 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] + ) + + XS_tmp[indE] = ( + omega + *np.pi *cnt.physical_constants['Bohr radius'][0]**2 + * 13.6 /engyEEDF[indE] + /(1 +2*w_upr) + ) *1e4 + + elif react == 'rr': + blah = 0 + + elif react == 'ci': + blah = 0 + + # Output + return XS_tmp + # Calculate Maxwellian energy distribution function def _get_Max( Te=None, # [eV], dim(ntemp,) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 75d76fb..6b0e76b 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -74,9 +74,9 @@ def read_FAC( 'en', # ASCII-format Energy levels 'tr', # ASCII-format Einstein coefficients 'ce', # ASCII-format Collisional excitation - #'rr', # ASCII-format Radiative recombination + 'rr', # ASCII-format Radiative recombination #'ai', # ASCII-format Autoionization/dielectronic recombination - #'ci', # ASCII-format Collision ionization + 'ci', # ASCII-format Collision ionization ] ######## -------- Reads data -------- ######## diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py index 95ecdd3..1175111 100644 --- a/colradpy/read_FAC_utils.py +++ b/colradpy/read_FAC_utils.py @@ -70,8 +70,9 @@ def _conv_ascii2colradpy( engyXS = XSdata[lwr][upr]['engy'], m = 0, dE = np.asarray([XSdata[lwr][upr]['dE']]), - Bethe = XSdata[lwr][upr]['limit'][None,:], + limit = XSdata[lwr][upr]['limit'][None,:], w_upr = np.asarray([XSdata[lwr][upr]['w_upr']]), + w_lwr = np.asarray([XSdata[lwr][upr]['w_lwr']]), verbose=verbose, use_rel = True, react = react, @@ -91,6 +92,7 @@ def _conv_ascii2colradpy( dE = np.asarray([XSdata[lwr][upr]['dE']]), Bethe = XSdata[lwr][upr]['limit'][None,:], w_upr = np.asarray([XSdata[lwr][upr]['w_upr']]), + w_lwr = np.asarray([XSdata[lwr][upr]['w_lwr']]), verbose=verbose, use_rel = True, react = react, @@ -302,25 +304,54 @@ def _read_ascii( # Stores cross-section data, [cm2] out[lwr][upr]['XS'] = data_fil[1][blk][data_lbl][trn,:]*1e-20 - # 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 transition energy, [eV] out[lwr][upr]['dE'] = data_fil[1][blk]['Delta E'][trn] - # Stores total angular momentum - out[lwr][upr]['w_upr'] = data_fil[1][blk]['upper_2J'][trn]/2 - - # Stores parameters for high-energy behavior 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 From df6257417a3e8c90e223a9baa2211dd17124f85e Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 25 Jan 2024 09:56:49 -0500 Subject: [PATCH 29/46] [#18] reduced redundancy in convolve_EEDF. Simplified high-energy cross-section for excit, ioniz --- colradpy/convolve_EEDF.py | 112 ++++++++++++++++++++++++++------------ 1 file changed, 78 insertions(+), 34 deletions(-) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index 28ad4ee..11017b4 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -136,9 +136,6 @@ def _calc_DC( use_rel=None, ): - # Useful constants - eV2electron = cnt.e /(cnt.m_e*cnt.c**2) # [1/eV] - # Init output ratec = np.zeros((EEDF.shape[1], len(XS))) # [cm3/s], dim(ntemp,ntrans) @@ -150,15 +147,10 @@ def _calc_DC( engy_tmp = engyXS[nn] + dE[nn] # Incident electron velocity - if use_rel: # relativistic form - vel = cnt.c *100 *np.sqrt( - 1 - - 1/(engyEEDF[:,tt] *eV2electron +1)**2 - ) # [cm/s], dim(ngrid,) - else: # classical form - vel = cnt.c *100 *np.sqrt( - 2*engyEEDF[:,tt] *eV2electron - ) # [cm/s], dim(ngrid,) + vel = _get_vel( + E_inc = engyEEDF[:,tt], + use_rel = use_rel, + ) # [cm/s], dim(ngrid,) # Loop over temperatures for tt in np.arange(EEDF.shape[1]): @@ -192,9 +184,6 @@ def _calc_ratec( react = None, ): - # Useful constants - eV2electron = cnt.e /(cnt.m_e*cnt.c**2) # [1/eV] - # Init output ratec = np.zeros((EEDF.shape[1], XS.shape[1])) # [cm3/s], dim(ntemp,ntrans) if verbose == 1: @@ -208,15 +197,10 @@ def _calc_ratec( # Loop over temperatures for tt in np.arange(EEDF.shape[1]): # Incident electron velocity - if use_rel: # relativistic form - vel = cnt.c *100 *np.sqrt( - 1 - - 1/(engyEEDF[:,tt] *eV2electron +1)**2 - ) # [cm/s], dim(ngrid,) - else: # classical form - vel = cnt.c *100 *np.sqrt( - 2*engyEEDF[:,tt] *eV2electron - ) # [cm/s], dim(ngrid,) + 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]): @@ -243,6 +227,7 @@ def _calc_ratec( limit = limit[nn,:], w_upr = w_upr[nn], w_lwr = w_lwr[nn], + react = react, ) # Preforms integration @@ -268,16 +253,23 @@ def _calc_ratec( ############################################################ # Calculate high-energy asymptotic behavior of cross-section +## NOTE: So far, FAC-specific def _get_limit( XS_tmp = None, - engy_tmp = None, - engyEEDF = None, - dE = 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, + 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: @@ -290,26 +282,54 @@ def _get_limit( 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] = ( - omega - *np.pi *cnt.physical_constants['Bohr radius'][0]**2 - * 13.6 /engyEEDF[indE] - /(1 +2*w_upr) - ) *1e4 + np.pi *omega + / k02 + / (1 +2*w_upr) + ) * a02 # [cm2] elif react == 'rr': blah = 0 elif react == 'ci': - blah = 0 + # 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_upr) + ) * a02 # [cm2] # Output return XS_tmp @@ -368,3 +388,27 @@ def _get_Max( # 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,) From e11cbd19c4f57b58387e4266aa9d8194cf51a5bd Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 25 Jan 2024 10:01:32 -0500 Subject: [PATCH 30/46] [#18] fill missing cross-section data right at threshold --- colradpy/convolve_EEDF.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index 11017b4..da9ce21 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -217,6 +217,13 @@ def _calc_ratec( fill_value = (-1e5,-1e5) )(np.log10(engyEEDF[:,tt])) # dim(ngrid,), [cm2] + # Fill missing data right at threshold + 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 if limit is not None: XS_tmp = _get_limit( From 7d04e8693f1b7a1ed6705a662b1efbad50400c4e Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 25 Jan 2024 10:44:35 -0500 Subject: [PATCH 31/46] [#18] reduced redundancy in loading FAC data --- colradpy/colradpy_class.py | 3 +- colradpy/read_FAC.py | 497 +++++++++++++++---------------------- 2 files changed, 200 insertions(+), 300 deletions(-) diff --git a/colradpy/colradpy_class.py b/colradpy/colradpy_class.py index c38d278..52ee0e7 100644 --- a/colradpy/colradpy_class.py +++ b/colradpy/colradpy_class.py @@ -349,7 +349,8 @@ def populate_data(self,fil): Zele=self.data['user']['Zele'], fil=fil, EEDF=self.data['user']['EEDF'], - physics=self.data['user']['atomic_physics'], + 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 diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 6b0e76b..90eb4fb 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -39,7 +39,7 @@ def read_FAC( fil = None, # Common path to FAC files, excluding physics extentions # Physics controls EEDF = None, # if None -> assumes Maxwell-averages rates from pfac.fac.MaxwellRate - physics = None, # if None -> looks for all file suffixes + reacts = None, # if None -> looks for all file suffixes Te = None, # if not using MaxwellRate files, [eV], dim(ntemp,) verbose = 1, ): @@ -56,28 +56,16 @@ def read_FAC( print('NON-MAXWELLIAN ELECTRONS NOT IMPLEMENTED YET') sys.exit(1) - # FAC data file suffixes to search for - if physics == 'incl_all': - # If already Maxwell-averaged rate coefficients - if use_mr: - physics = [ - 'en', # ASCII-format Energy levels - 'tr', # ASCII-format Einstein coefficients - 'ce.mr', # Collisional excitation - 'rr.mr', # Radiative recombination - #'ai', # Autoionization/dielectronic recombination - 'ci.mr', # Collision ionization - ] - # If using cross-section files - else: - physics = [ - 'en', # ASCII-format Energy levels - 'tr', # ASCII-format Einstein coefficients - 'ce', # ASCII-format Collisional excitation - 'rr', # ASCII-format Radiative recombination - #'ai', # ASCII-format Autoionization/dielectronic recombination - 'ci', # ASCII-format Collision ionization - ] + # 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 -------- ######## @@ -86,7 +74,7 @@ def read_FAC( FAC['rates'] = {} # Energy levels - if 'en' in physics: + if 'en' in reacts: FAC = _en( FAC=FAC, fil=fil, @@ -99,28 +87,29 @@ def read_FAC( print('NEED TO INCLUDE ENERGY LEVEL DATA IN MODELING!!!') sys.exit(1) - # Collisional excitation - if 'ce' in physics: - FAC = _ce( + # 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, - EEDF=EEDF, - Te=Te, verbose=verbose, + reacts=reacts, ) - elif 'ce.mr' in physics: - FAC = _ce_mr( + else: + # Use cross-section data files + FAC = _get_xs( FAC=FAC, fil=fil, + EEDF=EEDF, + Te=Te, verbose=verbose, + reacts=reacts, ) - # Error check - else: - print('NEED TO INCLUDE COLLISIONAL EXCITATION DATA IN MODELING!!!') - sys.exit(1) # Einstein coefficients - if 'tr' in physics: + if 'tr' in reacts: FAC = _tr( FAC=FAC, fil=fil, @@ -132,63 +121,17 @@ def read_FAC( print('NEED TO INCLUDE EINSTEIN COEFFICIENT DATA IN MODELING!!!') sys.exit(1) - # Radiative recombination - if 'rr' in physics: - FAC = _rr( - FAC=FAC, - fil=fil, - EEDF=EEDF, - Te=Te, - verbose=verbose, - ) - elif 'rr.mr' in physics: - FAC = _rr_mr( - FAC=FAC, - fil=fil, - ) - # If empty - else: - FAC['rates']['recomb'] = {} - FAC['rates']['recomb']['recomb_transitions'] = np.asarray([]) - FAC['rates']['recomb']['recomb_excit'] = np.asarray([]) - # Autoionization/dielectronic recombination - if 'ai' in physics: + if 'ai' in reacts: FAC = _ai( FAC=FAC, fil=fil, EEDF=EEDF, Te=Te, ) - elif 'ai.mr' in physics: - FAC = _ai_mr( - FAC=FAC, - fil=fil, - ) - - # Collisional ionization - if 'ci' in physics: - FAC = _ci( - FAC=FAC, - fil=fil, - EEDF=EEDF, - Te=Te, - verbose=verbose, - ) - elif 'ci.mr' in physics: - FAC = _ci_mr( - FAC=FAC, - fil=fil, - verbose=verbose, - ) - # If empty - else: - FAC['rates']['ioniz'] = {} - FAC['rates']['ioniz']['ion_transitions'] = np.asarray([]) - FAC['rates']['ioniz']['ion_excit'] = np.asarray([]) # Charge exchange - if 'cx' in physics: + if 'cx' in reacts: print('CHARGE EXCHANGE IS NOT IMPLEMENTED YET!!!') # If empty else: @@ -469,252 +412,208 @@ def _tr( ############################################################ # -# Cross-section Data Files +# Data Loading & Prep # ############################################################ -# Reads collisional excitation cross-section data files -def _ce( +# Reads cross-section data files +def _get_xs( EEDF=None, Te=None, # [eV], dim(ntemp,) FAC=None, fil=None, vebose = None, + reacts = None, ): + # Error check + if 'ce' not in reacts: + print('NEED TO INCLUDE COLLISIONAL EXCITATION DATA IN MODELING!!!') + sys.exit(1) + + # 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 - # Saves temperature grid data - # NOTE: ColRadPy assumes Te is in Kelvin within the data file (per ADF04 standard) - FAC['temp_grid'] = Te*eV2K # [K], dim(nt,) + # 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 = 'ce', - ) - - # 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 = 'ce', - verbose = verbose, + # Loads date from ascii file + XSdata = utils._read_ascii( + fil = fil, + react = react, ) - # Formats output - FAC['rates']['excit'] = {} - FAC['rates']['excit']['col_transitions'] = np.vstack( - (np.asarray(st_upr), np.asarray(st_lwr)) - ).T # dim(ntrans,2), (upr,lwr) states in ColRadPy indices - FAC['rates']['excit']['col_excit'] = np.asarray(data) # dim(ntrans, ntemp), [upsilon] - if verbose == 1: - FAC['rates']['excit']['col_trans_unfill'] = FAC['rates']['excit']['col_transitions'].copy() - FAC['rates']['excit']['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, ntemp), [cm3/s] - FAC['rates']['excit']['XS_cm2'] = np.asarray(XS) # dim(ntrans, nE), [cm2] - FAC['rates']['excit']['engy_eV'] = np.asarray(engy) # dim(ntrans, nE), [eV] + # 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), [cm2] + FAC['rates'][lbl]['engy_eV'] = np.asarray(engy) # dim(ntrans, nE), [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 radiative recombination cross-section data files -def _rr( - EEDF=None, - Te=None, # [eV], dim(ntemp,) - FAC=None, - fil=None, - verbose=None, +# Reads Maxwell-averaged data files +def _get_mr( + FAC = None, + fil = None, + verbose = None, + reacts = None, ): - # Loads date from ascii file - XSdata = utils._read_ascii( - fil = fil, - react = 'rr', - ) - - # Performs EEDF convolution to cross-sections - ( - data, - ion, state, - ratec, XS, engy - ) = utils._conv_ascii2colradpy( - FAC = FAC, - XSdata = XSdata, - EEDF = EEDF, - Te = Te, - react = 'rr', - verbose = verbose, - ) + # Error check + if 'ce' not in reacts: + print('NEED TO INCLUDE COLLISIONAL EXCITATION DATA IN MODELING!!!') + sys.exit(1) - # Formats output - FAC['rates']['recomb'] = {} - FAC['rates']['recomb']['recomb_transitions'] = np.vstack( - (np.asarray(ion), np.asarray(state)) - ).T # dim(ntrans,2), Z+1 state -> Z state - FAC['rates']['recomb']['recomb_excit'] = np.asarray(data) # dim(ntrans, nt), [cm3/s] - if verbose == 1: - FAC['rates']['recomb']['XS_cm2'] = np.asarray(XS) # dim(ntrans, nE), [cm2] - FAC['rates']['recomb']['engy_eV'] = np.asarray(engy) # dim(ntrans, nE), [eV] + # 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([]) - # Output - return FAC + # Useful constants + eV2K = 11604.5 + # Loop over reactions + for react in reacts: + # Skip + if react not in ['ce', 'rr', 'ci']: + continue -# Reads collisional ionization cross-section data files -def _ci( - EEDF=None, - Te=None, # [eV], dim(ntemp,) - FAC=None, - fil=None, - verbose=None, - ): + # Reads data file + mr = utils._read_mr( + fil = fil, + react = react, + ) - # Loads date from ascii file - XSdata = utils._read_ascii( - fil = fil, - react = 'ci', - ) + # 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,) - # Performs EEDF convolution to cross-sections - ( - data, - ion, state, - ratec, XS, engy - ) = utils._conv_ascii2colradpy( + # Converts data file to ColRadPy form + data, st_upr, st_lwr, ratec = utils._conv_mr2colradpy( FAC = FAC, - XSdata = XSdata, - EEDF = EEDF, - Te = Te, - react = 'ci', + mr = mr, + react = react, verbose = verbose, ) - # Formats output - FAC['rates']['ioniz'] = {} - FAC['rates']['ioniz']['ion_transitions'] = np.vstack( - (np.asarray(state), np.asarray(ion)) - ).T # dim(ntrans,2), Z state -> Z+1 state - FAC['rates']['ioniz']['ion_excit'] = np.asarray(data) # dim(ntrans, nt), [reduced rate] - if verbose == 1: - FAC['rates']['ioniz']['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, nt), [cm3/s] - FAC['rates']['ioniz']['XS_cm2'] = np.asarray(XS) # dim(ntrans, nE), [cm2] - FAC['rates']['ioniz']['engy_eV'] = np.asarray(engy) # dim(ntrans, nE), [eV] - - # Output - return FAC - -############################################################ -# -# fac.MaxwellRate Data Files -# -############################################################ - -# Reads Maxwell-averaged collisional excitation data files -def _ce_mr( - FAC=None, - fil=None, - verbose = None, - ): - - # Useful constants - eV2K = 11604.5 - - # Reads data file - mr = utils._read_mr( - fil=fil, - react='ce' - ) - - # Saves temperature grid data - # 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 = 'ce', - verbose=verbose, - ) - - # Formats output - FAC['rates']['excit'] = {} - FAC['rates']['excit']['col_transitions'] = np.vstack( - (np.asarray(st_upr), np.asarray(st_lwr)) - ).T # dim(ntrans,2), (upr,lwr) states in ColRadPy indices - FAC['rates']['excit']['col_excit'] = np.asarray(data) # dim(ntrans, nt), [cm3/s] - if verbose == 1: - FAC['rates']['excit']['col_trans_unfill'] = FAC['rates']['excit']['col_transitions'].copy() - FAC['rates']['excit']['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, nt), [cm3/s] + # 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 - -# Reads Maxwell-averaged Radiative recombination data -def _rr_mr( - FAC=None, - fil=None, - ): - - # Reads data file - mr = utils._read_mr( - fil=fil, - react='rr' - ) - - # Converts data file to ColRadPy form - data, ion, state, ratec = utils._conv_mr2colradpy( - FAC = FAC, - mr = mr, - react = 'rr', - ) - - # Formats output - FAC['rates']['recomb'] = {} - FAC['rates']['recomb']['recomb_transitions'] = np.vstack( - (np.asarray(ion), np.asarray(state)) - ).T # dim(ntrans,2), Z+1 state -> Z state - FAC['rates']['recomb']['recomb_excit'] = np.asarray(data) # dim(ntrans, nt), [cm3/s] - - # Ouput - return FAC - -# Reads Maxwell-averaged collisional ionization data -def _ci_mr( - FAC=None, - fil=None, - verbose=None, - ): - - # Reads data file - mr = utils._read_mr( - fil=fil, - react='ci' - ) - - # Converts data file to ColRadPy form - data, ion, state, ratec = utils._conv_mr2colradpy( - FAC = FAC, - mr = mr, - react = 'ci', - verbose=verbose, - ) - - # Formats output - FAC['rates']['ioniz'] = {} - FAC['rates']['ioniz']['ion_transitions'] = np.vstack( - (np.asarray(state), np.asarray(ion)) - ).T # dim(ntrans,2), Z state -> Z+1 state - FAC['rates']['ioniz']['ion_excit'] = np.asarray(data) # dim(ntrans, nt), [reduced rate] - if verbose == 1: - FAC['rates']['ioniz']['ratec_cm3/s'] = np.asarray(ratec) # dim(ntrans, nt), [cm3/s] - - # Ouput - return FAC + return FAC \ No newline at end of file From a29735b2e23949eda69e3bbdfce8000b57962b16 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 25 Jan 2024 12:31:24 -0500 Subject: [PATCH 32/46] [#18] added high-energy formula for rad recomb cross-section --- colradpy/convolve_EEDF.py | 65 +++++++++++++++++++++++++++++++++----- colradpy/read_FAC.py | 8 ++++- colradpy/read_FAC_utils.py | 2 ++ 3 files changed, 66 insertions(+), 9 deletions(-) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index da9ce21..d048268 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -5,6 +5,9 @@ cjperks, Jan 18th, 2024 +TO DO: + 1) high-energy asymptotic cross-section for radiative + recombination assumes you only have one ionized state ''' @@ -31,6 +34,7 @@ def convolve_EEDF( 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'] @@ -113,6 +117,7 @@ def convolve_EEDF( limit = limit, w_upr = w_upr, w_lwr = w_lwr, + ion_L = ion_L, verbose = verbose, use_rel=use_rel, react=react, @@ -179,6 +184,7 @@ def _calc_ratec( 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, @@ -218,11 +224,15 @@ def _calc_ratec( )(np.log10(engyEEDF[:,tt])) # dim(ngrid,), [cm2] # Fill missing data right at threshold - indE = np.where( - (engyEEDF[:,tt] >= dE[nn]) - & (engyEEDF[:,tt] <= engy_tmp[0]) - )[0] - XS_tmp[indE] = XS[0,nn] + ## 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, so hopefully this is just a small error !!! + 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 if limit is not None: @@ -234,6 +244,7 @@ def _calc_ratec( limit = limit[nn,:], w_upr = w_upr[nn], w_lwr = w_lwr[nn], + ion_L = ion_L[0], react = react, ) @@ -269,6 +280,7 @@ def _get_limit( limit = None, w_upr = None, w_lwr = None, + ion_L = None, react = None, ): @@ -307,11 +319,48 @@ def _get_limit( XS_tmp[indE] = ( np.pi *omega / k02 - / (1 +2*w_upr) + / (1 +2*w_lwr) ) * a02 # [cm2] elif react == 'rr': - blah = 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 @@ -335,7 +384,7 @@ def _get_limit( XS_tmp[indE] = ( omega / k02 - / (1 +2*w_upr) + / (1 +2*w_lwr) ) * a02 # [cm2] # Output diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 90eb4fb..41bd388 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -39,7 +39,7 @@ def read_FAC( 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 + reacts = None, # if None -> looks for all file suffixes Te = None, # if not using MaxwellRate files, [eV], dim(ntemp,) verbose = 1, ): @@ -189,6 +189,7 @@ def _en( 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 @@ -255,6 +256,10 @@ def _en( 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'], @@ -284,6 +289,7 @@ def _en( 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 diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py index 1175111..b0602f3 100644 --- a/colradpy/read_FAC_utils.py +++ b/colradpy/read_FAC_utils.py @@ -73,6 +73,7 @@ def _conv_ascii2colradpy( 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 = FAC['atomic']['ion_L'] verbose=verbose, use_rel = True, react = react, @@ -93,6 +94,7 @@ def _conv_ascii2colradpy( Bethe = 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 = FAC['atomic']['ion_L'] verbose=verbose, use_rel = True, react = react, From 2edfbfd97b01d219c794f985a5f242af7dfaee40 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 25 Jan 2024 13:21:37 -0500 Subject: [PATCH 33/46] [#18] added dielectronic capture rate coefficients. NOTE: Ignored autoionization rates --- colradpy/convolve_EEDF.py | 13 +---- colradpy/read_FAC.py | 110 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 112 insertions(+), 11 deletions(-) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index d048268..dbb59ee 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -100,8 +100,6 @@ def convolve_EEDF( engyEEDF=engyEEDF, XS=XS, engyXS=engyXS, - m=m, - dE=dE, use_rel=use_rel, ) # [cm3/s], dim(ntemp,ntrans) @@ -136,8 +134,6 @@ def _calc_DC( engyEEDF = None, # [eV], dim(ngrid, ntemp) XS = None, # [eV*cm2], dim(ntrans,) engyXS = None, # [eV], dim(ntrans,) - m = None, - dE = None, # [eV], dim(ntrans,) use_rel=None, ): @@ -146,14 +142,9 @@ def _calc_DC( # Loop over transitions for nn in np.arange(len(XS)): - if m == 0: - engy_tmp = engyXS[nn] - elif m == 1: - engy_tmp = engyXS[nn] + dE[nn] - # Incident electron velocity vel = _get_vel( - E_inc = engyEEDF[:,tt], + E_inc = engyXS[nn], use_rel = use_rel, ) # [cm/s], dim(ngrid,) @@ -165,7 +156,7 @@ def _calc_DC( np.log10(EEDF[:,tt]), bounds_error = False, fill_value = (-1e5,-1e5) - )(np.log10(engy_tmp)) # [1/eV] + )(np.log10(engyXS[nn])) # [1/eV] # Calculates rate coefficient, [cm3/s] ratec[tt,nn] = XS[nn] *vel *EEDF_tmp diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 41bd388..ffca2b4 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -422,6 +422,116 @@ def _tr( # ############################################################ +# Reads dielectronic recombination strength data files +## NOTE: Skips autoionization rates right now!!! +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] + + # Loop over blocks + for blk in np.arange(len(ai[1])): + # Loop over transitions + for tran 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 + ) + + 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 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(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) + + return FAC + # Reads cross-section data files def _get_xs( EEDF=None, From a7ecf7fc532c38501e4f10bde1e4720b5b7034af Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 25 Jan 2024 13:28:33 -0500 Subject: [PATCH 34/46] [#18] debugged having multiple ionized states for rad recomb cross-section high-energy formula --- colradpy/convolve_EEDF.py | 6 +----- colradpy/read_FAC.py | 3 +-- colradpy/read_FAC_utils.py | 10 ++++++++-- 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index dbb59ee..45726e7 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -5,10 +5,6 @@ cjperks, Jan 18th, 2024 -TO DO: - 1) high-energy asymptotic cross-section for radiative - recombination assumes you only have one ionized state - ''' # Modules @@ -235,7 +231,7 @@ def _calc_ratec( limit = limit[nn,:], w_upr = w_upr[nn], w_lwr = w_lwr[nn], - ion_L = ion_L[0], + ion_L = ion_L, react = react, ) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index ffca2b4..7e5467e 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -12,8 +12,7 @@ TO DO: 1) Add non-Maxwellian convolution 2) Missing energy level data (zpla, zpla1) - 3) Dielectronic recombination/autoionization data - 4) Charge exchange data + 3) Charge exchange data ''' diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py index b0602f3..10dc7ce 100644 --- a/colradpy/read_FAC_utils.py +++ b/colradpy/read_FAC_utils.py @@ -57,6 +57,12 @@ def _conv_ascii2colradpy( 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: ( @@ -73,7 +79,7 @@ def _conv_ascii2colradpy( 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 = FAC['atomic']['ion_L'] + ion_L = ion_L, verbose=verbose, use_rel = True, react = react, @@ -94,7 +100,7 @@ def _conv_ascii2colradpy( Bethe = 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 = FAC['atomic']['ion_L'] + ion_L = ion_L, verbose=verbose, use_rel = True, react = react, From 86f97aef9630e60b77a7c5754f704c16d0cdf716 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Tue, 20 Feb 2024 18:04:09 -0500 Subject: [PATCH 35/46] [#18] misc debug --- colradpy/convolve_EEDF.py | 10 +++++++--- colradpy/read_FAC.py | 19 ++++++++++++++----- colradpy/read_FAC_utils.py | 10 ++++++---- 3 files changed, 27 insertions(+), 12 deletions(-) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index 45726e7..9824be4 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -213,7 +213,8 @@ def _calc_ratec( # 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, so hopefully this is just a small error !!! + ## 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]) @@ -221,7 +222,7 @@ def _calc_ratec( )[0] XS_tmp[indE] = XS[0,nn] - # Fill values with high-energy asymptotic behavior if available + # 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, @@ -284,7 +285,7 @@ def _get_limit( indE = np.where(engyEEDF> engy_tmp[-1])[0] # If unnecessary - if len(indE) == 0: + if (len(indE) == 0 and react != 'rr'): return XS_tmp if react == 'ce': @@ -310,6 +311,9 @@ def _get_limit( ) * 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) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index 7e5467e..d076936 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -62,7 +62,7 @@ def read_FAC( 'tr', # Einstein coefficients 'ce', # Collisional excitation 'rr', # Radiative recombination - #'ai', # Autoionization/dielectronic recombination + 'ai', # Autoionization/dielectronic recombination 'ci', # Collision ionization ] @@ -123,10 +123,12 @@ def read_FAC( # Autoionization/dielectronic recombination if 'ai' in reacts: FAC = _ai( + nele=nele, FAC=FAC, fil=fil, EEDF=EEDF, Te=Te, + verbose=verbose, ) # Charge exchange @@ -422,7 +424,6 @@ def _tr( ############################################################ # Reads dielectronic recombination strength data files -## NOTE: Skips autoionization rates right now!!! def _ai( nele=None, EEDF=None, @@ -444,6 +445,9 @@ def _ai( dE = [] # [eV] DC = [] # [eV *cm2] + # Loads FAC AI data file + ai = rfac.read_ai(fil+'a.ai') + # Loop over blocks for blk in np.arange(len(ai[1])): # Loop over transitions @@ -528,6 +532,11 @@ def _ai( # Output FAC['rates']['recomb']['recomb_transitions'] = all_trans # dim(ntrans,2) FAC['rates']['recomb']['recomb_excit'] = all_recomb # [cm3/s], dim(ntrans, ntemp) + 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) return FAC @@ -537,7 +546,7 @@ def _get_xs( Te=None, # [eV], dim(ntemp,) FAC=None, fil=None, - vebose = None, + verbose = None, reacts = None, ): @@ -625,8 +634,8 @@ def _get_xs( 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), [cm2] - FAC['rates'][lbl]['engy_eV'] = np.asarray(engy) # dim(ntrans, nE), [eV] + 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() diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py index 10dc7ce..0c9db7a 100644 --- a/colradpy/read_FAC_utils.py +++ b/colradpy/read_FAC_utils.py @@ -85,8 +85,8 @@ def _conv_ascii2colradpy( react = react, ) data.append(tmp_data[:,0]) - XS.append(tmp_XS[:,0]) - engy.append(tmp_engy[:,0]) + XS.append(tmp_XS[:,:,0]) + engy.append(tmp_engy) else: data.append( @@ -283,12 +283,12 @@ def _read_ascii( upr_lbl = 'upper_index' data_lbl = 'crosssection' elif react == 'rr': - data_fil = rfac.read_rr(fil+'a.ce') + 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_rr(fil+'a.ce') + data_fil = rfac.read_ci(fil+'a.ci') lwr_lbl = 'bound_index' upr_lbl = 'free_index' data_lbl = 'crosssection' @@ -314,6 +314,8 @@ def _read_ascii( # Stores transition energy, [eV] out[lwr][upr]['dE'] = data_fil[1][blk]['Delta E'][trn] + if react == 'ce': + out[lwr][upr]['dE'] *= 1e3 if react == 'ce': # Stores energy grid in terms of incident electron energy, [eV] From d1a289a66b7713d8e56194204c52ec50b167ff38 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 29 Feb 2024 17:57:22 -0500 Subject: [PATCH 36/46] [#18] debug converting excitation rates to upsilon form --- colradpy/read_FAC_utils.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py index 0c9db7a..4203ff8 100644 --- a/colradpy/read_FAC_utils.py +++ b/colradpy/read_FAC_utils.py @@ -15,6 +15,7 @@ import numpy as np import copy from colradpy.convolve_EEDF import convolve_EEDF +import scipy.constants as cnt ############################################################ # @@ -247,14 +248,24 @@ def _conv_rate2upsilon( 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 return data * ( np.sqrt(np.asarray(Te_eV)/13.6058) - /2.1716e-8 - *(1 + 2*FAC['atomic']['w'][ind_upr]) + /factor + *(1 + 2*FAC['atomic']['w'][ind_lwr]) * np.exp( abs( FAC['atomic']['energy'][ind_upr] From 980eb31b0b4553ba944f49268adb0ca7a4292698 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 7 Mar 2024 16:22:37 -0500 Subject: [PATCH 37/46] [#18] added autoionization and various debug --- colradpy/colradpy_class.py | 42 +++++++++++++- colradpy/convolve_EEDF.py | 11 ++++ colradpy/read_FAC.py | 110 ++++++++++++++++++++++++++++--------- 3 files changed, 133 insertions(+), 30 deletions(-) diff --git a/colradpy/colradpy_class.py b/colradpy/colradpy_class.py index 52ee0e7..6f9cf35 100644 --- a/colradpy/colradpy_class.py +++ b/colradpy/colradpy_class.py @@ -886,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']),:] = \ @@ -944,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): @@ -984,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']),:,:] = \ @@ -1056,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] + @@ -2356,7 +2390,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. @@ -2390,7 +2424,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 @@ -2405,7 +2440,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) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index 9824be4..3dc0bf5 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -81,9 +81,17 @@ def convolve_EEDF( # 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]) else: print('NON-MAXWELLIAN ELECTRONS NOT IMPLEMENTED YET') @@ -378,6 +386,9 @@ def _get_limit( / (1 +2*w_lwr) ) * a02 # [cm2] + # Error check + XS_tmp[XS_tmp<0] = 0 + # Output return XS_tmp diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index d076936..e30a7de 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -22,6 +22,7 @@ import numpy as np import copy import colradpy.read_FAC_utils as utils +from colradpy.convolve_EEDF import convolve_EEDF ############################################################ # @@ -342,32 +343,6 @@ def _tr( # Transition array from TRTable() trans_ColRadPy = np.vstack((upr_ColRadPy, lwr_ColRadPy)).T # dim(ntrans,2) - # 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( - trans_ColRadPy == np.asarray([ - FAC['lvl_indices']['ColRadPy'][ind], 1 - ]), - axis = 1 - ) - )[0][0] - - a_val[ind2] += a_2photon - ## 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 @@ -409,6 +384,32 @@ def _tr( 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) @@ -444,6 +445,7 @@ def _ai( fre_ColRadPy = [] dE = [] # [eV] DC = [] # [eV *cm2] + AI = [] # [1/s] # Loads FAC AI data file ai = rfac.read_ai(fil+'a.ai') @@ -451,7 +453,7 @@ def _ai( # Loop over blocks for blk in np.arange(len(ai[1])): # Loop over transitions - for tran in np.arange(len(ai[1][blk]['Delta E'])): + for trn in np.arange(len(ai[1][blk]['Delta E'])): # Stores data dE.append( ai[1][blk]['Delta E'][trn] @@ -461,6 +463,10 @@ def _ai( ai[1][blk]['DC strength'][trn]*1e-20 ) + AI.append( + ai[1][blk]['AI rate'][trn] + ) + bnd.append( ai[1][blk]['bound_index'][trn] ) @@ -491,6 +497,48 @@ def _ai( # 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() @@ -519,6 +567,7 @@ def _ai( # Error check if len(ind_rr) != 1 and len(ind_dc) != 1: + print('DC Error') print(all_trans[tt,:]) print('xxxx') @@ -532,11 +581,18 @@ def _ai( # 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 From be9b5ba607a400871b66cf0e9d8567733e76a697 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 7 Mar 2024 18:10:14 -0500 Subject: [PATCH 38/46] [#18] debug random nans --- colradpy/read_FAC_utils.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py index 4203ff8..cd38408 100644 --- a/colradpy/read_FAC_utils.py +++ b/colradpy/read_FAC_utils.py @@ -262,7 +262,7 @@ def _conv_rate2upsilon( *cnt.hbar/cnt.m_e *1e4 # [cm^2/s] ) # [cm^3/s], 2.1716e-8 - return data * ( + tmp = data * ( np.sqrt(np.asarray(Te_eV)/13.6058) /factor *(1 + 2*FAC['atomic']['w'][ind_lwr]) @@ -275,6 +275,9 @@ def _conv_rate2upsilon( ) ) + tmp[np.isnan(tmp)] = 0 + return tmp + ############################################################ # # Reading Data Files From 13d29930beeae532467157b6a6891df8a97de7a4 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Fri, 12 Apr 2024 13:09:21 -0400 Subject: [PATCH 39/46] [#18] added some convenient information in write_adf15 --- colradpy/write_adf15.py | 121 +++++++++++++++++++++++++++++++++------- 1 file changed, 100 insertions(+), 21 deletions(-) diff --git a/colradpy/write_adf15.py b/colradpy/write_adf15.py index dcbc270..3673954 100644 --- a/colradpy/write_adf15.py +++ b/colradpy/write_adf15.py @@ -1,6 +1,6 @@ import numpy as np import itertools - +from datetime import date # group subset of an array # this is so there are 8 entries per line @@ -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(15, ' ') + +'-' + + ( + 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() From 8796d43b35f0c9cd190dd04fb2db5554da971c18 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Sat, 13 Apr 2024 11:01:37 -0400 Subject: [PATCH 40/46] [#18] debug for sim Kr --- colradpy/read_FAC_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py index cd38408..8dea518 100644 --- a/colradpy/read_FAC_utils.py +++ b/colradpy/read_FAC_utils.py @@ -329,7 +329,7 @@ def _read_ascii( # Stores transition energy, [eV] out[lwr][upr]['dE'] = data_fil[1][blk]['Delta E'][trn] if react == 'ce': - out[lwr][upr]['dE'] *= 1e3 + 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] From 99bba99338e3b1b02d4c606f4a6cf8e19ed8b6b0 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Sat, 13 Apr 2024 20:28:34 -0400 Subject: [PATCH 41/46] [#18] fixed overflow error in ioniz --- colradpy/read_FAC_utils.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py index 8dea518..1359790 100644 --- a/colradpy/read_FAC_utils.py +++ b/colradpy/read_FAC_utils.py @@ -123,13 +123,14 @@ def _conv_ascii2colradpy( elif react == 'ci': # Converts ionization data to adf04 reduced form - data[-1] = _conv_rate2reduct( + 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 @@ -194,13 +195,14 @@ def _conv_mr2colradpy( elif react == 'ci': # Converts ionization data to adf04 reduced form - data[-1] = _conv_rate2reduct( + 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 From 329ae167af107607f28f5169804d91d457b9e59a Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Mon, 13 May 2024 16:25:38 -0400 Subject: [PATCH 42/46] [#18] init including excitation-autoionization to SCD rates --- colradpy/colradpy_class.py | 67 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) diff --git a/colradpy/colradpy_class.py b/colradpy/colradpy_class.py index 6f9cf35..329e305 100644 --- a/colradpy/colradpy_class.py +++ b/colradpy/colradpy_class.py @@ -1302,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 @@ -1433,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', @@ -1595,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 From f24f4f4562b34a64d2af2eb4f1ef5a4848976172 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 16 May 2024 12:45:24 -0400 Subject: [PATCH 43/46] [#18] minor debug --- colradpy/read_FAC.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index e30a7de..c283b45 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -448,7 +448,11 @@ def _ai( AI = [] # [1/s] # Loads FAC AI data file - ai = rfac.read_ai(fil+'a.ai') + 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])): From b27a5bc91e7eefbf9e65cfec9458dc3b2031964a Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 16 May 2024 14:25:53 -0400 Subject: [PATCH 44/46] [#18] more padding in adf15 transition table --- colradpy/write_adf15.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/colradpy/write_adf15.py b/colradpy/write_adf15.py index 3673954..e8b7783 100644 --- a/colradpy/write_adf15.py +++ b/colradpy/write_adf15.py @@ -5,7 +5,7 @@ # 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: @@ -127,7 +127,7 @@ def write_adf15(fil_name,pec_inds, wave, pecs, pec_lvls, element_sym, +'('+str(atomic['S'][upr])+')' +str(atomic['L'][upr]) +'('+str(atomic['w'][upr])+')' - ).rjust(15, ' ') + ).rjust(17, ' ') +'-' + ( str(lwr) @@ -208,10 +208,10 @@ def write_adf15(fil_name,pec_inds, wave, pecs, pec_lvls, element_sym, f.write('C-----------------------------------------------------------------------\n') f.write('C\n') f.write( - "C ISEL WVLEN [A] TRANSITION lvl(S)L(w) TYPE\n" + "C ISEL WVLEN [A] TRANSITION lvl(S)L(w) TYPE\n" ) f.write( - "C ----- ------------ ------------------------------ -----\n" + "C ----- ------------ -------------------------------- -----\n" ) for tt in np.arange(len(trans)): f.write(trans[tt]) From f4cc8543fe45816fd3918b7a86823f0c4e536346 Mon Sep 17 00:00:00 2001 From: cjperks7 Date: Thu, 27 Jun 2024 11:10:39 -0400 Subject: [PATCH 45/46] [#18] added convolving Gaussian electron energy distribution functions to model beams such as EBIT plasmas --- colradpy/convolve_EEDF.py | 44 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 43 insertions(+), 1 deletion(-) diff --git a/colradpy/convolve_EEDF.py b/colradpy/convolve_EEDF.py index 3dc0bf5..a2d4d9f 100644 --- a/colradpy/convolve_EEDF.py +++ b/colradpy/convolve_EEDF.py @@ -43,7 +43,11 @@ def convolve_EEDF( if Te <10 keV = Maxwell-Boltzmann if Te >10 keV = Maxwell-Juttner (relativistic) - 2) NOT IMPLEMENTED YET (non-Maxwellians) + 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 @@ -93,6 +97,14 @@ def convolve_EEDF( 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) @@ -392,6 +404,36 @@ def _get_limit( # 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,) From 6e2f40ec0545b10ec22f044e5b867a227cd6ed30 Mon Sep 17 00:00:00 2001 From: cjperks Date: Thu, 27 Jun 2024 15:36:10 -0400 Subject: [PATCH 46/46] [#18] various minor debug --- colradpy/read_FAC.py | 13 +++++-------- colradpy/read_FAC_utils.py | 2 +- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/colradpy/read_FAC.py b/colradpy/read_FAC.py index c283b45..1e68bbd 100644 --- a/colradpy/read_FAC.py +++ b/colradpy/read_FAC.py @@ -18,7 +18,7 @@ # Module from pfac import rfac, crm -import os +import sys, os import numpy as np import copy import colradpy.read_FAC_utils as utils @@ -41,7 +41,7 @@ def read_FAC( 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 = 1, + verbose = 0, ): ######## -------- Determines which files to search for -------- ######## @@ -50,10 +50,10 @@ def read_FAC( if EEDF == 'Maxwellian_mr': # Use Maxwell-averaged rates from pfac.fac.MaxwellRate use_mr = True - elif EEDF == 'Maxwellian': + elif EEDF in ['Maxwellian', 'Gaussian']: use_mr = False else: - print('NON-MAXWELLIAN ELECTRONS NOT IMPLEMENTED YET') + print('REQUESTED ELECTRON ENERGY DISTRIBUTION NOT IMPLEMENTED YET') sys.exit(1) # FAC data files to search for @@ -85,7 +85,6 @@ def read_FAC( # Error check else: print('NEED TO INCLUDE ENERGY LEVEL DATA IN MODELING!!!') - sys.exit(1) # Gets rate coefficient data # Includes: collisional excit, radiative recomb, collisional ioniz @@ -119,7 +118,6 @@ def read_FAC( # Error check else: print('NEED TO INCLUDE EINSTEIN COEFFICIENT DATA IN MODELING!!!') - sys.exit(1) # Autoionization/dielectronic recombination if 'ai' in reacts: @@ -374,6 +372,7 @@ def _tr( # Error check if len(inda) != 1 and len(indc) != 1: + print('TR error') print(all_trans[tt,:]) print('xxxx') @@ -613,7 +612,6 @@ def _get_xs( # Error check if 'ce' not in reacts: print('NEED TO INCLUDE COLLISIONAL EXCITATION DATA IN MODELING!!!') - sys.exit(1) # Fill blank if 'rr' not in reacts: @@ -717,7 +715,6 @@ def _get_mr( # Error check if 'ce' not in reacts: print('NEED TO INCLUDE COLLISIONAL EXCITATION DATA IN MODELING!!!') - sys.exit(1) # Fill blank if 'rr' not in reacts: diff --git a/colradpy/read_FAC_utils.py b/colradpy/read_FAC_utils.py index 1359790..bd236d6 100644 --- a/colradpy/read_FAC_utils.py +++ b/colradpy/read_FAC_utils.py @@ -98,7 +98,7 @@ def _conv_ascii2colradpy( engyXS = XSdata[lwr][upr]['engy'], m = 0, dE = np.asarray([XSdata[lwr][upr]['dE']]), - Bethe = XSdata[lwr][upr]['limit'][None,:], + 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,