interoplate contamination onto a dense grid of wavelengths

xrtpy/response/effective_area.py

@@ -11,6 +11,7 @@
 
 import astropy.time
 import numpy as np
+import scipy.interpolate
 import scipy.io
 import sunpy.io.special
 import sunpy.time
@@ -523,29 +524,34 @@ def _CCD_contamination_transmission(self):
         return np.array([abs(transmittance[i] ** 2) for i in range(4000)])
 
     @property
-    def channel_wavelength(self):
+    def wavelength(self):
         """Array of wavelengths for every X-ray channel in Angstroms (Å)."""
-        return Channel(self.name).wavelength
+        _wave = self._channel.wavelength.to_value("AA")
+        delta_wave = 0.01
+        return np.arange(_wave[0], _wave[-1], delta_wave) * u.Angstrom
 
     @property
     def channel_geometry_aperture_area(self):
         """XRT flight model geometry aperture area."""
-        return Channel(self.name).geometry.geometry_aperture_area
+        return self._channel.geometry.geometry_aperture_area
 
     @property
     def channel_transmission(self):
         """XRT channel transmission."""
-        return Channel(self.name).transmission
+        return np.interp(
+            self.wavelength, self._channel.wavelength, self._channel.transmission
+        )
+
+    def _contamination_interpolator(self, x, y):
+        return np.interp(self.wavelength.to_value("Angstrom"), x, y)
 
     @property
     def _interpolated_CCD_contamination_transmission(self):
         """Interpolate filter contam transmission to the wavelength."""
-        CCD_contam_transmission = np.interp(
-            self.channel_wavelength.to_value("AA"),
+        return self._contamination_interpolator(
             self.n_DEHP_wavelength,
             self._CCD_contamination_transmission,
         )
-        return CCD_contam_transmission
 
     @cached_property
     def _filter_contamination_transmission(self):
@@ -587,12 +593,10 @@ def _filter_contamination_transmission(self):
     @property
     def _interpolated_filter_contamination_transmission(self):
         """Interpolate filter contam transmission to the wavelength."""
-        Filter_contam_transmission = np.interp(
-            self.channel_wavelength.to_value("AA"),
+        return self._contamination_interpolator(
             self.n_DEHP_wavelength,
             self._filter_contamination_transmission,
         )
-        return Filter_contam_transmission
 
     @u.quantity_input
     def effective_area(self) -> u.cm**2:

xrtpy/response/temperature_from_filter_ratio.py

@@ -530,7 +530,7 @@ def calculate_TE_errors(map1, map2, T_e, EM, model_ratio, tresp1, tresp2, Trange
         Narukage's K factor for image 2
     """
 
-    wvl = tresp1.channel_wavelength
+    wvl = tresp1.wavelength
     eVe = tresp1.ev_per_electron
     gain = tresp1.ccd_gain_right
     # (h*c/lambda) * 1/(eV per electron) * 1/gain

xrtpy/response/tests/test_effective_area.py

@@ -105,18 +105,28 @@ def get_IDL_data_files():
     return sorted(filter_data_files)
 
 
+# NOTE: This is marked as xfail because the IDL results that this test compares against
+# are incorrect due to the use of quadratic interpolation in the contamination curves
+# which leads to ringing near the edges in the contamination curve.
+# See https://github.com/HinodeXRT/xrtpy/pull/284#issuecomment-2334503108
+@pytest.mark.xfail
 @pytest.mark.parametrize("filename", get_IDL_data_files())
 def test_effective_area_compare_idl(filename):
     with Path.open(filename) as f:
         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_effective_area = np.loadtxt(filename, skiprows=3)[:, 1] * u.cm**2
+    IDL_data = np.loadtxt(filename, skiprows=3)
+    IDL_wavelength = IDL_data[:, 0] * u.AA
+    IDL_effective_area = IDL_data[:, 1] * u.cm**2
     instance = EffectiveAreaFundamental(filter_name, filter_obs_date)
     actual_effective_area = instance.effective_area()
+    IDL_effective_area = np.interp(
+        instance.wavelength, IDL_wavelength, IDL_effective_area
+    )
     assert u.allclose(
         actual_effective_area,
         IDL_effective_area,
-        atol=1e-2 * u.cm**2,
+        rtol=1e-6,
     )