-
Notifications
You must be signed in to change notification settings - Fork 15
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
Add continuum contributions to IonCollection.radiative_loss
#282
base: main
Are you sure you want to change the base?
Conversation
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## main #282 +/- ##
==========================================
+ Coverage 92.77% 93.20% +0.42%
==========================================
Files 37 37
Lines 2907 3017 +110
==========================================
+ Hits 2697 2812 +115
+ Misses 210 205 -5 ☔ View full report in Codecov by Sentry. |
Hmm do you have a comparison of the IDL and Python cases? I definitely agree that you should see the uptick in the losses at high temperatures. I had forgotten that IDL computes the radiative loss contribution from the continuum in a slightly different way. I would think these should be consistent though... Also, if you just plot the wavelength-averaged ff contribution, does that show an increase at high temperatures? |
I haven't had much time to come back to this since I posted. The CHIANTI User's guide shows the uptick, but that plot is from v6 of the database. I don't know if it's the Gaunt factor causing the discrepancy, but I've tried changing the wavelength range and binning without any change. If there's a programming mistake here, it's a subtle one! |
I wonder if this is a consequence of not including the lower wavelengths that make the dominant contributions at high temperatures? For example, computing just the free-free contribution and averaging over the wavelength dimension, import astropy.units as u
import fiasco
import numpy as np
import matplotlib.pyplot as plt
from astropy.visualization import quantity_support
temperature = np.logspace(5,8,100) * u.K
all_ions = fiasco.IonCollection(*[fiasco.Ion(i, temperature) for i in fiasco.list_ions()])
wavelength = np.logspace(-2,2,100) * u.AA
ff_short_wave = all_ions.free_free(wavelength).mean(axis=1)
wavelength = np.logspace(1,2,100) * u.AA
ff_long_wave = all_ions.free_free(wavelength).mean(axis=1)
with quantity_support():
plt.plot(temperature, ff_short_wave, label=r'$\lambda_{min}=0.01$ Å')
plt.plot(temperature, ff_long_wave, label=r'$\lambda_{min}=10$ Å')
plt.xscale('log')
plt.yscale('log')
plt.legend() gives |
That plot almost looks like the FF and FB are missing some scaling factor. I'm surprised that they are so far below the bound-bound losses as well as being comparable to the two-photon losses at most temperatures. |
I think I figured it out. The discrepancy might be the integration. Changing the code to ignore the bin width
This isn't correct (fb > ff), but there's something about the integration that's not working. In IDL ff_rad_loss.pro, the calculation is
which is just a straight sum over all the ions (no integration in wavelength space). So, I guess the question is why does integration not seem to work? Shouldn't the total radiation for e.g. free-free be given by something like |
The f-f emissivity according to technical report 11 is ![]() CHIANTI's technical report #9 addresses the radiative losses. They neglect two-photon's contributions to the total losses. The equation for free-free (2) doesn't appear to be a direct integration of the emissivity. ![]() Note the exponential factor has disappeared. This is perhaps why a direct numerical integration is not producing the correct result. Edit: The Gaunt factor depends on wavelength as well. Sutherland 1998 gives full derivations of both quantities. See section 2.4. |
@wtbarnes, this should be good to go! I've followed the implementation in the IDL version, which doesn't integrate the f-f and f-b functions directly, but implements separate functions for radiative losses. Two-photon emission is neglected in this. I pinned the numpy version < 2 for now. One small note: the implementation of the f-b Gaunt factor in IDL uses a polynomial approximation to an infinite sum (Eq. 14 of Mewe et al 1986). Scipy has the appropriate function for the analytic solution, so I've used that instead. The plot below compares the approximation to the analytic solution. It should be a negligible difference but could cause a small discrepancy with IDL's output. We could make the approximation an option, if wanted. ![]() |
Thanks Jeff! This is great and sorry for being so unresponsive on the review front. I will try to get to this either today or early next week. Re: your last comment, I think using the scipy implementation is better. I did something similar for the Gauss-Laguerre quadrature in the case of the ionization rates. I'm happy to add an IDL comparison test here so that we can check the output as we do any refactoring. Would you mind pasting here snippet of IDL code you were running when you were comparing against the CHIANTI IDL results? |
No worries -- enjoy Tenerife! I'll add the IDL code when I'm in the office tomorrow. |
The calculation is extremely slow, but seems to speed up significantly without the warnings. Perhaps it would be useful to have a keyword to suppress the warnings (or vice versa)? |
Fixes #232
I've added f-f and f-b
, and 2-pcontinuum emission to the radiative loss calculation.It appears to be working, but I'd like to do a bit more testing.One potential issue: the calculation currently assumes a fixed bin width. It might be nice to generalize and do the integral properly to allow for non-fixed bins, particularly since it can be integrated over a very large range.