chianti-atomic / ChiantiPy

ChiantiPy is a python package to calculate the radiative properties of astrophysical plasmas based on the CHIANTI atomic database
63 stars 32 forks source link

Ability to give custom abundance files #350

Closed MadFisa closed 2 years ago

MadFisa commented 2 years ago

Hi,

I might have missed it in the documentation, but is there a way to give custom abundances as input for calculations? Maybe from a file in the same format as those provided by default? In particular, I need to change the elemental abundances in steps and get the corresponding spectra.

kdere commented 2 years ago

Hi, that you do the trick regards, Ken

MadFisa commented 2 years ago

hi @kdere ,

I think your comment got clipped, probably the code didn't make it through.

kdere commented 2 years ago

Hi, yes, it should have read " that should do the trick" so just copy one of the files in the $XUVTOP/abund directory and edit it.

depending on what you are doing you can specify the abundance file you want to use

for example, if you look in the module page in the documentation, it shows the for the bunch and spectrum classes you can specify the abundance keyword to the new file that you can create.

hope this helps,

sorry for the typo. Ken

MadFisa commented 2 years ago

Thanks for the reply. This is what I have been doing so far. I created a file called custom.abund in the $XUVTOP/abund and use that in the spectrum using abundance keyword. But the problem is that I need to change custom.abund repeatedly since I have to step through various abundances.

If I rewrite the custom.abund file and run again it just loads the original custom.abund without reading the new file. It will read the newfile if I reload the chiantiPy package by exiting out of the python and running it again.

This is a bit of problem for me as I have to run the code for a lot of different abundance values and unable to do so without exiting python makes it a bit difficult. Is there anyway to get around this problem?

MadFisa commented 2 years ago

Here is a sample code that will reproduce the problem. It uses a template to make new abundances files and then try to create spectrum for 3 He abundance values. But for some reason, all 3 runs ends up using the same abundance value.

abund_spec.py

import ChiantiPy.core as ch
import numpy as np
import matplotlib.pyplot as plt
import ChiantiPy.tools.filters as chfilters
import pickle
from string import Template
import os

logT = 7.20
temp = 10**7.2
dens=1.0e+9
em = 6.5e26
conversion_factor = 12.398 # wvl in Armstrong to energy in keV
energy_kev = np.linspace(0.4,10.,4000) # Wavelength in keV
wvl=conversion_factor/energy_kev #in nm

try:
    if os.environ['XUVTOP']:
        chianti_dir = os.environ['XUVTOP']
        print(f"Chianti directory is {chianti_dir}")
except KeyError:
    print("Chianti environment variable not set")
abund_file_name = 'custom'

#%% Create abundance file

#These are from sun_coronal_2012_schmelz_ext.abund
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"]
default_abundances =[12.00,10.78, 1.36, 1.66, 3.03, 8.35, 7.71,
                    8.61, 4.41, 7.90, 6.62, 7.87, 6.79, 7.86, 5.49,
                    7.23, 5.10, 6.35, 5.44, 6.64, 3.43, 5.26, 4.33, 
                    5.98, 5.70, 7.85, 5.25, 6.55, 4.54, 4.92]

abund_dictionary = {e:ab for e,ab in zip(elements_list,default_abundances)}

He_list = [11,10,9]
spectrum_list =[]
#Loading a template file of abundance
with open('template.abund','r') as fh:
    templ = Template(fh.read())

#Creating abundance files for various values and trying to run with that file
for He_i in He_list:
    # abund_dictionary['Al'] = 6.2
    abund_dictionary['He'] = He_i
    abund_file_cont = templ.substitute(abund_dictionary)
    with open(f"{chianti_dir}/abundance/{abund_file_name}.abund",'w') as fh:
        fh.write(abund_file_cont)
    #Calculating continuum
    s = ch.spectrum(temp, dens, wvl, doContinuum=1, em=em,minAbund=1.e-4,verbose=1,ionList=['fe_25','ca_19','ca_20','ar_18','si_15','si_14','si_13'],abundance= abund_file_name)
    spectrum_list.append(s)

#%% Printing 
print('The abundances for each run is')
for i,spectrum_i in enumerate(spectrum_list):
    print(f'The values for run: {i} is \n{spectrum_i.AbundAll}')

template.abund

 1   $H    H
 2   $He   He
 3   $Li   Li
 4   $Be   Be
 5   $B    B
 6   $C    C
 7   $N    N
 8   $O    O
 9   $F    F
10   $Ne   Ne
11   $Na   Na
12   $Mg   Mg
13   $Al   Al
14   $Si   Si
15   $P    P
16   $S    S
17   $Cl   Cl
18   $Ar   Ar
19   $K    K
20   $Ca   Ca
21   $Sc   Sc
22   $Ti   Ti
23   $V    V
24   $Cr   Cr
25   $Mn   Mn
26   $Fe   Fe
27   $Co   Co
28   $Ni   Ni
29   $Cu   Cu
30   $Zn   Zn
 -1
%filename: unity.abund

when read, all of the abundances will be 1.0.  This may be helpful in some cases.

% produced for the CHIANTI atomic database by Ken Dere , May 2015
 -1

You can see that if you run the above code, the abundances used in the all three runs is the same.

Am I making this more complicated than it should be? Is there something I am missing here?

MadFisa commented 2 years ago

Going through ChiantiPy code, it seems that the a list of abundance files is created when the ChiantiPy.tools.data is imported https://github.com/chianti-atomic/ChiantiPy/blob/b74d06d883f865922be0f10747faead96561f9d7/ChiantiPy/tools/data.py#L56-L59

My understanding is that, this list will only have the contents when it was imported and will not reflect any changes made after the import. This might be what is causing the problem in the above code. Is there any way to force the list to update again?

MadFisa commented 2 years ago

I managed to get it working by forcing a reload of chdata after modifying the abund files using importlib.reload(chdata). But this seems a bit hacky way of doing it. My suggestion providing a way for the user to give abundance file by giving a arbitrary path would be a better way to do it. There are all needed code is already seems to be in place. Maybe make it in such a way that if the abund file is not in the abund list, also look at the string as path to a file before bringing up the gui menu?

kdere commented 2 years ago

Hi, this is a good idea. I will see what can be done. It should not be too difficult Ken


From: MadFisa @.> Sent: Wednesday, June 1, 2022 8:04 AM To: chianti-atomic/ChiantiPy @.> Cc: Kenneth P Dere @.>; Mention @.> Subject: Re: [chianti-atomic/ChiantiPy] Ability to give custom abundance files (Issue #350)

I managed to get it working by forcing a reload of chdata after modifying the abund files using importlib.reload(chdata). But this seems a bit hacky way of doing it. My suggestion providing a way for the user to give abundance file by giving a arbitrary path would be a better way to do it. There are all needed code is already seems to be in place. Maybe make it in such a way that if the abund file is not in the abund list, also look at the string as path to a file before bringing up the gui menu?

— Reply to this email directly, view it on GitHubhttps://secure-web.cisco.com/1TJEqNRswSNBVaj-Efm0mB9CjMwWvXqTHxHTbGEb7bUv65ztbnb45frSV6tTatSLVVrOwFKNxNTF2yiiFNbvzK8evubFfI4sDajlSqCYyzkUob3z2ceFSIi2bjN3vJPohOk1h7u4ZUw54m90dynWe0DE49D_31Vpwb2s9TCPvtQ90d4ax-EgTLQIJTCW0Am1OwZkyfJ-ilSaEMOGRIBOUcLGWcphDMv9qVg8sy2fLp9rikcvNm5uivyG_TN028Nnmm9jW5glO-bQ_TSXMKc6JnpCCedcN_hA51843sFAhkfPCwsvrGR8kyX_IOtt7jD9sD1HXDSueRX_K1gdFQ09nOssovMcmOCFoyk39j5f2Y3dbTDmIFKSiyGByGh-BqttbY6D7aCDi6ZvDv3Ebps5UD0FDOI67l8E_jHNmLMPYUPsVXSGBHKpgpCJ4mtLTAODv/https%3A%2F%2Fgithub.com%2Fchianti-atomic%2FChiantiPy%2Fissues%2F350%23issuecomment-1143519252, or unsubscribehttps://secure-web.cisco.com/1OJfTEvqY5ZPxTWC-wQLkHuna2N3I72M2atRlfIj8MYAoeHV0w1mAaYUP-9-7Y-XpCSFtu-l5JepxjDqQgV7D_oa42Ir0WJ4SBSr9iyNDjTFy0YBdbRcPwd6MjgXGo8Oh9GE9P_l4vt-rKfbnULck3EL348uPwO_VK8h5viupSQBbyNJTqRumFQf_IMWDge53fcoI7AcAzJ9pjmxZG0-m7g_tiPsWhWLI5abP! %20ExJKq-4au0hNpVqHHwzf0GBgd32hBl0D_zUOHtGlAl2rvi87SyGSMtmfd5A2-vCH-bOFJTNi-3GXZACpUTlUg_d9u14O0MxMV1W4O9xLirha5-a2zumouN5qH1XEDjNC6ZO7aaeacnITO7xQX4eIYkw2NwMKuqvNv3pxVFM2L76VqjqsemyUicuPJtzjdqZEatrMz1yT4-XTfDQrqMlXKx27GvQk/https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAAFUQVTL67YN5UD6WDF5DV3VM5GT3ANCNFSM5XMISXNQ. You are receiving this because you were mentioned.Message ID: @.***>