Replace ioneq with ionization_fraction and add setter (#337)

* first pass at replacing ioneq

* test fixes -- ._ioneq not removed

* copy_arrays -> memmap

* revert ioneq in examples -- asdf files still have ioneq key

* fix property names

* add setter for ionization fraction

* add setter test; remove ionization_filename

* fix tests

* fix bug with generate_hash_table

* fix test list to found other ioneq file

* add one new test

* Update fiasco/ions.py

Co-authored-by: Will Barnes <will.t.barnes@gmail.com>

* Update fiasco/ions.py

Co-authored-by: Will Barnes <will.t.barnes@gmail.com>

* Update fiasco/ions.py

Co-authored-by: Will Barnes <will.t.barnes@gmail.com>

* Update fiasco/ions.py

Co-authored-by: Will Barnes <will.t.barnes@gmail.com>

* Update fiasco/ions.py

Co-authored-by: Will Barnes <will.t.barnes@gmail.com>

* undo change to add_property and update setter test

* undo change to util

* fix util test

* remove redundant ValueError; fix test accordingly

* update example for new aiapy version

---------

Co-authored-by: Will Barnes <will.t.barnes@gmail.com>

Author: jwreep

docs/quick_start.rst

@@ -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>

examples/idl_comparisons/continuum_comparison.py

@@ -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'])
@@ -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'])
@@ -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()

examples/idl_comparisons/goft_comparison.py

@@ -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']))

examples/idl_comparisons/ioneq_comparison.py

@@ -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',
@@ -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()

examples/user_guide/aia_response.py

@@ -12,6 +12,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
+from aiapy.calibrate.util import get_correction_table
 from aiapy.response import Channel
 from astropy.visualization import quantity_support
 from scipy.interpolate import interp1d
@@ -26,7 +27,9 @@
 # First, create the `~aiapy.response.Channel` object for the
 # 94 angstrom channel and compute the wavelength response.
 ch = Channel(94*u.angstrom)
-response = ch.wavelength_response() * ch.plate_scale
+correction_table = get_correction_table("jsoc")
+response = ch.wavelength_response(correction_table=correction_table)
+response *= ch.plate_scale
 
 #############################################################
 # Plot the wavelength response

examples/user_guide/plot_ioneq.py examples/user_guide/plot_ionization_fraction.py

@@ -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()
@@ -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()

fiasco/base.py

@@ -142,6 +142,11 @@ 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):
     """
@@ -162,4 +167,3 @@ def property_template(self):
     add_property(IonBase, filetype)
     add_property(IonBase, '/'.join(['dielectronic', filetype]))
 add_property(IonBase, 'ip')
-add_property(IonBase, 'ioneq')

fiasco/collections.py

@@ -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
@@ -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
@@ -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
@@ -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
@@ -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

fiasco/elements.py

@@ -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
@@ -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

fiasco/fiasco.py

@@ -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)
@@ -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,