diff --git a/changelog/231.bugfix.rst b/changelog/231.bugfix.rst new file mode 100644 index 00000000..4ccc53d4 --- /dev/null +++ b/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. diff --git a/sunkit_instruments/response/tests/test_channel.py b/sunkit_instruments/response/tests/test_channel.py index 69a95f80..9dc3aaa7 100644 --- a/sunkit_instruments/response/tests/test_channel.py +++ b/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")) diff --git a/sunkit_instruments/response/thermal.py b/sunkit_instruments/response/thermal.py index 7de91ee5..8e8e0b89 100644 --- a/sunkit_instruments/response/thermal.py +++ b/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)