apply 0.83 scaling factor to fiasco goft, use iron (#288)

docs/conf.py

@@ -15,7 +15,7 @@
 
 # -- Project information -----------------------------------------------------
 project = 'fiasco'
-copyright = f'{datetime.datetime.now(datetime.UTC).year}, Will Barnes'
+copyright = f'{datetime.datetime.utcnow().year}, Will Barnes'
 author = 'Will Barnes'
 # The full version, including alpha/beta/rc tags
 from fiasco import __version__

examples/idl_comparisons/goft_comparison.py

@@ -85,7 +85,8 @@ def plot_idl_comparison(x, y_idl, y_python, fig, n_rows, i_row, quantity_name, t
     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']))
-    goft = contribution_func[:, 0, idx]
+    # NOTE: Multiply by 0.83 because the fiasco calculation does not include the n_H/n_e ratio
+    goft = contribution_func[:, 0, idx] * 0.83
     axes = plot_idl_comparison(ion.temperature, idl_result['contribution_function'], goft,
                                fig, len(goft_files), 3*i, f'{ion.ion_name_roman}')
     axes[0].legend()

examples/user_guide/nei_populations.py

@@ -39,14 +39,14 @@
 # to :math:`10^4` K.
 t = np.arange(0, 30, 0.01) * u.s
 Te_min = 1e4 * u.K
-Te_max = 2e6 * u.K
+Te_max = 1.5e7 * u.K
 t_mid = 15 * u.s
 Te = Te_min + (Te_max - Te_min) / t_mid * t
 Te[t>t_mid] = Te_max - (Te_max - Te_min) / t_mid * (t[t>t_mid] - t_mid)
 
 ###############################################################
 # Similarly, we'll model the density as increasing from
-# :math:`10^8\,\mathrm{cm}^{-3}` to :math:`1.5^10\,\mathrm{cm}^{-3}`
+# :math:`10^8\,\mathrm{cm}^{-3}` to :math:`10^{10}\,\mathrm{cm}^{-3}`
 # and then back down to :math:`10^8\,\mathrm{cm}^{-3}`.
 ne_min = 1e8 * u.cm**(-3)
 ne_max = 1e10 * u.cm**(-3)
@@ -69,19 +69,19 @@
 # As shown in previous examples, we can get the ion population
 # fractions of any element in the CHIANTI database over some
 # temperature array with `~fiasco.Element.equilibrium_ionization`.
-# First, let's model the ionization fractions of C for the above
+# First, let's model the ionization fractions of Fe for the above
 # temperature model, assuming ionization equilibrium.
 temperature_array = np.logspace(4, 8, 1000) * u.K
-carbon = Element('carbon', temperature_array)
-func_interp = interp1d(carbon.temperature.to_value('K'), carbon.equilibrium_ionization.value,
+iron = Element('iron', temperature_array)
+func_interp = interp1d(iron.temperature.to_value('K'), iron.equilibrium_ionization.value,
                        axis=0, kind='cubic', fill_value='extrapolate')
-carbon_ieq = u.Quantity(func_interp(Te.to_value('K')))
+iron_ieq = u.Quantity(func_interp(Te.to_value('K')))
 
 ###############################################################
 # We can plot the population fractions as a function of time.
 plt.figure(figsize=(12, 4))
-for ion in carbon:
-    plt.plot(t, carbon_ieq[:, ion.charge_state], label=ion.ion_name_roman)
+for ion in iron:
+    plt.plot(t, iron_ieq[:, ion.charge_state], label=ion.ion_name_roman)
 plt.xlim(t[[0,-1]].value)
 plt.legend(ncol=4, frameon=False)
 plt.show()
@@ -103,47 +103,48 @@
 # is the time step.
 #
 # First, we can set our initial state to the equilibrium populations.
-carbon_nei = np.zeros(t.shape + (carbon.atomic_number + 1,))
-carbon_nei[0, :] = carbon_ieq[0,:]
+iron_nei = np.zeros(t.shape + (iron.atomic_number + 1,))
+iron_nei[0, :] = iron_ieq[0,:]
 
 ###############################################################
 # Then, interpolate the rate matrix to our modeled temperatures
-func_interp = interp1d(carbon.temperature.to_value('K'), carbon._rate_matrix.value,
+func_interp = interp1d(iron.temperature.to_value('K'), iron._rate_matrix.value,
                        axis=0, kind='cubic', fill_value='extrapolate')
-fe_rate_matrix = func_interp(Te.to_value('K')) * carbon._rate_matrix.unit
+fe_rate_matrix = func_interp(Te.to_value('K')) * iron._rate_matrix.unit
 
 ###############################################################
 # Finally, we can iteratively compute the non-equilibrium ion
 # fractions using the above equation.
-identity = u.Quantity(np.eye(carbon.atomic_number + 1))
+identity = u.Quantity(np.eye(iron.atomic_number + 1))
 for i in range(1, t.shape[0]):
     dt = t[i] - t[i-1]
     term1 = identity - ne[i] * dt/2. * fe_rate_matrix[i, ...]
     term2 = identity + ne[i-1] * dt/2. * fe_rate_matrix[i-1, ...]
-    carbon_nei[i, :] = np.linalg.inv(term1) @ term2 @ carbon_nei[i-1, :]
-    carbon_nei[i, :] = np.fabs(carbon_nei[i, :])
-    carbon_nei[i, :] /= carbon_nei[i, :].sum()
+    iron_nei[i, :] = np.linalg.inv(term1) @ term2 @ iron_nei[i-1, :]
+    iron_nei[i, :] = np.fabs(iron_nei[i, :])
+    iron_nei[i, :] /= iron_nei[i, :].sum()
 
-carbon_nei = u.Quantity(carbon_nei)
+iron_nei = u.Quantity(iron_nei)
 
 ###############################################################
 # And plot the non-equilibrium populations as a function of time
 plt.figure(figsize=(12,4))
-for ion in carbon:
-    plt.plot(t, carbon_nei[:, ion.charge_state], ls='-', label=ion.ion_name_roman,)
+for ion in iron:
+    plt.plot(t, iron_nei[:, ion.charge_state], ls='-', label=ion.ion_name_roman,)
 plt.xlim(t[[0,-1]].value)
 plt.legend(ncol=4, frameon=False)
 plt.show()
 
 ###############################################################
 # We can compare the two equilibrium and non equilbrium cases
-# directly by overplotting the C V population fractions for
+# directly by overplotting the Fe XVIII population fractions for
 # the two cases.
 plt.figure(figsize=(12,4))
-plt.plot(t, carbon_ieq[:, 4], label='IEQ')
-plt.plot(t, carbon_nei[:, 4], label='NEI',)
+plt.plot(t, iron_ieq[:, 17], label='IEQ')
+plt.plot(t, iron_nei[:, 17], label='NEI',)
 plt.xlim(t[[0,-1]].value)
 plt.legend(frameon=False)
+plt.title(f'IEQ versus NEI for {iron[17].ion_name_roman}')
 plt.show()
 
 ###############################################################
@@ -156,8 +157,8 @@
 # temperature that one would infer if they observed the
 # non-equilibrium populations and assumed the populations
 # were in equilibrium.
-Te_eff =  carbon.temperature[[(np.fabs(carbon.equilibrium_ionization - carbon_nei[i, :])).sum(axis=1).argmin()
-                          for i in range(carbon_nei.shape[0])]]
+Te_eff =  iron.temperature[[(np.fabs(iron.equilibrium_ionization - iron_nei[i, :])).sum(axis=1).argmin()
+                          for i in range(iron_nei.shape[0])]]
 plt.plot(t, Te.to('MK'), label=r'$T$')
 plt.plot(t, Te_eff.to('MK'), label=r'$T_{\mathrm{eff}}$')
 plt.xlim(t[[0,-1]].value)

examples/user_guide/plot_ioneq.py

@@ -18,8 +18,8 @@
 
 ################################################
 # First, create the `~fiasco.Element` object for carbon.
-temperature = 10**np.arange(3.9, 6.5, 0.01) * u.K
-el = Element('C', temperature)
+temperature = 10**np.arange(4, 8, 0.01) * u.K
+el = Element('Fe', temperature)
 
 ################################################
 # The ionization fractions in equilibrium can be determined by calculating the
@@ -47,20 +47,20 @@
 # Note that these population fractions returned by `~fiasco.Ion.ioneq` 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[3].ioneq)
+plt.plot(el.temperature, el[11].ioneq)
 plt.xscale('log')
-plt.title(f'{el[3].ion_name_roman} Equilibrium Ionization')
+plt.title(f'{el[11].ion_name_roman} Equilibrium Ionization')
 plt.show()
 
 ################################################
 # We can then compare tabulated and calculated results for a single ion.
 # Note that the two may not be equal due to differences in the rates when
 # the tabulated results were calculated or due to artifacts from the
 # interpolation.
-plt.plot(el.temperature, el.equilibrium_ionization[:, el[3].charge_state], label='calculated')
-plt.plot(el.temperature, el[3].ioneq, label='interpolated')
-plt.xlim(4e4, 3e5)
+plt.plot(el.temperature, el.equilibrium_ionization[:, el[11].charge_state], label='calculated')
+plt.plot(el.temperature, el[11].ioneq, label='interpolated')
+plt.xlim(5e5, 5e6)
 plt.xscale('log')
 plt.legend()
-plt.title(f'{el[3].ion_name_roman} Equilibrium Ionization')
+plt.title(f'{el[11].ion_name_roman} Equilibrium Ionization')
 plt.show()