Fold the temperature response on the spectra wavelengths

changelog/231.bugfix.rst

@@ -0,0 +1 @@
+`~sunkit_instruments.response.SourceSpectra.temperature_response` now interpolates the wavelength response onto the wavelengths of the spectra, instead of the spectra onto the channel wavelengths, so that spectra sampled more finely than the channel keep their flux.

sunkit_instruments/response/tests/test_channel.py

@@ -90,3 +90,17 @@ def test_temperature_response(fake_channel, fake_spectra, obstime):
     temp_response = fake_spectra.temperature_response(fake_channel, obstime=obstime)
     assert isinstance(temp_response, u.Quantity)
     assert temp_response.shape == fake_spectra.temperature.shape
+
+
+@pytest.mark.parametrize("offset", [0, 0.5])
+def test_temperature_response_line(fake_channel, offset):
+    # A unit-integral line on (0) or between (0.5) channel wavelengths gives the response at the line
+    line = fake_channel.wavelength[50] + offset * np.diff(fake_channel.wavelength)[50]
+    delta = 0.02 * u.AA
+    wavelength = line + np.arange(-100, 101) * delta
+    data = np.zeros((1, wavelength.size)) * u.Unit("photon cm3 s-1 sr-1 Angstrom-1")
+    data[0, 100] = u.Quantity(1, "photon cm3 s-1 sr-1") / delta
+    expected = np.interp(line, fake_channel.wavelength, fake_channel.wavelength_response())
+    for order in [slice(None), slice(None, None, -1)]:
+        spectra = SourceSpectra([1e6] * u.K, wavelength[order], data[:, order])
+        assert u.allclose(spectra.temperature_response(fake_channel), expected * u.Quantity(1, "photon cm3 s-1 sr-1"))

sunkit_instruments/response/thermal.py

@@ -165,13 +165,14 @@ def temperature_response(
                                    attrs={'unit': 'Angstrom'}),],
             attrs={'unit': wave_response.unit.to_string()}
         )
-        spec_interp = self._da.interp(
-            wavelength=da_response.wavelength,
+        # Integrate on the spectra wavelengths so that narrow lines keep their flux
+        response_interp = da_response.interp(
+            wavelength=self._da.wavelength,
             kwargs={'bounds_error': False, 'fill_value': 0.0},
         )
-        spec_interp_weighted = spec_interp * da_response
-        temp_response = spec_interp_weighted.integrate(coord='wavelength')
-        final_unit = (u.Unit(spec_interp.unit)
+        spec_weighted = (self._da * response_interp).sortby('wavelength')
+        temp_response = spec_weighted.integrate(coord='wavelength')
+        final_unit = (u.Unit(self._da.unit)
                       * u.Unit(da_response.attrs['unit'])
-                      * u.Unit(spec_interp_weighted.wavelength.unit))
+                      * u.Unit(self._da.wavelength.unit))
         return u.Quantity(temp_response.data, final_unit)