Skip to content

Ionization balance with charge exchange using ADF11 rates #27

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dmkriete wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Ionization balance with charge exchange using ADF11 rates #27

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b6dff9e
Can now pull charge exchange rate coefficients from ADF11 files for p…
Jul 2, 2024
2e8643d
Energy balance can now include CX in the calculation of the underlyin…
Jul 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 29 additions & 3 deletions colradpy/energy_balance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,11 @@ class energy_balance(object):
densities/temperatures at a time (it is not vectorized across a grid of
temperatures/densities like the ionization balance).

Currently, if charge exchange effects are turned on, they only directly
affect the ionization balance. The energy balance is then indirectly
affected through the ionization balance. The direct effects of charge
exchange on the energy balance are not modelled.

:Args:
:param fils: Array of the input adf04 or adf11 files
:type fils: string array
Expand Down Expand Up @@ -54,7 +59,16 @@ class energy_balance(object):

:param ion_charge_number: Main ion species charge number
:type ion_charge_number: int


:param use_cx: Include thermal charge exchange in the ionization balance
:param use_cx: bool

:param hydrogen_temp: Neutral hydrogen temperature in (eV)
:type hydrogen_temp: float

:param hydrogen_dens: Neutral hydrogen density in (cm^-3)
:type hydrogen_dens: float

:param init_temp: Initial temperatures in (eV) of the charge states for the modelled species
:type init_temp: float array

Expand All @@ -63,12 +77,19 @@ class energy_balance(object):

def __init__(self, fils, metas, mass, polarizability, electron_temp,
electron_dens, ion_temp, ion_dens, ion_mass, ion_charge_number,
use_cx=False, hydrogen_temp=None, hydrogen_dens=None,
init_temp=np.array([]), **kwargs):

# Set up ionization balance
self.ion_balance = ionization_balance(
fils, metas, temp_grid=np.array([electron_temp]),
dens_grid=np.array([electron_dens]), **kwargs,
fils,
metas,
temp_grid=np.array([electron_temp]),
dens_grid=np.array([electron_dens]),
htemp_grid=np.array([hydrogen_temp]),
hdens_grid=np.array([hydrogen_dens]),
use_cx=use_cx,
**kwargs,
)

# Initialize data
Expand All @@ -83,6 +104,10 @@ def __init__(self, fils, metas, mass, polarizability, electron_temp,
self.data["user"]["ion_mass"] = ion_mass
self.data["user"]["ion_charge_number"] = ion_charge_number
self.data["user"]["init_temp"] = init_temp
self.data["user"]["use_cx"] = use_cx
if use_cx:
self.data["user"]["hydrogen_temp"] = hydrogen_temp
self.data["user"]["hydrogen_dens"] = hydrogen_dens


def solve(self, n0=np.array([]), T0=np.array([]), td_t=np.array([])):
Expand Down Expand Up @@ -169,6 +194,7 @@ def _evolve_energy(self, time, energies, pops_fun):

Returns a vector of the time derivative of the energy in each charge state
"""
# TODO: Implement effects of charge exchange reactions on energy balance
# Get fractional abundance for this time and calculate state temperatures
pops = pops_fun(np.array([time]))[:, 0, 0, 0]
# Handle special cases where populations are zero (typically at beginning)
Expand Down
28 changes: 16 additions & 12 deletions colradpy/ionization_balance_class.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,16 +35,18 @@ class ionization_balance():
:param metas: List of arrays for the metastable levels in a charge state
:type metas: list

:param temp_grid: Array of the temperatures in (eV)
:param temp_grid: Array of electron temperatures in (eV)
:type temp_grid: float array

:param dens_grid: Array of the densities in (cm-3)
:param dens_grid: Array of electron densities in (cm-3)
:type dens_grid: float array

:param htemp_grid: Temperature grid of thermal CX hydrogen (eV)
:param htemp_grid: Temperature grid of thermal CX hydrogen (eV). Must be
aligned with the electron temperature grid.
:type htemp_grid: float array

:param hdens_grid: Density grid of the thermal CX hydrogen densities in (cm-3)
:param hdens_grid: Density grid of the thermal CX hydrogen densities in
(cm^-3). Must be aligned with the electron density grid.
:type hdens_grid: float array

:param soln_times: Times to calculate the solution for the time dependent solutions (s)
Expand Down Expand Up @@ -120,7 +122,9 @@ def __init__(self,fils,metas,temp_grid, dens_grid, htemp_grid=np.array([]), hden
self.data['cr_data']['gcrs'] = {}
self.data['user'] = {}
self.data['user']['temp_grid'] = np.asarray(temp_grid) #eV
self.data['user']['dens_grid'] = np.asarray(dens_grid)#cm-3
self.data['user']['dens_grid'] = np.asarray(dens_grid) #cm^-3
self.data['user']['htemp_grid'] = np.asarray(htemp_grid) #eV
self.data['user']['hdens_grid'] = np.asarray(hdens_grid) # cm^-3
self.data['user']['fils'] = np.asarray(fils)
self.data['user']['init_abund'] = np.asarray(init_abund)
self.data['user']['soln_times'] = np.asarray(soln_times)
Expand All @@ -147,9 +151,9 @@ def __init__(self,fils,metas,temp_grid, dens_grid, htemp_grid=np.array([]), hden


if(type(use_cx) == bool):
self.data['user']['use_cx'] = np.ones(len(adf11['input_file']['metas'])-1,dtype=bool)
self.data['user']['use_cx'][:] = False
self.data['user']['use_cx'] = np.ones(len(adf11['input_file']['metas'])-1, dtype=bool)
self.data['user']['use_cx'][:] = use_cx


for j in range(0,len(adf11['input_file']['metas'])-1):
if( str(j) not in self.data['cr_data']['gcrs']):
Expand All @@ -168,20 +172,20 @@ def __init__(self,fils,metas,temp_grid, dens_grid, htemp_grid=np.array([]), hden
self.data['input_file']['acd'] = adf11['input_file']

if( 'qcd' in fils[i]):
self.data['cr_data']['gcrs'][str(j)]['qcd']= interp_rates_adf11(adf11['input_file']['temp_grid'],
self.data['cr_data']['gcrs'][str(j)]['qcd'] = interp_rates_adf11(adf11['input_file']['temp_grid'],
adf11['input_file']['dens_grid'],
temp_grid,dens_grid,adf11['input_file'][str(j)])
self.data['input_file']['qcd'] = adf11['input_file']

if( 'xcd' in fils[i]):
self.data['cr_data']['gcrs'][str(j)]['xcd']= interp_rates_adf11(adf11['input_file']['temp_grid'],
self.data['cr_data']['gcrs'][str(j)]['xcd'] = interp_rates_adf11(adf11['input_file']['temp_grid'],
adf11['input_file']['dens_grid'],
temp_grid,dens_grid,adf11['input_file'][str(j)])
self.data['input_file']['xcd'] = adf11['input_file']
if( 'ccd' in fils[i]):
self.data['cr_data']['gcrs'][str(j)]['ccd']= interp_rates_adf11(adf11['input_file']['temp_grid'],
self.data['cr_data']['gcrs'][str(j)]['ccd'] = interp_rates_adf11(adf11['input_file']['temp_grid'],
adf11['input_file']['dens_grid'],
temp_grid,dens_grid,adf11['input_file'][str(j)])
htemp_grid,hdens_grid,adf11['input_file'][str(j)])
self.data['input_file']['ccd'] = adf11['input_file']

if('qcd' not in self.data['cr_data']['gcrs'][str(j)]):
Expand Down
16 changes: 8 additions & 8 deletions colradpy/read_adf11.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,8 @@ def read_adf11(fil):
adf11['input_file']['charge_max'] = int(tmp[4])

f.readline() #reading '-------------'
if('r' in re.split('_', os.path.split(fil)[-1])[0]):
file_type = re.split('_', os.path.split(fil)[-1])[0]
if('r' in file_type):
adf11['input_file']['metas'] = np.array(list( #metastables
map(int,re.findall('(\d+)',f.readline()))))
f.readline() #reading '---------------'
Expand All @@ -49,32 +49,32 @@ def read_adf11(fil):
num_stages = len(adf11['input_file']['metas'])
for i in range(0,num_stages):
adf11['input_file'][str(i)] = {}
if( 'scd' in fil or 'acd' in fil):
if any(x in file_type for x in ['scd', 'acd', 'ccd']):
adf11['input_file'][str(i)] = np.zeros((adf11['input_file']['metas'][i],
adf11['input_file']['metas'][i+1],
len(adf11['input_file']['temp_grid']),
len(adf11['input_file']['dens_grid'])))
if( 'qcd' in fil ):
if( 'qcd' in file_type ):
adf11['input_file'][str(i)] = np.zeros((adf11['input_file']['metas'][i],
adf11['input_file']['metas'][i],
len(adf11['input_file']['temp_grid']),
len(adf11['input_file']['dens_grid'])))
if('xcd' in fil):
if('xcd' in file_type):
adf11['input_file'][str(i)] = np.zeros((adf11['input_file']['metas'][i+1],
adf11['input_file']['metas'][i+1],
len(adf11['input_file']['temp_grid']),
len(adf11['input_file']['dens_grid'])))
if('plt' in fil):
if('plt' in file_type):
adf11['input_file'][str(i)] = np.zeros((adf11['input_file']['metas'][i],
adf11['input_file']['metas'][i],
len(adf11['input_file']['temp_grid']),
len(adf11['input_file']['dens_grid'])))
if('prb' in fil):
if('prb' in file_type):
adf11['input_file'][str(i)] = np.zeros((adf11['input_file']['metas'][i],
adf11['input_file']['metas'][i],
len(adf11['input_file']['temp_grid']),
len(adf11['input_file']['dens_grid'])))

#Reading the GCR value portion
gcr_line = f.readline()
ii = 0
Expand Down