.. 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 <sphx_glr_download_generated_gallery_idl_comparisons_goft_comparison.py>`
        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-21

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from astropy.visualization import quantity_support
    from packaging.version import Version

    import fiasco

    from fiasco.tests.idl.helpers import read_idl_test_output
    from fiasco.util.setup_db import LATEST_VERSION

    quantity_support()






.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <astropy.visualization.units.quantity_support.<locals>.MplQuantityConverter object at 0x7ce9dd7c0770>



.. GENERATED FROM PYTHON SOURCE LINES 22-25

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 25-64

.. 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('Contribution Function')
        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 65-71

Next, plot the comparison between the contribution function 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 contribution function. 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 71-112

.. code-block:: Python

    goft_files = [
        'goft_20_15_200.972',
        'goft_26_9_171.073',
        'goft_26_9_188.496',
        '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')
    for i, name in enumerate(goft_files):
        idl_result = read_idl_test_output(name, LATEST_VERSION)
        ion = fiasco.Ion((idl_result['Z'], idl_result['iz']),
                         idl_result['temperature'],
                         abundance_filename=idl_result['abundance'],
                         ionization_fraction=idl_result['ionization_fraction'])
        contribution_func = ion.contribution_function(idl_result['density'])
        transitions = ion.transitions.wavelength[ion.transitions.is_bound_bound]
        idx = np.argmin(np.abs(transitions - idl_result['wavelength']))
        # NOTE: fiasco does not include the n_H/n_e ratio
        if Version(idl_result['chianti_idl_version']) <  Version('9'):
            # Prior to v9, the CHIANTI IDL software assumed this ratio was a constant 0.83
            n_H_n_e = 0.83
        else:
            # Later versions use the actual temperature-dependent proton-to-electron ratio
            n_H_n_e = ion.proton_electron_ratio
        goft = contribution_func[:, 0, idx] * n_H_n_e
        line_label = f'{ion.ion_name_roman} {idl_result["wavelength"]:latex_inline}'
        axes = plot_idl_comparison(
            ion.temperature,
            idl_result['contribution_function'],
            goft,
            fig,
            len(goft_files),
            3*i,
            line_label,
        )
        axes[0].legend()
        print(f"CHIANTI database {idl_result['database_version']}")
        print(f"CHIANTI IDL {idl_result['chianti_idl_version']}")
        print(f'IDL code to produce {line_label} contribution function result:')
        print(idl_result['idl_script'])



.. image-sg:: /generated/gallery/idl_comparisons/images/sphx_glr_goft_comparison_001.png
   :alt: Contribution Function, 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

    WARNING: No autoionization or level-resolved radiative recombination data available for Ca 15. Using single-ion model for level populations calculation. To silence this warning, set use_two_ion_model=False. [fiasco.ions]
    /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: divide by zero encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: invalid value encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    CHIANTI database 11.0.2
    CHIANTI IDL 11.0.2
    IDL code to produce Ca XV $200.972 \; \mathrm{\mathring{A}}$ contribution function result:

    
        abundance_subdirs = ['abundance', 'archive']
    
        abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR=abundance_subdirs)
        ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq')

    
        ; Set ioneq_file this way to get around a bug that always causes the GUI picker to pop up
        ; even when the file is specified. Seems to only happen in v9.
        defsysv,'!ioneq_file',ioneq_file
    

        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,$
                                       wrange=[wave_min, wave_max])
        ; Call this function to get the temperature array
        read_ioneq,ioneq_file,temperature,ioneq,ref

    
        defsysv,'!ioneq_file',''
    
    
    WARNING: No autoionization or level-resolved radiative recombination data available for Fe 9. Using single-ion model for level populations calculation. To silence this warning, set use_two_ion_model=False. [fiasco.ions]
    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/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: divide by zero encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: invalid value encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    CHIANTI database 11.0.2
    CHIANTI IDL 11.0.2
    IDL code to produce Fe IX $171.073 \; \mathrm{\mathring{A}}$ contribution function result:

    
        abundance_subdirs = ['abundance', 'archive']
    
        abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR=abundance_subdirs)
        ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq')

    
        ; Set ioneq_file this way to get around a bug that always causes the GUI picker to pop up
        ; even when the file is specified. Seems to only happen in v9.
        defsysv,'!ioneq_file',ioneq_file
    

        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,$
                                       wrange=[wave_min, wave_max])
        ; Call this function to get the temperature array
        read_ioneq,ioneq_file,temperature,ioneq,ref

    
        defsysv,'!ioneq_file',''
    
    
    WARNING: No autoionization or level-resolved radiative recombination data available for Fe 9. Using single-ion model for level populations calculation. To silence this warning, set use_two_ion_model=False. [fiasco.ions]
    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/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: divide by zero encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: invalid value encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    CHIANTI database 11.0.2
    CHIANTI IDL 11.0.2
    IDL code to produce Fe IX $188.496 \; \mathrm{\mathring{A}}$ contribution function result:

    
        abundance_subdirs = ['abundance', 'archive']
    
        abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR=abundance_subdirs)
        ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq')

    
        ; Set ioneq_file this way to get around a bug that always causes the GUI picker to pop up
        ; even when the file is specified. Seems to only happen in v9.
        defsysv,'!ioneq_file',ioneq_file
    

        density = 10d / 1d
        wave_min = 6596929029167645d / 35184372088832d
        wave_max = 6667297773345309d / 35184372088832d

        contribution_function = g_of_t(26,$
                                       9,$
                                       dens=density,$
                                       abund_file=abund_file,$
                                       ioneq_file=ioneq_file,$
                                       wrange=[wave_min, wave_max])
        ; Call this function to get the temperature array
        read_ioneq,ioneq_file,temperature,ioneq,ref

    
        defsysv,'!ioneq_file',''
    
    
    WARNING: No autoionization or level-resolved radiative recombination data available for Fe 11. Using single-ion model for level populations calculation. To silence this warning, set use_two_ion_model=False. [fiasco.ions]
    /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: divide by zero encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: invalid value encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    CHIANTI database 11.0.2
    CHIANTI IDL 11.0.2
    IDL code to produce Fe XI $188.497 \; \mathrm{\mathring{A}}$ contribution function result:

    
        abundance_subdirs = ['abundance', 'archive']
    
        abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR=abundance_subdirs)
        ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq')

    
        ; Set ioneq_file this way to get around a bug that always causes the GUI picker to pop up
        ; even when the file is specified. Seems to only happen in v9.
        defsysv,'!ioneq_file',ioneq_file
    

        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,$
                                       wrange=[wave_min, wave_max])
        ; Call this function to get the temperature array
        read_ioneq,ioneq_file,temperature,ioneq,ref

    
        defsysv,'!ioneq_file',''
    
    
    WARNING: No autoionization or level-resolved radiative recombination data available for Fe 14. Using single-ion model for level populations calculation. To silence this warning, set use_two_ion_model=False. [fiasco.ions]
    /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: divide by zero encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: invalid value encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    CHIANTI database 11.0.2
    CHIANTI IDL 11.0.2
    IDL code to produce Fe XIV $197.862 \; \mathrm{\mathring{A}}$ contribution function result:

    
        abundance_subdirs = ['abundance', 'archive']
    
        abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR=abundance_subdirs)
        ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq')

    
        ; Set ioneq_file this way to get around a bug that always causes the GUI picker to pop up
        ; even when the file is specified. Seems to only happen in v9.
        defsysv,'!ioneq_file',ioneq_file
    

        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,$
                                       wrange=[wave_min, wave_max])
        ; Call this function to get the temperature array
        read_ioneq,ioneq_file,temperature,ioneq,ref

    
        defsysv,'!ioneq_file',''
    
    
    WARNING: No proton data available for Fe 16. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions]
    WARNING: No proton data available for Fe 17. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions]
    /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: divide by zero encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    /home/docs/checkouts/readthedocs.org/user_builds/fiasco/conda/stable/lib/python3.12/site-packages/astropy/units/quantity.py:648: RuntimeWarning: invalid value encountered in divide
      result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    CHIANTI database 11.0.2
    CHIANTI IDL 11.0.2
    IDL code to produce Fe XVI $262.984 \; \mathrm{\mathring{A}}$ contribution function result:

    
        abundance_subdirs = ['abundance', 'archive']
    
        abund_file = FILEPATH('sun_coronal_1992_feldman_ext.abund', ROOT_DIR=!xuvtop, SUBDIR=abundance_subdirs)
        ioneq_file = FILEPATH('chianti.ioneq', ROOT_DIR=!xuvtop, SUBDIR='ioneq')

    
        ; Set ioneq_file this way to get around a bug that always causes the GUI picker to pop up
        ; even when the file is specified. Seems to only happen in v9.
        defsysv,'!ioneq_file',ioneq_file
    

        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,$
                                       wrange=[wave_min, wave_max])
        ; Call this function to get the temperature array
        read_ioneq,ioneq_file,temperature,ioneq,ref

    
        defsysv,'!ioneq_file',''
    
    





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (2 minutes 39.328 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 <goft_comparison.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: goft_comparison.py <goft_comparison.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: goft_comparison.zip <goft_comparison.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_