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.

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

Changing CHIANTI files for thermal model? #99

Open
ianan opened this issue May 30, 2023 · 10 comments
Open

Changing CHIANTI files for thermal model? #99

ianan opened this issue May 30, 2023 · 10 comments

Comments

@ianan
Copy link

ianan commented May 30, 2023

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

ianan commented Jun 7, 2023

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

settwi commented Sep 11, 2023

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

@ianan
Copy link
Author

ianan commented Sep 22, 2023

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

settwi commented Sep 23, 2023

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

-WS

@settwi
Copy link
Contributor

settwi commented Sep 23, 2023

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

@KriSun95
Copy link
Collaborator

KriSun95 commented Oct 4, 2023

@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
Copy link
Author

ianan commented Nov 7, 2023

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

ianan commented Jan 11, 2024

@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
Copy link
Collaborator

@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
Copy link
Author

ianan commented Feb 14, 2024

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants