Temperature response functions

[nabobalis]

In GitLab by @wtbarnes on Sep 11, 2019, 08:42

Now that we have a way to calculate the wavelength response functions, for each channel, we should think about how to implement the temperature response functions.

To get the temperature response requires calculating emission lines from CHIANTI. From what I can gather from SSW, the emissivity data is stored as a look up table and then read in when one calls aia_get_response.pro. Ideally, we would include the capability for the user to re-compute all of these emission lines with different abundances, ionization equilibriums, version of CHIANTI, etc. However, this requires a reliable Python interface to CHIANTI which does not yet exist.

Perhaps in the meantime, we could just use the tables in SSW and swap them out once a better solution comes along.

There was an effort a while back to implement all of this in sunpy that ultimately failed. These efforts are summarized in this now very old PR to the main sunpy repo. It is likely that bits and pieces of this old code can be salvaged and reshaped into the response subpackage here.

[nabobalis]

In GitLab by @markcheung on Sep 11, 2019, 15:31

I think for the time being we would just need to reply on those precomputed emissivity tables.

[nabobalis]

In GitLab by @wtbarnes on Sep 12, 2019, 21:08

I agree. At a later point, those tables can be substituted for something calculated in Python.

[nabobalis]

In GitLab by @wtbarnes on Feb 28, 2020, 12:42

As far as implementing this in pure Python, I have a few ideas I'll record here so they don't get forgotten.

We could have an EmissionModel object which is a subclass of fiasco.IonCollection. This would hold all of the ions that a user wants to include when computing the temperature response functions. This also provides a flexible way for users to specify different abundance sets, ionization equilibria, etc.

The EmissionModel could optionally create a table of level populations for each ion and save them (e.g. in an HDF5 file) as this calculation is very expensive. In pseudocode, maybe this looks something like

class EmissionModel(fiasco.IonCollection):

    def __init__(self, density, ions=None):
        if ions is None:  # All CHIANTI ions by default
            ions = fiasco.get_all_ions()
        self.density = density
        super().__init__(*ions)

    def get_level_populations(self, ion):
        if ion in level_populations_file:
            # read T-by-n-wave array from file
        else:
            pop = ion.level_populations(self.density)
        return pop

[nabobalis]

In GitLab by @markcheung on Jan 26, 2022, 15:05

I'm assuming there is no standard format for an emission model yet. But having this type of interface would make sense eventually.

[nabobalis]

In GitLab by @markcheung on Jan 26, 2022, 16:11

As a test, I have split of aia_V9_fullemiss.genx into two files at https://drive.google.com/drive/folders/1Cuf8qK7H2UYB00cGvpY3msc9ORlPi42D?usp=sharing

1. aia_V9_fullemiss.nc, and
2. aia_V9_lines.nc

The first file has the emission model (summed over all ~800k lines stored in 2., and over continuum) for a range of logTe and wavelengths.

Computing the temperature response function would be simple:

import xarray as xr # or replace with some netcdf-aware module
from aiapy.response import Channel
emiss = xr.load_dataset('aia_V9_fullemiss.nc')
tresp = []
channels = [94,131,171,193,211,304,335]
for ch in channels:
    c = Channel(ch*u.angstrom)
    emiss_matrix = emiss.interp(wave=c.wavelength, method='linear').total.data*u.Unit(emiss.total.units)
    dlambda = (c.wavelength[1:]-c.wavelength[0:-1]).mean()
    tresp.append(emiss_matrix@c.wavelength_response()*dlambda*pixel_solid_angle)
tresp = np.stack(tresp).T

# Plot response functions
plt.figure(figsize=(10,6))
plt.plot(emiss.logte, tresp, linewidth=2, label=['{0:d} $\AA$'.format(ch) for ch in channels])
plt.title('AIA Temperature Response Functions', fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.yscale('log')
plt.ylim(1e-28,1e-23)
plt.xlim(5.0,7.8)
plt.xlabel(emiss.logte.units, fontsize=20)
plt.ylabel(tresp.unit, fontsize=20)
plt.legend(fontsize=20)

We can extend the spectral response gallery example, or have a new one. Thoughts?