sunpy / sunkit-spex

:construction: Under Development :construction: A package for solar X-ray spectroscopy
BSD 3-Clause "New" or "Revised" License
22 stars 26 forks source link

Changing CHIANTI files for thermal model? #99

Open ianan opened 1 year ago

ianan commented 1 year ago

Describe the feature

The CHIANTI files used for the thermal model are v7.1 sav files, i.e. sunxspex/sunxspex/thermal.py, whereas the OSPEX model has been using V9.0.1 geny files for a while. The latter files does change the lines a bit, i.e comparison model spectra but given the energy range and resolution of some of the spectra fitted by sunxspex it might not make a huge difference. The change to the current CHIANTI v10.0.2 is more minor update to v9.0.1.

Note also that the soon to be released CHIANTI 10.2 is changing the default coronal abundances and moving the older abundance files. That might cause problems in the generation of the files in ssw, and possibly the running of ospex/fvth model in ssw (not a sunxspex problem). But longterm would be good to have latest CHIANTI and updated abundances (possibly as the default?). Both sunxspex and ospex have these CHIANTI files saved with unity abundance so it can be chosen at time of model call.

Proposed solution

I can try and generate v10.0.2 files in the same format as the current v7.1 ones in ssw, and see how they change the model in sunxspex itself.

Once v10.2 is out can explore what want to do with those files - will likely impact OSPEX fvth model, so might want to replicate what it does.

ianan commented 1 year ago

Got the generation working and created an example using the thermal model in sunxspex with default v7.1 and new v10.0.2 files, i.e. https://github.com/ianan/fvth_stuff/blob/main/python/test_thermal_sav.ipynb

For resolution of instruments using (just sub-keV) at worst ~20% different at some lines depending on T.

Note that new v10.0.2 file generated using fewer temperature bins (41 instead of 300 or 750) and over narrower range (logT 6-8 instead of 6-9) as otherwise line file would be too big (100s MB instead of 10s MB).

settwi commented 1 year ago

@ianan can you make a pull request that updates the .sav file download in the code?

ianan commented 1 year ago

@settwi Can do the PR but tried to use the new files for fitting and sunxspex crashes with IndexError: index 84 is out of bounds for axis 0 with size 84 in sunxspex/thermal.py looks like how the lines are weighted and added to the continuum, line 677. Something no doubt to do with the different temperature binning (and maybe more lines).

Example of thermal using existing v7.1 is here: https://github.com/ianan/fvth_stuff/blob/main/python/test_v7_fit.ipynb Example of thermal using my v10 (which breaks) is here: https://github.com/ianan/fvth_stuff/blob/main/python/test_v10_fit.ipynb

settwi commented 1 year ago

I'll take a peek @ianan. Maybe @DanRyanIrish can provide some support too. I think he wrote the thermal emission code.

-WS

settwi commented 1 year ago

Ok, dumb question: how do you get that pair of .sav files?

KriSun95 commented 1 year ago

@ianan @DanRyanIrish

Hacked fix to get fit working

Change Line 4753 in sunxspex.sunxspex_fitting.fitter.py from _test_e_range = np.arange(1.6, 5.01, 0.04)[:, None] to _test_e_range = np.arange(1.6, 15.01, 0.04)[:, None].

Details

This appears to be a bug in the thermal code that is there without changing the CHIANTI files as well.

Running the thermal model with too few energy bins will cause this index error for some reason. I found this out a while back, I can't remember if this was ever brought up in the past, effectively I try to test that a given model is usable by running a quick test. I run the model with parameter values of 1 and energies of _test_e_range = np.arange(1.6, 5.01, 0.04)[:, None] where I pass np.concatenate((_test_e_range[:-1], _test_e_range[1:]), axis=1) to the model function (see these lines).

This isn't great, I know, but it's there as a check for now and can be dealt with in the re-factor. Anyway, the energy range was chosen to mimic NuSTAR energy bin resolution and starting at 1.6 keV but had to be up to 5.01 keV otherwise the same index error that Iain is seeing now would occur. I.e., even if I just changed the definition of _test_e_range to go to 5.00 keV then the index error would result.

I also managed to produce the same error with the new files and sunxspex.thermal.thermal_emission by changing the engs definition in Iain's Notebook (cell three) from engs=np.arange(1.01,21,0.01) to engs=np.arange(1.01,5,0.01). So it seems like the index issue is in the thermal model and not the fitting where the lowest number of bins it can handle is tied to the grid, I am unsure how though.

ianan commented 1 year ago

Can confirm @KriSun95 fix works - but as discussed might be other issues with f_vth that need to think about (i.e. max energy range shouldn't be limited by continuum file? But min energy should be limited by the lines/continuum files?).

As expected though the fit results for v7 vs v10 are very similar due to the energy resolution of the data we are working with.

Also, @settwi if you haven't found it already the sswidl code to produce the sav files is in the same repo but here https://github.com/ianan/fvth_stuff/blob/main/idl/make_new_ch_files_sunxspex.pro

ianan commented 10 months ago

@KriSun95 Just to check but the temporary bug fix for this hasn't been implemented yet? As L4756 is still _test_e_range = np.arange(1.6, 5.01, 0.04)[:, None]

KriSun95 commented 10 months ago

@ianan, you are correct and apologies this quick fix has been held up. I've now opened a PR with this fix applied in #134.

ianan commented 9 months ago

fyi, I have better versions of the chianti database files to use when fitting the thermal model (either in ospex or sunkit-spex) available here https://github.com/ianan/fvth_stuff/tree/main/better_chxdb The issue is that the default sunkit-spex one uses the old v7 files (missing stuff) and the newer v9/v10 ones (v9 default in ospex) don't have enough temperature bins so can get "clustering" in the T results. These new ones are v10.1 and have enough temperature bins to avoid issues (I think?), and aren't too large a file.