Skip to content
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.

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

Replace ioneq with ionization_fraction and add setter #337

Open
wants to merge 21 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 12 commits
Commits
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
2 changes: 1 addition & 1 deletion docs/quick_start.rst
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ Now that you have your database, you can use your ion object that you created ab
<BLANKLINE>
HDF5 Database: ...chianti_dbase.h5
Using Datasets:
ioneq: chianti
ionization_fraction: chianti
abundance: sun_coronal_1992_feldman_ext
ip: chianti
<BLANKLINE>
Expand Down
6 changes: 3 additions & 3 deletions examples/idl_comparisons/continuum_comparison.py
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ def plot_idl_comparison(wavelength, temperature, result_fiasco, result_idl):
# continuum emission, i.e. that emission produced by
# thermal bremsstrahlung.
idl_result_freefree = read_idl_test_output('freefree_all_ions', '8.0.7')
ion_kwargs = {'abundance': idl_result_freefree['abundance'], 'ioneq_filename': idl_result_freefree['ioneq']}
ion_kwargs = {'abundance': idl_result_freefree['abundance'], 'ionization_fraction': idl_result_freefree['ioneq']}
all_ions = [fiasco.Ion(ion_name, idl_result_freefree['temperature'], **ion_kwargs) for ion_name in fiasco.list_ions()]
all_ions = fiasco.IonCollection(*all_ions)
free_free = all_ions.free_free(idl_result_freefree['wavelength'])
Expand All @@ -82,7 +82,7 @@ def plot_idl_comparison(wavelength, temperature, result_fiasco, result_idl):
# Next, let's compare the outputs for the free-bound
# continuum emission.
idl_result_freebound = read_idl_test_output('freebound_all_ions', '8.0.7')
ion_kwargs = {'abundance': idl_result_freebound['abundance'], 'ioneq_filename': idl_result_freebound['ioneq']}
ion_kwargs = {'abundance': idl_result_freebound['abundance'], 'ionization_fraction': idl_result_freebound['ioneq']}
all_ions = [fiasco.Ion(ion_name, idl_result_freebound['temperature'], **ion_kwargs) for ion_name in fiasco.list_ions()]
all_ions = fiasco.IonCollection(*all_ions)
free_bound = all_ions.free_bound(idl_result_freebound['wavelength'])
Expand All @@ -99,7 +99,7 @@ def plot_idl_comparison(wavelength, temperature, result_fiasco, result_idl):
# Finally, let's compare the outputs for the two-photon
# continuum emission.
idl_result_twophoton = read_idl_test_output('twophoton_all_ions', '8.0.7')
ion_kwargs = {'abundance': idl_result_twophoton['abundance'], 'ioneq_filename': idl_result_twophoton['ioneq']}
ion_kwargs = {'abundance': idl_result_twophoton['abundance'], 'ionization_fraction': idl_result_twophoton['ioneq']}
all_ions = [fiasco.Ion(ion_name, idl_result_twophoton['temperature'], **ion_kwargs) for ion_name in fiasco.list_ions()]
all_ions = fiasco.IonCollection(*all_ions)
two_photon = all_ions.two_photon(idl_result_twophoton['wavelength'], idl_result_twophoton['density']).squeeze()
Expand Down
2 changes: 1 addition & 1 deletion examples/idl_comparisons/goft_comparison.py
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ def plot_idl_comparison(x, y_idl, y_python, fig, n_rows, i_row, quantity_name, t
ion = fiasco.Ion((idl_result['Z'], idl_result['iz']),
idl_result['temperature'],
abundance_filename=idl_result['abundance'],
ioneq_filename=idl_result['ioneq'])
ionization_fraction=idl_result['ioneq'])
contribution_func = ion.contribution_function(idl_result['density'])
transitions = ion.transitions.wavelength[~ion.transitions.is_twophoton]
idx = np.argmin(np.abs(transitions - idl_result['wavelength']))
Expand Down
18 changes: 9 additions & 9 deletions examples/idl_comparisons/ioneq_comparison.py
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ def plot_idl_comparison(x, y_idl, y_python, fig, n_rows, i_row, quantity_name, t
# as differences due to the interpolation schemes may show large deviations between
# the two approaches. However, the ionization fraction in these regions does not
# meaningfully contribute to any other derived quantities.
ioneq_files = [
ionization_files = [
'ioneq_1_1',
'ioneq_6_1',
'ioneq_6_2',
Expand All @@ -81,19 +81,19 @@ def plot_idl_comparison(x, y_idl, y_python, fig, n_rows, i_row, quantity_name, t
'ioneq_26_20',
'ioneq_26_27',
]
fig = plt.figure(figsize=(9,3*len(ioneq_files)), layout='constrained')
for i, name in enumerate(ioneq_files):
fig = plt.figure(figsize=(9,3*len(ionization_files)), layout='constrained')
for i, name in enumerate(ionization_files):
idl_result = read_idl_test_output(name, '8.0.7')
ion = fiasco.Ion((idl_result['Z'], idl_result['iz']),
idl_result['temperature'],
ioneq_filename=idl_result['ioneq_filename'])
ionization_filename=idl_result['ioneq_filename'])
element = fiasco.Element(ion.atomic_symbol, ion.temperature)
ioneq = element.equilibrium_ionization
print(f'IDL code to produce ioneq result for {ion.ion_name_roman}:')
ionization_fraction = element.equilibrium_ionization
print(f'IDL code to produce ionization_fraction result for {ion.ion_name_roman}:')
print(Template(idl_result['idl_script']).render(**idl_result))
axes = plot_idl_comparison(ion.temperature, idl_result['ioneq'], ion.ioneq,
fig, len(ioneq_files), 3*i, f'{ion.ion_name_roman}')
axes[0].plot(element.temperature, ioneq[:, ion.charge_state],
axes = plot_idl_comparison(ion.temperature, idl_result['ioneq'], ion.ionization_fraction,
fig, len(ionization_files), 3*i, f'{ion.ion_name_roman}')
axes[0].plot(element.temperature, ionization_fraction[:, ion.charge_state],
label='fiasco (rates)', color='C1', ls='-')
axes[0].legend()
plt.show()
Original file line number Diff line number Diff line change
Expand Up @@ -30,24 +30,24 @@
# We can use this to plot the population fraction of each ion as a
# function of temperature.
for ion in el:
ioneq = el.equilibrium_ionization[:, ion.charge_state]
imax = np.argmax(ioneq)
plt.plot(el.temperature, ioneq)
plt.text(el.temperature[imax], ioneq[imax], ion.ionization_stage_roman,
ionization_fraction = el.equilibrium_ionization[:, ion.charge_state]
imax = np.argmax(ionization_fraction)
plt.plot(el.temperature, ionization_fraction)
plt.text(el.temperature[imax], ionization_fraction[imax], ion.ionization_stage_roman,
horizontalalignment='center')
plt.xscale('log')
plt.title(f'{el.atomic_symbol} Equilibrium Ionization')
plt.show()

################################################
# The CHIANTI database also includes tabulated ionization equilibria for
# all ions in the database. The ``ioneq`` attribute on each
# all ions in the database. The ``ionization_fraction`` attribute on each
# `~fiasco.Ion` object returns the tabulated population
# fractions interpolated onto the ``temperature`` array.
# Note that these population fractions returned by `~fiasco.Ion.ioneq` are
# Note that these population fractions returned by `~fiasco.Ion.ionization_fraction` are
# stored in the CHIANTI database and therefore are set to NaN
# for temperatures outside of the tabulated temperature data given in CHIANTI.
plt.plot(el.temperature, el[11].ioneq)
plt.plot(el.temperature, el[11].ionization_fraction)
plt.xscale('log')
plt.title(f'{el[11].ion_name_roman} Equilibrium Ionization')
plt.show()
Expand All @@ -58,7 +58,7 @@
# the tabulated results were calculated or due to artifacts from the
# interpolation.
plt.plot(el.temperature, el.equilibrium_ionization[:, el[11].charge_state], label='calculated')
plt.plot(el.temperature, el[11].ioneq, label='interpolated')
plt.plot(el.temperature, el[11].ionization_fraction, label='interpolated')
plt.xlim(5e5, 5e6)
plt.xscale('log')
plt.legend()
Expand Down
13 changes: 10 additions & 3 deletions fiasco/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -142,8 +142,13 @@ def _abund(self):
data_path = '/'.join([self.atomic_symbol.lower(), 'abundance'])
return DataIndexer.create_indexer(self.hdf5_dbase_root, data_path)

@property
def _ion_fraction(self):
data_path = '/'.join([self.atomic_symbol.lower(), self._ion_name, 'ioneq'])
return DataIndexer.create_indexer(self.hdf5_dbase_root, data_path)


def add_property(cls, filetype):
def add_property(cls, filetype, property_name=None):
"""
Dynamically add filetype properties to base data access class
"""
Expand All @@ -152,7 +157,10 @@ def property_template(self):
return DataIndexer.create_indexer(self.hdf5_dbase_root, data_path)

property_template.__doc__ = f'Data in {filetype} type file'
property_template.__name__ = f'_{"_".join(filetype.split("/"))}'
if property_name:
property_template.__name__ = f'_{"_".join(property_name.split("/"))}'
else:
property_template.__name__ = f'_{"_".join(filetype.split("/"))}'
jwreep marked this conversation as resolved.
Show resolved Hide resolved
setattr(cls, property_template.__name__, property(property_template))


Expand All @@ -162,4 +170,3 @@ def property_template(self):
add_property(IonBase, filetype)
add_property(IonBase, '/'.join(['dielectronic', filetype]))
add_property(IonBase, 'ip')
add_property(IonBase, 'ioneq')
18 changes: 9 additions & 9 deletions fiasco/collections.py
Original file line number Diff line number Diff line change
Expand Up @@ -118,12 +118,12 @@ def free_free(self, wavelength: u.angstrom):
ff = ion.free_free(wavelength)
try:
abundance = ion.abundance
ioneq = ion.ioneq
ionization_fraction = ion.ionization_fraction
except MissingDatasetException as e:
self.log.warning(f'{ion.ion_name} not included in free-free emission. {e}')
continue
else:
free_free += ff * abundance * ioneq[:, np.newaxis]
free_free += ff * abundance * ionization_fraction[:, np.newaxis]
return free_free

@u.quantity_input
Expand Down Expand Up @@ -164,12 +164,12 @@ def free_bound(self, wavelength: u.angstrom, **kwargs):
# NOTE: the free-bound emissivity gets multiplied by the population
# fraction of the recombining ion, that is, the ion with one higher
# charge state.
ioneq = ion.next_ion().ioneq
ionization_fraction = ion.next_ion().ionization_fraction
except MissingDatasetException as e:
self.log.warning(f'{ion.ion_name} not included in free-bound emission. {e}')
continue
else:
free_bound += fb * abundance * ioneq[:, np.newaxis]
free_bound += fb * abundance * ionization_fraction[:, np.newaxis]
return free_bound

@u.quantity_input
Expand Down Expand Up @@ -211,7 +211,7 @@ def two_photon(self, wavelength: u.angstrom, electron_density: u.cm**-3, **kwarg
self.log.warning(f'{ion.ion_name} not included in two-photon emission. {e}')
continue
else:
two_photon += tp * ion.abundance * ion.ioneq[:, np.newaxis, np.newaxis]
two_photon += tp * ion.abundance * ion.ionization_fraction[:, np.newaxis, np.newaxis]
return two_photon

@u.quantity_input
Expand Down Expand Up @@ -399,12 +399,12 @@ def free_free_radiative_loss(self, use_itoh=False) -> u.Unit('erg cm3 s-1'):
try:
ff = ion.free_free_radiative_loss(use_itoh=use_itoh)
abundance = ion.abundance
ioneq = ion.ioneq
ionization_fraction = ion.ionization_fraction
except MissingDatasetException as e:
self.log.warning(f'{ion.ion_name} not included in free-free radiative loss. {e}')
continue
else:
free_free += ff * abundance * ioneq
free_free += ff * abundance * ionization_fraction
return free_free

@u.quantity_input
Expand All @@ -423,10 +423,10 @@ def free_bound_radiative_loss(self) -> u.Unit('erg cm3 s-1'):
try:
fb = ion.free_bound_radiative_loss()
abundance = ion.abundance
ioneq = ion.ioneq
ionization_fraction = ion.ionization_fraction
except MissingDatasetException as e:
self.log.warning(f'{ion.ion_name} not included in free-bound radiative loss. {e}')
continue
else:
free_bound += fb * abundance * ioneq
free_bound += fb * abundance * ionization_fraction
return free_bound
10 changes: 5 additions & 5 deletions fiasco/elements.py
Original file line number Diff line number Diff line change
Expand Up @@ -106,8 +106,8 @@ def equilibrium_ionization(self):
--------
>>> temperature = 10**np.arange(3.9, 6.5, 0.01) * u.K
>>> carbon = Element('C', temperature)
>>> carbon_ioneq = carbon.equilibrium_ionization
>>> carbon_ioneq[:, 4].max() # max populuation fraction of C V as a function of temperature
>>> carbon_ionization = carbon.equilibrium_ionization
>>> carbon_ionization[:, 4].max() # max population fraction of C V as a function of temperature
<Quantity 0.99776769>

See Also
Expand All @@ -120,10 +120,10 @@ def equilibrium_ionization(self):
# Select columns of V with smallest eigenvalues (returned in descending order)
# NOTE: must take the absolute value as the SVD solution is only accurate up
# to the sign. We require that the solutions must be positive.
ioneq = np.fabs(V[:, -1, :])
ioneq /= ioneq.sum(axis=1)[:, np.newaxis]
ionization_fraction = np.fabs(V[:, -1, :])
ionization_fraction /= ionization_fraction.sum(axis=1)[:, np.newaxis]

return u.Quantity(ioneq)
return u.Quantity(ionization_fraction)

def __getitem__(self, value):
if isinstance(value, (str, tuple)): # NOQA: UP038
Expand Down
21 changes: 11 additions & 10 deletions fiasco/fiasco.py
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,7 @@ def proton_electron_ratio(temperature: u.K, **kwargs):
# Import here to avoid circular imports
from fiasco import log
h_2 = fiasco.Ion('H +1', temperature, **kwargs)
numerator = h_2.abundance * h_2._ioneq[h_2._instance_kwargs['ioneq_filename']]['ionization_fraction']
numerator = h_2.abundance * h_2._ion_fraction[h_2._instance_kwargs['ionization_fraction']]['ionization_fraction']
denominator = u.Quantity(np.zeros(numerator.shape))
for el_name in list_elements(h_2.hdf5_dbase_root):
el = fiasco.Element(el_name, temperature, **h_2._instance_kwargs)
Expand All @@ -138,21 +138,22 @@ def proton_electron_ratio(temperature: u.K, **kwargs):
f'Not including {el.atomic_symbol}. Abundance not available from {abund_file}.')
continue
for ion in el:
ioneq_file = ion._instance_kwargs['ioneq_filename']
# NOTE: We use ._ioneq here rather than .ioneq to avoid doing an interpolation to the
# temperature array every single time and instead only interpolate once at the end.
# It is assumed that the ioneq temperature array for each ion is the same.
ionization_file = ion._instance_kwargs['ionization_fraction']
# NOTE: We use ._ion_fraction here rather than .ionization_fraction to avoid
# doing an interpolation to the temperature array every single time and instead only
# interpolate once at the end.
# It is assumed that the ionization_fraction temperature array for each ion is the same.
try:
ioneq = ion._ioneq[ioneq_file]['ionization_fraction']
t_ioneq = ion._ioneq[ioneq_file]['temperature']
ionization_fraction = ion._ion_fraction[ionization_file]['ionization_fraction']
t_ionization_fraction = ion._ion_fraction[ionization_file]['temperature']
except KeyError:
log.warning(
f'Not including {ion.ion_name}. Ionization fraction not available from {ioneq_file}.')
f'Not including {ion.ion_name}. Ionization fraction not available from {ionization_file}.')
continue
denominator += ioneq * abundance * ion.charge_state
denominator += ionization_fraction * abundance * ion.charge_state

ratio = numerator / denominator
f_interp = interp1d(t_ioneq.to(temperature.unit).value,
f_interp = interp1d(t_ionization_fraction.to(temperature.unit).value,
ratio.value,
kind='linear',
bounds_error=False,
Expand Down
Loading
Loading