clean up temperature response and tests

xrtpy/response/temperature_response.py

@@ -4,19 +4,15 @@
 
 from pathlib import Path
 
+import astropy.constants as const
 import numpy as np
 import scipy.io
 from astropy import units as u
-from astropy.constants import c, h
 from scipy import interpolate
 
 from xrtpy.response.channel import Channel, resolve_filter_name
 from xrtpy.response.effective_area import EffectiveAreaFundamental
 
-_c_Å_per_s = c.to(u.angstrom / u.second).value
-_h_eV_s = h.to(u.eV * u.s).value
-
-
 _abundance_model_file_path = {
     "coronal_abundance_path": Path(__file__).parent.absolute()
     / "data/chianti_emission_models"
@@ -158,31 +154,31 @@ def file_spectra(self):
 
     @property
     @u.quantity_input
-    def wavelength(self):
+    def _wavelength_spectra(self):
         """Emission model file wavelength values in Å."""
         return u.Quantity(self._get_abundance_data["wavelength"] * u.Angstrom)
 
     @property
     @u.quantity_input
-    def channel_wavelength(self):
+    def wavelength(self):
         """Array of wavelengths for every X-ray channel in Å."""
-        return u.Quantity((Channel(self.filter_name).wavelength[:3993]) * u.photon)
+        return self._effective_area_fundamental.wavelength
 
     @property
     def focal_len(self):
         """Focal length of the telescope in units of cm."""
-        return Channel(self.filter_name).geometry.geometry_focal_len
+        return self._channel.geometry.geometry_focal_len
 
     @property
     def ev_per_electron(self):
         """Amount of energy it takes to dislodge 1 electron in the CCD."""
-        return Channel(self.filter_name).ccd.ccd_energy_per_electron
+        return self._channel.ccd.ccd_energy_per_electron
 
     @property
     @u.quantity_input
     def pixel_size(self) -> u.cm:
         """CCD pixel size. Units converted from μm to cm."""
-        ccd_pixel_size = Channel(self.filter_name).ccd.ccd_pixel_size
+        ccd_pixel_size = self._channel.ccd.ccd_pixel_size
         return ccd_pixel_size.to(u.cm)
 
     @property
@@ -204,14 +200,12 @@ def spectra(self) -> u.photon * u.cm**3 / (u.sr * u.s * u.Angstrom):
         spectra_interpolate = []
         for i in range(61):
             interpolater = interpolate.interp1d(
-                self.wavelength,
+                self._wavelength_spectra.to_value("AA"),
                 self.file_spectra[i],
                 kind="linear",
             )
-            spectra_interpolate.append(interpolater(self.channel_wavelength))
-        return spectra_interpolate * (
-            u.photon * u.cm**3 * (1 / u.sr) * (1 / u.s) * (1 / u.Angstrom)
-        )
+            spectra_interpolate.append(interpolater(self.wavelength.to_value("AA")))
+        return spectra_interpolate * u.Unit("photon cm3 sr-1 s-1 Angstrom-1")
 
     @u.quantity_input
     def effective_area(self) -> u.cm**2:
@@ -235,27 +229,24 @@ def integration(self) -> u.electron * u.cm**5 / (u.s * u.pix):
         astropy.units.Quantity
             Integrated temperature response in electron cm^5 / (s pix).
         """
-        wavelength = (self.channel_wavelength).value
-        constants = (_c_Å_per_s * _h_eV_s / self.channel_wavelength).value
-        factors = (self.solid_angle_per_pixel / self.ev_per_electron).value
-        effective_area = (self.effective_area()).value
-        dwvl = wavelength[1:] - wavelength[:-1]
-        dwvl = np.append(dwvl, dwvl[-1])
+        constants = const.h * const.c / self.wavelength / u.photon
+        constants *= self.solid_angle_per_pixel / self.ev_per_electron
         # Simple summing like this is appropriate for binned data like in the current
         # spectrum file. More recent versions of Chianti include the line width,
         # which then makes the previous version that uses Simpson's method
         # to integrate more appropriate (10/05/2022)
-        temp_resp_w_u_c = (
-            self.spectra().value * effective_area * constants * factors * dwvl
+        return (
+            self.spectra()
+            * self.effective_area()
+            * constants
+            * np.gradient(self.wavelength)
         ).sum(axis=1)
 
-        return temp_resp_w_u_c * (u.electron * u.cm**5 * (1 / u.s) * (1 / u.pix))
-
     @property
     @u.quantity_input
     def ccd_gain_right(self) -> u.electron / u.DN:
         """Provide the camera gain in electrons per data number."""
-        return Channel(self.filter_name).ccd.ccd_gain_right
+        return self._channel.ccd.ccd_gain_right
 
     @u.quantity_input
     def temperature_response(self) -> u.DN * u.cm**5 / (u.s * u.pix):

xrtpy/response/tests/test_temperature_response.py

@@ -1,89 +1,67 @@
-from datetime import datetime
 from pathlib import Path
 
+import astropy.units as u
 import numpy as np
 import pytest
+from astropy.utils.data import get_pkg_data_filenames
 
 from xrtpy.response.temperature_response import TemperatureResponseFundamental
 
 
-def get_IDL_data_files(model):
-    path = (
-        Path(__file__).parent.parent.absolute()
-        / "tests"
-        / "data"
-        / f"temperature_response_{model}_IDL_testing_files"
-    )
-    filter_data_files = list(path.glob("**/*.txt"))
+def get_IDL_data_files(abundance):
+    filter_data_files = []
+    for dir in get_pkg_data_filenames(
+        f"data/temperature_response_{abundance}_IDL_testing_files",
+        package="xrtpy.response.tests",
+    ):
+        filter_data_files += list(Path(dir).glob("*.txt"))
     return sorted(filter_data_files)
 
 
-def _IDL_raw_data_list(filename):
-    with open(filename) as filter_file:  # noqa: PTH123
-        IDL_data_list = []
-        for line in filter_file:
-            stripped_line = line.strip()
-            line_list = stripped_line.split()
-            IDL_data_list.append(line_list)
-    return IDL_data_list
-
-
-def IDL_test_abundance_name(IDL_data_list):
-    return str(IDL_data_list[1][1])
-
-
-def IDL_test_filter_name(IDL_data_list):
-    return str(IDL_data_list[2][1])
-
-
-def IDL_test_date(IDL_data_list):
-    obs_date = str(IDL_data_list[3][1])
-    obs_time = str(IDL_data_list[3][2])
-
-    day = int(obs_date[:2])
-    month_datetime_object = datetime.strptime(obs_date[3:6], "%b")
-    month = month_datetime_object.month
-    year = int(obs_date[8:12])
-
-    hour = int(obs_time[:2])
-    minute = int(obs_time[3:5])
-    second = int(obs_time[6:8])
-    return datetime(year, month, day, hour, minute, second)
-
-
-def _IDL_temperature_response_raw_data(filename):
-    with open(filename) as filter_file:  # noqa: PTH123
-        IDL_data_list = []
-        for line in filter_file:
-            stripped_line = line.strip()
-            line_list = stripped_line.split()
-            IDL_data_list.append(line_list)
-
-    new_IDL_data_list = [IDL_data_list[i][1] for i in range(5, len(IDL_data_list))]
-    return [float(i) for i in new_IDL_data_list]
-
+filenames = (
+    get_IDL_data_files("coronal")
+    + get_IDL_data_files("hybrid")
+    + get_IDL_data_files("photospheric")
+)
+
+
+@pytest.mark.parametrize("filename", filenames)
+def test_temperature_response(filename):
+    with Path.open(filename) as f:
+        _ = f.readline()
+        abundance = f.readline().split()[1]
+        filter_name = f.readline().split()[1]
+        filter_obs_date = " ".join(f.readline().split()[1:])
+    # NOTE: Annoyingly the date strings use "Sept" instead of "Sep" for "September"
+    filter_obs_date = filter_obs_date.replace("Sept", "Sep")
+    IDL_data = np.loadtxt(filename, skiprows=5)
+    IDL_temperature = IDL_data[:, 0] * u.K
+    IDL_temperature_response = IDL_data[:, 1] * u.Unit("DN cm5 pix-1 s-1")
+
+    instance = TemperatureResponseFundamental(
+        filter_name, filter_obs_date, abundance_model=abundance
+    )
 
-@pytest.mark.parametrize("abundance_model", ["coronal", "hybrid", "photospheric"])
-def test_temperature_response(
-    abundance_model,
-):
-    filenames = get_IDL_data_files(abundance_model)
-    for filename in filenames:
-        IDL_data = _IDL_raw_data_list(filename)
-        filter_name = IDL_test_filter_name(IDL_data)
-        filter_obs_date = IDL_test_date(IDL_data)
-        IDL_temperature_response = _IDL_temperature_response_raw_data(filename)
+    IDL_temperature_response = np.interp(
+        instance.CHIANTI_temperature, IDL_temperature, IDL_temperature_response
+    )
 
-        instance = TemperatureResponseFundamental(
-            filter_name, filter_obs_date, abundance_model=abundance_model
-        )
+    actual_temperature_response = instance.temperature_response()
 
-        actual_temperature_response = instance.temperature_response()
-        atol = actual_temperature_response.value.max() * 0.013
+    # NOTE: there may be small deviations where the response function is very small, likely
+    # due to differences in the interpolation schemes. These are not critical as the response
+    # is effectively zero in these regions anyway.
+    i_valid = np.where(
+        actual_temperature_response > 1e-8 * actual_temperature_response.max()
+    )
 
-        assert np.allclose(
-            actual_temperature_response.value,
-            IDL_temperature_response,
-            rtol=0.028,
-            atol=atol,
-        )
+    # NOTE: The relative tolerance is set comparatively high here because the CCD right gain
+    # values in the IDL and xrtpy codes are explicitly different. See https://github.com/HinodeXRT/xrtpy/pull/76.
+    # If the IDL results files are corrected to account for this updated gain, then this
+    # relative tolerance can be set significantly lower. Setting the gains to be the same,
+    # nearly all cases match within less than 1%.
+    assert u.allclose(
+        actual_temperature_response[i_valid],
+        IDL_temperature_response[i_valid],
+        rtol=5e-2,
+    )