Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

IDL and ChiantiPy gives different results. #420

Open
MadFisa opened this issue Dec 15, 2022 · 4 comments
Open

IDL and ChiantiPy gives different results. #420

MadFisa opened this issue Dec 15, 2022 · 4 comments

Comments

@MadFisa
Copy link

MadFisa commented Dec 15, 2022 •
edited

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.

Continnuum_flux

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.

Continnuum_ratio

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

lines_flux

Here also we can see some differences. Next plot is the ratio.
lines_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
Copy link
Author

MadFisa commented Dec 15, 2022

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
Copy link
Member

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

@MadFisa
Copy link
Author

MadFisa commented Dec 18, 2022 •
edited

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
Copy link
Contributor

kdere commented Dec 19, 2022 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@kdere @wtbarnes @MadFisa