CHIANTI IDL Comparison: Contribution Functions

Compare the contribution function calculation to that in the CHIANTI IDL routines for a few select ions.

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

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

Define a function for comparing and plotting the outputs from fiasco and IDL. Note that we have precomputed the IDL outputs.

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

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 \(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.

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

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',''