Add continuum contributions to `IonCollection.radiative_loss`

[jwreep]

Fixes #232

I've added f-f, f-b, and 2-p continuum 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.

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 96.79144% with 6 lines in your changes missing coverage. Please review.

Project coverage is 92.69%. Comparing base (7cba8b8) to head (1e97da8). Report is 31 commits behind head on main.

Files with missing lines	Patch %	Lines
fiasco/collections.py	82.85%	6 Missing :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #282      +/-   ##
==========================================
- Coverage   92.77%   92.69%   -0.09%     
==========================================
  Files          37       37              
  Lines        2907     3079     +172     
==========================================
+ Hits         2697     2854     +157     
- Misses        210      225      +15     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[jwreep]

The output from all ions looks reasonable, except for log T > 7. There should be an increase with temperature from free-free emission (I thought?).

Edit: The calculation of radiative losses in CHIANTI uses a wavelength-averaged Gaunt factor. Perhaps this is the cause of the discrepancy.

[wtbarnes]

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?

[jwreep]

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!

[wtbarnes]

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

[jwreep]

I tried to redo the function this morning, and while the f-f does have an uptick, I think the b-b is dominating the total losses.

Here's how I changed the continuum calculation:

        wavelength = np.logspace(-2,4,200) * u.angstrom

        ff = self.free_free(wavelength, **kwargs)
        fb = self.free_bound(wavelength, **kwargs)
        tp = self.two_photon(wavelength, density, **kwargs)

        for i in range(len(density)):
            for j in range(len(self.temperature)):
                rad_loss[j,i] += np.trapz(ff[j,:], wavelength)
                rad_loss[j,i] += np.trapz(fb[j,:], wavelength)
                rad_loss[j,i] += np.trapz(tp[:,j,i], wavelength)

and I broke it apart by components as well:

It looks like f-f is still small compared to the total, which might be because b-b is too large or one or all of the continua are too small. Alternatively, maybe this is right?

[wtbarnes]

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.

[jwreep]

I think I figured it out. The discrepancy might be the integration. Changing the code to ignore the bin width

        ff = self.free_free(wavelength, **kwargs).sum(axis=1)*u.angstrom
        fb = self.free_bound(wavelength, **kwargs).sum(axis=1)*u.angstrom
        tp = self.two_photon(wavelength, density, **kwargs).sum(axis=0)*u.angstrom

I get:

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

rad_loss[k]=rad_loss[k]+this_abund*this_ioneq[k]*factor*sqrt(t[k])*z^2*gaunt[k]

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 $R_{ff} = \int_{\lambda} I_{ff}(\lambda) d\lambda$ ?

[jwreep]

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. ~~I have no idea where the exponential went, though!~~

Edit: The Gaunt factor depends on wavelength as well. Sutherland 1998 gives full derivations of both quantities. See section 2.4.

[jwreep]

Tentatively fixed. Still needs a few unit tests to be added, but the losses compare well to IDL.

This uses the same implementation as the IDL formulation, which is based on Mewe et al 1986 for f-b losses and Sutherland 1998 for f-f. Two-photon is neglected.

[jwreep]

@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.

[wtbarnes]

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?

[jwreep]

No worries -- enjoy Tenerife! I'll add the IDL code when I'm in the office tomorrow.

[jwreep]

Using sun_coronal_2021_chianti.abund:

IDL> ff_rad_loss,temperature,ff 
IDL> fb_rad_loss,temperature,fb      
IDL> rad_loss,temperature,total   
IDL> plot,temperature,total,chars=2,xtit='Temperature',ytit='Loss Rate',/ylog,/xlog,yr=[1e-25,1e-21]
IDL> oplot,temperature,fb,line=1,thick=2     
IDL> oplot,temperature,ff,line=2    
IDL> oplot,temperature,(total-ff-fb),line=3     

[jwreep]

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)?

[wtbarnes]

Apologies again that this has taken so long. I've just added a few tests that compare the IDL and Python outputs and made some minor tweaks to one of the gaunt factor functions as I was getting overflow warnings due to large terms in the exponential at low temperatures. These should now be handled properly.

[jwreep]

I've added the doc strings now and some tests to keep Codecov happy. I would suggest keeping each new collection method public -- sometimes it's helpful to break apart the contributions to different emission processes.

[jwreep]

Codecov is still showing a gap for the try-except statements in the f-f and f-b radiative loss collection functions. I'm not sure if it's possible for either to currently throw an exception. The f-b has some error handling in the corresponding ion class function (returns zeroes), while the f-f doesn't use any files particular to a given ion.

So, I think the options are to remove the try-except statements or keep the gap in codecov. I don't think there's anything wrong with keeping the error handling in case something changes in the future, but it shouldn't matter at present.

[wtbarnes]

Don't worry about the code coverage for those parts. I'm not so worried about that "failure".

[wtbarnes]

Ok I did some minor fixing of the formatting for the new topic guide page. The equations weren't rendering quite right. Once the docs build preview finishes, give that a look over and let me know if it's good to go. Other than that, once the tests pass, I think this is ready.

[jwreep]

The math looks like it rendered fine this time. It should be ready to merge once it finishes the last few tests, then!

[wtbarnes]

Ok I found some more nitpicky formatting fixes, almost all caused by me. I also moved the two-photon continuum explanation to a topic guide as well as I found it a bit overly long for a docstring. I'm declaring this the last round of fixes and actually merging once the tests pass this time!

[wtbarnes]

Thanks as always! And sorry that this once again took forever to get merged.