IDL and ChiantiPy gives different results.

[MadFisa]

Hi,

I was generating a lookup table for spectra of individual elements for creating a look up table. While doing the analysis, I found that my results were different from my collaborator. On further investigation, we found that the spectra generated by the IDL scripts he used was different from the ones I generated using the ChiantiPy. I will attach the results below. The plots are generated using mspectrum class with unity abundance. If images looks too small in github, downloading and viewing should help.

This is the continuum for each element plotted with IDL using solid lines and Python with dotted lines. Different colors are different temperatures.

We can see there is a difference in the values after first row of elements. Next plot shows the ratio IDL/Python for the elements. Note that ratio for Cu, V and Sc is in 1e10s . These particular elements show biggest difference and as far as I know, do not have ions in CHIANTI.

Next plot shows line calculation for each elements for a particular temperature with IDL in solid blue and python ones in dotted orange.

Here also we can see some differences. Next plot is the ratio.

I have also attached the codes used to generate the IDL and python spectra.

IDL code

pro getSpec

!PATH = expand_path('+/home/mithun/Bin/chianti/idl/')+':'+ !PATH
use_chianti, '/home/mithun/Bin/chianti/dbase'

temperature=10^[6.5,6.8,7.0,7.2]

edensity=1d10

allelements=['h','he','li','be','b','c','n','o','f','ne','na','mg','al','si','p','s','cl','ar','k','ca','sc','ti','v','cr','mn','fe','co','ni','cu','zn']
ents_Z=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30]

masterlist_file='/home/mithun/Bin/chianti/dbase/masterlist/masterlist.ions'
read_masterlist,masterlist_file,all_ions

n_elem=n_elements(allelements)
n_temperature=n_elements(temperature)

wmin=0.8
wmax=124
wstep=0.01

nwbins=long((wmax-wmin)/wstep)+1
allLineSpec=dblarr(nwbins,n_temperature,n_elem)
allContSpec=dblarr(nwbins,n_temperature,n_elem)

abund_file='/home/mithun/Bin/chianti/dbase/abundance/unity.abund'

restore,'lineSpec_CHIANTI_10p0p2.sav'

for j=0,n_elem-1 do begin

    ion_list=all_ions[where(STRMATCH(all_ions, allelements[j]+'_'+'*') eq 1)]
isothermal, wmin,wmax,wstep,temperature,lmbda,spectrum,list_wvl,list_ident,$
      edensity=edensity,/ergs,sngl_ion=ion_list,abund_name=abund_file,ioneq_name=!ioneq_file,$
      /all  

    allLineSpec(*,*,j)=spectrum     

em_int=make_array(n_temperature,value=1,/double)
temp=temperature
ioneq_name=!ioneq_file
freebound, temp,lmbda,fb, $
           em_int=em_int, abund_file=abund_file, $
           ioneq_file=ioneq_name,Element=allelements[j]
freefree,temp,lmbda,ff, $
         em_int=em_int, abund_file=abund_file, $
         ioneq_file=ioneq_name,Element=allelements[j]
two_photon, temp,  lmbda, two_phot, $
            edensity=edensity, em_int=em_int, abund_file=abund_file, $
            ioneq_file=ioneq_name,Element=allelements[j]
totcont=(fb+ff+two_phot)/1d40*wstep

    allContSpec(*,*,j)=totcont 

endfor

save,lmbda,temperature,allContSpec,file='contSpec_CHIANTI_10p0p2.sav'
save,lmbda,temperature,allLineSpec,file='lineSpec_CHIANTI_10p0p2.sav'

stop
end

Python code

#!/usr/bin/env python3

import ChiantiPy.core as ch
import matplotlib.pyplot as plt
import numpy as np
from scipy.io import readsav

#%% Define params
CALCULATE = True  # Flag to decide whether to do calcuation again or use stored ones.

elements_list = [
    "H",
    "He",
    "Li",
    "Be",
    "B",
    "C",
    "N",
    "O",
    "F",
    "Ne",
    "Na",
    "Mg",
    "Al",
    "Si",
    "P",
    "S",
    "Cl",
    "Ar",
    "K",
    "Ca",
    "Sc",
    "Ti",
    "V",
    "Cr",
    "Mn",
    "Fe",
    "Co",
    "Ni",
    "Cu",
    "Zn",
    "Ga",
    "Ge",
    "As",
    "Se",
    "Br",
    "Kr",
    "Rb",
    "Sr",
    "Y",
    "Zr",
    "Nb",
    "Mo",
    "Tc",
    "Ru",
    "Rh",
    "Pg",
    "Ag",
    "Cd",
    "In",
    "Sn",
]

elements = elements_list[:30]

test_elements = elements[:5]

#%% Load idl data
dLine = readsav("lineSpec_CHIANTI_10p0p2.sav") # Line spectra (nelem,ntemperature,nwave) 
dCont = readsav("contSpec_CHIANTI_10p0p2.sav") # Continuum spectra  (nelem,ntemperature,nwave) 

temperature = np.asarray(dLine["temperature"], dtype=np.float64)
temp_label = [f"logT={np.log10(temperature_i):.2f}" for temperature_i in temperature]
wvl = np.asarray(dLine["lmbda"], dtype=np.float64)
dwvl = np.concatenate([np.diff(wvl), [0.01]])  # wavelength bin size

#%% Create spectra
# Doing the same calcuation as idl in python
cont_file_name = "ChiantiPy0p15p0ContSpec_CHIANTI10p02_unfiltered.npy"
line_file_name = "ChiantiPy0p15p0LineSpec_CHIANTI10p02_unfiltered.npy"
if CALCULATE:
    density = 1e10
    py_ContSpec = []
    py_LineSpec = []
    for element in elements:
        s_cont = ch.mspectrum(
            temperature,
            density,
            wvl,
            verbose=1,
            abundance="unity",
            elementList=[element],
            doContinuum=1,
            doLines=0,
            # filter=(dirac_delt,1000)
        )
        s_line = ch.mspectrum(
            temperature,
            density,
            wvl,
            verbose=1,
            abundance="unity",
            elementList=[element],
            doContinuum=0,
            doLines=1,
            # filter=(dirac_delt,1000)
        )

        # Multiply by wavelenght bin size to convert from per AA to  to  total flux in the bin.
        py_ContSpec.append(np.multiply(dwvl, s_cont.Spectrum["intensity"]))
        py_LineSpec.append(np.multiply(dwvl, s_line.Spectrum["intensity"]))

    #%% Save to text
    py_ContSpec = np.array(py_ContSpec)
    py_LineSpec = np.array(py_LineSpec)
    np.save(cont_file_name, py_ContSpec)
    np.save(line_file_name, py_LineSpec)

#%% Load from text
py_ContSpec = np.load(cont_file_name)
py_LineSpec = np.load(line_file_name)
#%% Plot Continuum
fig1, ax1 = plt.subplots(ncols=5, nrows=6, figsize=(16, 9), sharex=True)
fig1.set_tight_layout(True)
fig1.suptitle(r"flux$(erg \; s^{-1} sr^{-1} em^{-1})$ vs wavelength$(\AA^{-1})$")
fig2, ax2 = plt.subplots(ncols=5, nrows=6, figsize=(16, 9), sharex=True)
fig2.set_tight_layout(True)
fig2.suptitle("ratio vs wavelength")
ax1 = ax1.flat
ax2 = ax2.flat
for i, el in enumerate(elements):
    ax1[i].set_title(el)
    ax1[i].plot(
        wvl,
        dCont.allContSpec[i, :, :].T,
    )
    ax1[i].set_yscale("log")
    ax1[i].set_xscale("log")
    ax1[i].set_prop_cycle(None)  # Resetting the color cycle
    ax1[i].plot(
        wvl,
        py_ContSpec[i, :, :].T,
        "--",
    )
    ax1[i].set_ylim([1e-32, 1e-20])
    ratio = dCont.allContSpec[i, :, :] / py_ContSpec[i, :, :]
    ax2[i].plot(
        wvl,
        ratio.T,
    )
    ax2[i].set_title(f"{el}_ratio")
    # ax2[i].set_xscale("log")

fig1.supxlabel("Wavelength (Armstrong)")
fig1.supylabel("flux")
fig1.legend(temp_label)
fig2.supxlabel("Wavelength (Armstrong)")
fig2.supylabel("IDL/Py ratio")
fig2.legend(temp_label)
plt.show()

#%% Plot Lines
fig1, ax1 = plt.subplots(ncols=5, nrows=6, figsize=(16, 9))
fig1.set_tight_layout(True)
fig2, ax2 = plt.subplots(ncols=5, nrows=6, figsize=(16, 9))
fig2.set_tight_layout(True)
ax1 = ax1.flat
ax2 = ax2.flat
temp_idx = 2  # Only plotting temperature with this index for clarity
for i, el in enumerate(elements):
    ax1[i].set_title(el)
    ax1[i].plot(wvl, dLine.alllinespec[i, temp_idx, :].T, label=temperature)
    ax1[i].set_yscale("log")
    ax1[i].set_xscale("log")
    ax1[i].plot(wvl, py_LineSpec[i, temp_idx, :].T, "--", label=temperature)
    ax1[i].set_ylim([1e-32, 1e-20])
    ax1[i].set_prop_cycle(None)  # Resetting the color cycle
    ratio = dLine.alllinespec[i, temp_idx, :] / py_LineSpec[i, temp_idx, :]
    ax2[i].plot(wvl, ratio.T, label=temperature)
    ax2[i].set_title(f"{el}_ratio")
    ax2[i].set_xscale("log")
fig1.suptitle(
    r"flux$(erg \; s^{-1} sr^{-1} em^{-1})$ vs wavelength$(\AA^{-1})$ "
    + f"{temp_label[temp_idx]}"
)
fig1.supylabel("flux")
fig1.legend(["IDL", "Python"])
fig2.supxlabel("Wavelength (Armstrong)")
fig2.supylabel("IDL/Py ratio")
# fig2.legend(np.log10(temperature))
plt.show()

I am not sure what is the reason for this discrepancy, especially in the continuum. I am sorry if its as simple as missing some keyword arguments. Is this discrepancy known?

[MadFisa]

I have also made a link with all the related files and codes, if you want to look into it.
https://drive.google.com/file/d/1V6CpEY4q0SzM_KwWm8_kpFK-6ImmYePJ/view?usp=sharing

[wtbarnes]

Re the continuum results, this also may be related to #118

[MadFisa]

Seems like this is beyond my level of understanding, but please let me know if I can do something to help for investigating this. Also, I don't know if it is relevant, but we first noticed the discrepancy while working in the soft X-ray wavelengths (approximately 0.1 - 10 nm range).

[kdere]

I have done some testing for oxygen ions but it is difficult to sort out a spectrum consisting of many elements and ions. Also, I do not have an IDL license to I have to get someone to do IDL calculations for me. At this point, I am not exactly sure what I will do.

…

________________________________ From: MadFisa ***@***.***> Sent: Sunday, December 18, 2022 4:25 AM To: chianti-atomic/ChiantiPy ***@***.***> Cc: Subscribed ***@***.***> Subject: Re: [chianti-atomic/ChiantiPy] IDL and ChiantiPy gives different results. (Issue #420) Seems like this is beyond my level of understanding, but please let me know if I can do something to help for investigating this. Also, I don't know if it is relevant, but we first noticed the discrepancy while working in hte soft X-ray wavelengths (approximately 0.1 - 10 nm range). — Reply to this email directly, view it on GitHub<https://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fchianti-atomic%2FChiantiPy%2Fissues%2F420%23issuecomment-1356752773&data=05%7C01%7Ckdere%40gmu.edu%7C367810833a9745a2441b08dae0d9be88%7C9e857255df574c47a0c00546460380cb%7C0%7C0%7C638069523053586216%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=HJ%2BqgNbXtyh%2FFKHqwNo91DdItUEm5Y%2F1k85YLUmwI3M%3D&reserved=0>, or unsubscribe<https://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAAFUQVVB2GA252R7LW6I4XLWN3J63ANCNFSM6AAAAAAS7OHLWE&data=05%7C01%7Ckdere%40gmu.edu%7C367810833a9745a2441b08dae0d9be88%7C9e857255df574c47a0c00546460380cb%7C0%7C0%7C638069523053742479%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=FYgOrf9gIhF7oxmtYqbp1TGV1mVN%2BR2t2JO7S7rT2cE%3D&reserved=0>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>