ch_ss like synthetic spectra with chiantipy

[jacobdparker]

Hello all,

I am attempting to reproduce the results from ch_ss using Chiantipy under the following conditions:

Constant Pressure = 1e+15
DEM = 'quiet_sun.dem'
abundance = 'sun_coronal_2012_schmelz'
wavelength_range = 580 - 630 Angstroms

So far I have tried the bunch and spectrum classes, but I am struggling to get the information I need out of those objects (pouring over the documentation and code now).

A few things I have noticed so far:
I can't find a constant pressure option, so I have tried to calculate it and input a density array,
I can't find a DEM option, so I have loaded in the temperatures and em manually,
ion_list keyword works great for bunch, but throws and error in spectrum,
bunch.Intensity seems to ignore the wavelength range input into ch.bunch

So far bunch.intensityList() seems to be the closest to what I want, but I'm not sure the best way to integrate over temperature since it takes the index input. I've sorted bunch.Intensity to only contain the lines within wavelength the range I want and get the same number of lines as ch_ss. The brightest 10 lines in Intensity['integrated'] seem to at least be the correct lines, but the intensity and relative intensities are different from ch_ss.

I feel like I am close but am missing a few key things about these objects, any help is appreciated. Here is my test so far:

import ChiantiPy.core as ch
import numpy as np
import matplotlib.pyplot as plt
import pandas
import astropy.units as u

# The goal of this test file is to match the results of Chiantipy to that of ch_ss in IDL to make sure we know how
# everything works.

if __name__ == '__main__':
    dem_file = '/home/jake/chianti/dem/quiet_sun.dem'
    dem = pandas.read_csv(dem_file, sep=' ', skipinitialspace=True, skipfooter=9, names=['logT', 'EM'])
    wvl_range = [580, 630]
    temperature = 10 ** dem['logT'].to_numpy()
    pressure = 1e15
    dens = pressure / temperature
    em = 10 ** dem['EM'].to_numpy()
    ion_list = ['o_3', 'o_4', 'o_5', 'he_1', 'mg_10']
    abund_file = 'sun_coronal_2012_schmelz'
    minabund = 7.4e-5

    test_bunch = ch.bunch(temperature, dens, wvl_range, em=em, abundance=abund_file,
                          allLines=0,
                          verbose=True, ionList=ion_list,
                          # minAbund=minabund,  # Note, minAbund will override ionList
                          )

    # test_bunch.intensityList(wvlRange=wvl_range,top=5,)
    # print(test_bunch.Intensity)
    wvl = test_bunch.Intensity['wvl']
    int = test_bunch.Intensity['integrated']
    mask = (wvl < wvl_range[1]) & (wvl > wvl_range[0])
    wvl = wvl[mask]
    print(wvl.shape)
    int = int[mask]
    sort = int.argsort()
    sort = sort[::-1]
    print(int[sort][:10])
    print(wvl[sort][:10])
    
    # test_spectrum = ch.spectrum(temperature, dens, wavelength=wvl_range, abundance=abund_file, verbose=1, allLines=0,
    #                             # ionList=ion_list, #ionList did not work, threw an error
    #                             minAbund=minabund, doContinuum=False)

    # plot_test = test_spectrum.spectrumPlot(integrated=True) # Seems to do nothing
    # plt.show()

Thanks,
Jake

[jacobdparker]

Update:
It appears that Intensity['integrated'] is a simple sum along the temperature axis of Intensity['intensity']. When using an integration method like np.trapz() that accounts for the uneven spacing in temperature I get answers matching ch_ss.

So I guess my question now is, am I doing this the hard way or is there a functionality in ChiantiPy that uses the DEM files that I am missing? Also, is there a different way I should be inputting the wavelength range so it is used by bunch?

[kdere]

ChiantiPy does not currently use the .dem files. If you convert the to emission measures, the can be used a the em input.
You might have to do some digging into the IDL code to figure this out. I believe that what IDL uses to multiply the emissivity by is the dem times dT. dT is often 10**0.1

I don't understand why you need a different way to input the temperature into bunch

hope this is of some help
Ken

[jacobdparker]

Thanks for the response Ken. If ChiantiPy isn't yet set up for DEMs then multiplying by dT prior to summing is what I'll do.

Would it make sense to update Intensity['integrated'] to account for non-uniform temperature spacing?

I don't understand why you need a different way to input the temperature into bunch

Sorry, I was asking about the wavelength input. I still don't understand why bunch.Intensity ignores the wavelength range input?

[kdere]

I will check into bunch for the wavelength range input. Thanks for pointing this out

[kdere]

I think I have fixed this in the github repository.
Ken