.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/gallery/idl_comparisons/goft_comparison.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_generated_gallery_idl_comparisons_goft_comparison.py: CHIANTI IDL Comparison: Contribution Functions =============================================== Compare the contribution function calculation to that in the CHIANTI IDL routines for a few select ions. .. GENERATED FROM PYTHON SOURCE LINES 7-20 .. code-block:: Python import hissw import matplotlib.pyplot as plt import numpy as np from astropy.visualization import quantity_support import fiasco from fiasco.tests.idl.helpers import read_idl_test_output quantity_support() .. rst-class:: sphx-glr-script-out .. code-block:: none .MplQuantityConverter object at 0x7fc3a25256d0> .. GENERATED FROM PYTHON SOURCE LINES 21-24 Define a function for comparing and plotting the outputs from fiasco and IDL. Note that we have precomputed the IDL outputs. .. GENERATED FROM PYTHON SOURCE LINES 24-63 .. code-block:: Python def plot_idl_comparison(x, y_idl, y_python, fig, n_rows, i_row, quantity_name, thresh=1e-6): # Highlight region where value of goft is within some value of the peak idx = np.where(y_idl>=y_idl.max()*thresh) x_thresh_min, x_thresh_max = x[idx[0][[0,-1]]] # Direct comparison ax1 = fig.add_subplot(n_rows, 3, i_row+1) ax1.plot(x, y_python, label='fiasco') ax1.plot(x, y_idl, color='k', marker='o', ls='--', markevery=2, label='IDL') ax1.axvspan(x[0], x_thresh_min, color='r', alpha=0.25) ax1.axvspan(x_thresh_max, x[-1], color='r', alpha=0.25) ax1.set_xscale('log') ax1.set_yscale('log') if i_row == 0: ax1.set_title('Ionization Fraction') ax1.set_xlim(1e5, 1e8) ax1.set_ylim(y_python.max()*np.array([thresh,2])) ax1.set_ylabel(quantity_name) # Ratio ax2 = fig.add_subplot(n_rows, 3, i_row+2, sharex=ax1) ax2.plot(x, (y_python/y_idl).decompose()) ax2.axvspan(x[0], x_thresh_min, color='r', alpha=0.25) ax2.axvspan(x_thresh_max, x[-1], color='r', alpha=0.25) ax2.axhline(y=1, color='k', ls=':') if i_row == 0: ax2.set_title('fiasco / IDL') diff_limit = .05 ax2.set_ylim(1-diff_limit, 1+diff_limit) # Normalized difference ax3 = fig.add_subplot(n_rows, 3, i_row+3, sharex=ax1) ax3.plot(x, ((y_python - y_idl)/y_idl).decompose()) ax3.axvspan(x[0], x_thresh_min, color='r', alpha=0.25) ax3.axvspan(x_thresh_max, x[-1], color='r', alpha=0.25) ax3.axhline(y=0, color='k', ls=':') if i_row == 0: ax3.set_title('(fiasco - IDL)/IDL') ax3.set_ylim(-diff_limit, diff_limit) return ax1, ax2, ax3 .. GENERATED FROM PYTHON SOURCE LINES 64-70 Next, plot the comparison between the ionization fraction results for a selected number of ions. The regions highlighted in red denote places where the contribution function is less than :math:`10^{-6}` times the peak of the ionization fraction. In these regions, the comparison is between the two results is not as critical as the contribution function is comparatively small in these regions. .. GENERATED FROM PYTHON SOURCE LINES 70-94 .. code-block:: Python goft_files = [ 'goft_20_15_200.972', 'goft_26_9_171.073', 'goft_26_11_188.497', 'goft_26_14_197.862', 'goft_26_16_262.984', ] fig = plt.figure(figsize=(9,3*len(goft_files)), layout='constrained') template_env = hissw.Environment(ssw_home='', idl_home='').env for i, name in enumerate(goft_files): idl_result = read_idl_test_output(name, '8.0.7') ion = fiasco.Ion((idl_result['Z'], idl_result['iz']), idl_result['temperature'], abundance_filename=idl_result['abundance'], ioneq_filename=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'])) goft = contribution_func[:, 0, idx] 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() print(f'IDL code to produce {ion.ion_name_roman} contribution function result:') print(template_env.from_string(idl_result['idl_script']).render(**idl_result)) .. image-sg:: /generated/gallery/idl_comparisons/images/sphx_glr_goft_comparison_001.png :alt: Ionization Fraction, fiasco / IDL, (fiasco - IDL)/IDL :srcset: /generated/gallery/idl_comparisons/images/sphx_glr_goft_comparison_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/astropy/units/format/utils.py:219: UnitsWarning: The unit 'erg' has been deprecated in the VOUnit standard. Suggested: cm**2.g.s**-2. warnings.warn(message, UnitsWarning) /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/astropy/units/format/utils.py:219: UnitsWarning: The unit 'Angstrom' has been deprecated in the VOUnit standard. Suggested: 0.1nm. warnings.warn(message, UnitsWarning) /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/asdf/_asdf.py:359: AsdfWarning: File 'file:///home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/fiasco/tests/idl/data/goft_20_15_200.972_v8.0.7.asdf' was created with extension URI 'asdf://astropy.org/core/extensions/core-1.5.0' (from package asdf-astropy==0.5.0), which is not currently installed warnings.warn(msg, AsdfWarning) /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/astropy/units/quantity.py:671: RuntimeWarning: invalid value encountered in divide result = super().__array_ufunc__(function, method, *arrays, **kwargs) IDL code to produce Ca XV contribution function result: abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR='abundance') ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq') density = 10d / 1d wave_min = 7035889255347913d / 35184372088832d wave_max = 7106257999525577d / 35184372088832d contribution_function = g_of_t(20,$ 15,$ dens=density,$ abund_file=abund_file,$ ioneq_file=ioneq_file,$ index=297,/quiet,$ wrange=[wave_min, wave_max]) ; Call this function to get the temperature array read_ioneq,ioneq_file,temperature,ioneq,ref /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/asdf/_asdf.py:359: AsdfWarning: File 'file:///home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/fiasco/tests/idl/data/goft_26_9_171.073_v8.0.7.asdf' was created with extension URI 'asdf://astropy.org/core/extensions/core-1.5.0' (from package asdf-astropy==0.5.0), which is not currently installed warnings.warn(msg, AsdfWarning) WARNING: No proton data available for Fe 9. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions] /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/astropy/units/quantity.py:671: RuntimeWarning: divide by zero encountered in divide result = super().__array_ufunc__(function, method, *arrays, **kwargs) IDL code to produce Fe IX contribution function result: abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR='abundance') ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq') density = 10d / 1d wave_min = 5983911714263925d / 35184372088832d wave_max = 6054280458441589d / 35184372088832d contribution_function = g_of_t(26,$ 9,$ dens=density,$ abund_file=abund_file,$ ioneq_file=ioneq_file,$ index=13946,/quiet,$ wrange=[wave_min, wave_max]) ; Call this function to get the temperature array read_ioneq,ioneq_file,temperature,ioneq,ref /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/asdf/_asdf.py:359: AsdfWarning: File 'file:///home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/fiasco/tests/idl/data/goft_26_11_188.497_v8.0.7.asdf' was created with extension URI 'asdf://astropy.org/core/extensions/core-1.5.0' (from package asdf-astropy==0.5.0), which is not currently installed warnings.warn(msg, AsdfWarning) IDL code to produce Fe XI contribution function result: abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR='abundance') ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq') density = 10d / 1d wave_min = 3298482106769867d / 17592186044416d wave_max = 3333666478858699d / 17592186044416d contribution_function = g_of_t(26,$ 11,$ dens=density,$ abund_file=abund_file,$ ioneq_file=ioneq_file,$ index=22527,/quiet,$ wrange=[wave_min, wave_max]) ; Call this function to get the temperature array read_ioneq,ioneq_file,temperature,ioneq,ref /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/asdf/_asdf.py:359: AsdfWarning: File 'file:///home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/fiasco/tests/idl/data/goft_26_14_197.862_v8.0.7.asdf' was created with extension URI 'asdf://astropy.org/core/extensions/core-1.5.0' (from package asdf-astropy==0.5.0), which is not currently installed warnings.warn(msg, AsdfWarning) IDL code to produce Fe XIV contribution function result: abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR='abundance') ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq') density = 10d / 1d wave_min = 6926465858151645d / 35184372088832d wave_max = 6996834602329309d / 35184372088832d contribution_function = g_of_t(26,$ 14,$ dens=density,$ abund_file=abund_file,$ ioneq_file=ioneq_file,$ index=27964,/quiet,$ wrange=[wave_min, wave_max]) ; Call this function to get the temperature array read_ioneq,ioneq_file,temperature,ioneq,ref /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/asdf/_asdf.py:359: AsdfWarning: File 'file:///home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/latest/lib/python3.11/site-packages/fiasco/tests/idl/data/goft_26_16_262.984_v8.0.7.asdf' was created with extension URI 'asdf://astropy.org/core/extensions/core-1.5.0' (from package asdf-astropy==0.5.0), which is not currently installed warnings.warn(msg, AsdfWarning) WARNING: No proton data available for Fe 16. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions] IDL code to produce Fe XVI contribution function result: abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR='abundance') ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq') density = 10d / 1d wave_min = 4608871268660281d / 17592186044416d wave_max = 4644055640749113d / 17592186044416d contribution_function = g_of_t(26,$ 16,$ dens=density,$ abund_file=abund_file,$ ioneq_file=ioneq_file,$ index=1415,/quiet,$ wrange=[wave_min, wave_max]) ; Call this function to get the temperature array read_ioneq,ioneq_file,temperature,ioneq,ref .. rst-class:: sphx-glr-timing **Total running time of the script:** (5 minutes 40.509 seconds) .. _sphx_glr_download_generated_gallery_idl_comparisons_goft_comparison.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: goft_comparison.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: goft_comparison.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_