Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 92.52%. Comparing base (f7e2560) to head (764e29b).

@@            Coverage Diff             @@
##             main     #282      +/-   ##
==========================================
+ Coverage   92.38%   92.52%   +0.13%     
==========================================
  Files          38       37       -1     
  Lines        2904     2915      +11     
==========================================
+ Hits         2683     2697      +14     
+ Misses        221      218       -3     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

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.

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

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?

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

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