Add a CHIANTI 11 source spectrum

aiapy/response/temperature.py

@@ -16,51 +16,80 @@
 
 __all__ = ["get_source_spectra", "temperature_response"]
 
-_EMISSIVITY_URL = (
-    "https://raw.githubusercontent.com/LM-SAL/backup-files/refs/heads/aia/static/aiapy/aia_V9_emissivity.fits"
-)
-_EMISSIVITY_SHA256 = "b926fe19f7daabd596f5b42d8997114223ebf7a5f434ba09d086a9f8cf2bd24b"
+_EMISSIVITY_URL = "https://raw.githubusercontent.com/LM-SAL/backup-files/refs/heads/aia/static/aiapy/"
+_EMISSIVITY_FILES = {"solarsoft": "emissivity_v9", "chianti11": "emissivity_chianti11"}
 
 
-@manager.require("emissivity_v9", [_EMISSIVITY_URL], _EMISSIVITY_SHA256)
-def get_source_spectra():
+@manager.require(
+    "emissivity_v9",
+    [_EMISSIVITY_URL + "aia_V9_emissivity.fits"],
+    "b926fe19f7daabd596f5b42d8997114223ebf7a5f434ba09d086a9f8cf2bd24b",
+    defer_download=True,
+)
+@manager.require(
+    "emissivity_chianti11",
+    [_EMISSIVITY_URL + "aia_chianti11_emissivity.fits"],
+    "90dbd44f455a8a6e9c2f0aca2f32474ad5dd8ea8167851aad6f751ed546c252f",
+    defer_download=True,
+)
+def get_source_spectra(source="solarsoft"):
     r"""
-    The emissivity that SolarSoft uses for the AIA temperature responses.
-
-    This is the summed line and continuum emissivity of ``aia_V9_fullemiss.genx``
-    (CHIANTI 9.0.1, Feldman 1992 extended coronal abundances, ``chianti.ioneq``,
-    constant pressure of :math:`10^{15}` K cm\\ :sup:`-3`, He II 303.786 Å boosted by 20),
-    on a 0.1 Å grid from 1 to 901 Å.
-    SolarSoft converts it to DN with :math:`hc` = 12398 eV Å, so it is rescaled here by
-    12398 eV Å / :math:`hc`; folding it with the exact :math:`hc` that
-    `aiapy.response.Channel` uses then reproduces SolarSoft.
-    ``meta`` records these settings and the ``aia_V9_chiantifix.genx`` correction
-    used by `temperature_response`.
+    The emissivity of an optically thin plasma for the AIA temperature responses.
+
+    Both sources are summed line and continuum emissivities on a 0.1 Å grid from 1 to 901 Å
+    and log T from 4 to 9 in steps of 0.05, with the plasma of SolarSoft's
+    ``aia_V9_fullemiss.genx``: Feldman 1992 extended coronal abundances, the CHIANTI
+    ionization equilibrium, a constant pressure of :math:`10^{15}` K cm\\ :sup:`-3`,
+    and the He II 304 Å doublet boosted by 20 for its optical depth.
+    SolarSoft's V9 file boosts only the weak member of the doublet, because its exact
+    wavelength match for the strong one failed on CHIANTI 9 (V6 and V8 boosted both);
+    ``"chianti11"`` boosts both, so its 304 Å response is about 2.7 times SolarSoft's
+    at He II temperatures.
+    ``meta`` records these settings, and for ``"solarsoft"`` the ``aia_V9_chiantifix.genx``
+    correction used by `temperature_response`.
+
+    Parameters
+    ----------
+    source : {"solarsoft", "chianti11"}, optional
+        ``"solarsoft"`` is the emissivity that SolarSoft uses, from ``aia_V9_fullemiss.genx``
+        (CHIANTI 9.0.1), with each line moved to the nearest grid wavelength.
+        SolarSoft converts it to DN with :math:`hc` = 12398 eV Å, so it is rescaled here by
+        12398 eV Å / :math:`hc`; folding it with the exact :math:`hc` that
+        `aiapy.response.Channel` uses then reproduces SolarSoft.
+        ``"chianti11"`` is computed from CHIANTI 11.0.2 with fiasco, with each line shared
+        between the two grid wavelengths around it, so folding it with a linearly interpolated
+        effective area gives the response at the line's exact wavelength.
+        Between CHIANTI 9 and 11 the Fe XIV line at 93.614 Å moved to 91.273 Å, out of the
+        94 Å passband, which lowers the 94 Å response at log T 6.3 by two thirds; the
+        chiantifix does not apply to it.
 
     Returns
     -------
     `sunkit_instruments.response.SourceSpectra`
     """
-    with fits.open(manager.get("emissivity_v9")) as hdul:
+    with fits.open(manager.get(_EMISSIVITY_FILES[source])) as hdul:
         header = hdul[0].header
         # astype gives native byte order, which xarray needs
         temperature = 10 ** hdul["TEMPERATURE"].data["logte"].astype(float) * u.K
         wavelength = hdul["WAVELENGTH"].data["wavelength"].astype(float) * u.angstrom
         emissivity = hdul["EMISSIVITY"].data.astype(float) * u.Unit("ph cm3 s-1 sr-1 Angstrom-1")
-        emissivity *= (12398 * u.eV * u.angstrom / (h * c)).to_value(u.dimensionless_unscaled)
-        chiantifix = hdul["CHIANTIFIX"].data
         meta = {
             "chianti_version": header["CHIANTI"],
             "abundance": header["ABUNDANC"],
             "ioneq": header["IONEQ"],
             "pressure": header["PRESSURE"] * u.K / u.cm**3,
             "fixes": list(header["HISTORY"]),
-            "source": header["EMISFILE"],
-            "chiantifix": {
+        }
+        if "EMISFILE" in header:  # SolarSoft's file: its hc and its chiantifix
+            emissivity *= (12398 * u.eV * u.angstrom / (h * c)).to_value(u.dimensionless_unscaled)
+            chiantifix = hdul["CHIANTIFIX"].data
+            meta["source"] = header["EMISFILE"]
+            meta["chiantifix"] = {
                 str(channel): row.astype(float) * u.Unit("cm5 DN / (pix s)")
                 for channel, row in zip(chiantifix["channel"], chiantifix["chiantifix"], strict=True)
-            },
-        }
+            }
+        else:
+            meta["source"] = f"fiasco {header['FIASCO']}"
     return SourceSpectra(temperature, wavelength, emissivity, density=meta["pressure"] / temperature, meta=meta)
 
 
@@ -114,7 +143,7 @@ def temperature_response(channel, spectra, obstime=None, *, chiantifix=False):
     response = spectra.temperature_response(channel)
     if chiantifix:
         if "chiantifix" not in spectra.meta:
-            msg = "chiantifix only applies to the SolarSoft emissivity from get_source_spectra()"
+            msg = 'chiantifix only applies to the SolarSoft emissivity from get_source_spectra("solarsoft")'
             raise ValueError(msg)
         if not channel.include_eve_correction:
             msg = "chiantifix was derived for the EVE-normalized response; use include_eve_correction=True"

aiapy/response/tests/test_temperature.py

@@ -1,6 +1,5 @@
 import numpy as np
 import pytest
-from sunkit_instruments.response import SourceSpectra
 
 import astropy.time
 import astropy.units as u
@@ -18,6 +17,11 @@ def spectra():
     return get_source_spectra()
 
 
+@pytest.fixture(scope="module")
+def spectra_chianti11():
+    return get_source_spectra("chianti11")
+
+
 @pytest.mark.parametrize(
     ("channel", "corrected", "expected"),
     # SSW aia_get_response(/temperature,/dn) at log T = 5.6, 6.0, 6.4 and 7.0 with evenorm=0,
@@ -44,10 +48,34 @@ def test_chiantifix_requires_eve(spectra) -> None:
         temperature_response(Channel(94 * u.angstrom), spectra, chiantifix=True)
 
 
-def test_chiantifix_requires_ssw_spectra(spectra) -> None:
-    other = SourceSpectra(spectra.temperature, spectra.wavelength, spectra.data)
+@pytest.mark.parametrize(
+    ("channel", "index", "expected", "rtol"),
+    # sdo.aia.temperature_response's aia_chianti11_temperature_response.ecsv (sun-data/solar-dynamics-observatory#28,
+    # the same plasma, fiasco 0.8.2 and CHIANTI 11.0.2), raw, at log T 6.0, 6.4 and 7.0. Two known differences:
+    # sdo boosts only one member of the He II doublet (see get_source_spectra); He II reaches 335 Å through the
+    # short-wavelength wing of its filter and 94 Å through the 304 Å crosstalk, worth 1.8e-3 of 335 Å at log T 6.0,
+    # 3.9e-4 at 6.4 and 7.9e-4 at 7.0 (and 1.1e-4 of 94 Å at 5.6), so 335 Å is compared at 6.4 and 7.0 only. And sdo
+    # computed the continua on a grid starting at 25 Å, where CHIANTI's free-bound routine (and fiasco's port of it)
+    # leaves out recombination onto the excited levels of the H-like and He-like ions (mostly O VIII, O VII and
+    # N VII here), worth 2.4e-3 of 94 Å at log T 6.4 and under 1e-4 elsewhere. The rest is the fold, about 4e-5.
+    [
+        (94, [40, 60], [1.60677102e-27, 1.91762518e-27], 2e-4),
+        (94, [48], [1.34141879e-28], 3e-3),
+        (335, [48, 60], [4.95033912e-27, 4.59150350e-28], 1.5e-3),
+    ],
+)
+def test_temperature_response_chianti11(spectra_chianti11, channel, index, expected, rtol) -> None:
+    assert spectra_chianti11.meta["chianti_version"] == "11.0.2"
+    assert "chiantifix" not in spectra_chianti11.meta
+    temperature, response = temperature_response(Channel(channel * u.angstrom), spectra_chianti11)
+    assert_quantity_allclose(np.log10(temperature[[40, 48, 60]].to_value(u.K)), [6.0, 6.4, 7.0])
+    assert_quantity_allclose(response[index], expected * u.Unit("cm5 DN / (pix s)"), rtol=rtol)
+
+
+def test_chiantifix_requires_ssw_spectra(spectra_chianti11) -> None:
+    channel = Channel(94 * u.angstrom, include_eve_correction=True)
     with pytest.raises(ValueError, match="get_source_spectra"):
-        temperature_response(Channel(94 * u.angstrom, include_eve_correction=True), other, chiantifix=True)
+        temperature_response(channel, spectra_chianti11, chiantifix=True)
 
 
 def test_temperature_response_obstime_array(spectra) -> None:

changelog/445.feature.rst

@@ -0,0 +1 @@
+`aiapy.response.get_source_spectra` takes ``source="chianti11"`` for an emissivity computed from CHIANTI 11.0.2 with fiasco, with the plasma of the SolarSoft emissivity and each line folded at its exact wavelength.

examples/temperature_response.py

@@ -49,9 +49,18 @@
 temperature, response = temperature_response(channel, spectra, obstime)
 _, response_chiantifix = temperature_response(channel, spectra, obstime, chiantifix=True)
 
+###############################################################################
+# The emissivity can also come from CHIANTI 11, computed with fiasco. Between
+# CHIANTI 9 and 11 the Fe XIV line at 93.6 Å moved out of the 94 Å passband,
+# so the warm 94 Å response is lower still. The chiantifix does not apply to it.
+
+spectra_chianti11 = get_source_spectra("chianti11")
+_, response_chianti11 = temperature_response(channel, spectra_chianti11, obstime)
+
 fig, ax = plt.subplots()
-ax.plot(temperature, response, label="94 Å")
-ax.plot(temperature, response_chiantifix, label="94 Å with chiantifix")
+ax.plot(temperature, response, label="94 Å, SolarSoft (CHIANTI 9)")
+ax.plot(temperature, response_chiantifix, label="94 Å, SolarSoft with chiantifix")
+ax.plot(temperature, response_chianti11, label="94 Å, CHIANTI 11")
 ax.set_xscale("log")
 ax.set_yscale("log")
 ax.set_xlim(1e5, 1e8)