Merge pull request #284 from wtbarnes/test-refactoring

pyproject.toml

@@ -53,8 +53,9 @@ dev = [
   "nox >= 2022.8.7",
 ]
 tests = [
-  "pytest-allclose >= 1.0.0",
+  "pytest >= 8.0.0",
   "pytest-astropy",
+  "pytest-xdist >= 3.6.1",
 ]
 docs = [
   "imageio >= 2.20.0",

xrtpy/image_correction/tests/test_remove_lightleak.py

@@ -1,61 +1,34 @@
 from pathlib import Path
 
+import numpy as np
 import pytest
+from astropy.utils.data import get_pkg_data_path
 from sunpy.map import Map
 
 from xrtpy.image_correction.remove_lightleak import remove_lightleak
 
-
-def get_IDL_data_file():
-    """
-    The XRT composite fits file that been lightleak corrected in
-    IDL have been rename to begin with "ll" referring to lightleak.
-    """
-    directory = (
-        Path(__file__).parent.absolute()
-        / "data"
-        / "light_leak_testing_data_files"
-        / "IDL_lightleak_corrected_data_test_files"
+data_dir = Path(
+    get_pkg_data_path(
+        "data/light_leak_testing_data_files", package="xrtpy.image_correction.tests"
     )
-    data_files = directory.glob("ll_comp_XRT*.fits")
-    return sorted(data_files)
-
-
-IDL_filenames = get_IDL_data_file()
-
-
-def get_composite_data_files():
-    """
-    The XRT composite fits file are no corrected in IDL.
-    These files will be corrected using XRTpy.
-    """
-    directory = (
-        Path(__file__).parent.absolute()
-        / "data"
-        / "light_leak_testing_data_files"
-        / "xrtpy_lightleak_data_test_files"
-    )
-    data_file = directory.glob("comp_XRT*.fits")
-    return sorted(data_file)
-
-
-composite_filenames = get_composite_data_files()
-
-# Using zip as an iterator to pair the data files together. Trouble-free method to use in pytest-parametrize
-data_files = list(zip(IDL_filenames, composite_filenames, strict=False))
-
-
-@pytest.mark.parametrize(("idlfile", "compfile"), data_files)
-def test_lightleak(idlfile, compfile, allclose):
+)
+composite_filenames = sorted(
+    (data_dir / "xrtpy_lightleak_data_test_files").glob("comp_XRT*.fits")
+)
+IDL_filenames = sorted(
+    (data_dir / "IDL_lightleak_corrected_data_test_files").glob("ll_comp_XRT*.fits")
+)
+
+
+@pytest.mark.parametrize(
+    ("idlfile", "compfile"), list(zip(IDL_filenames, composite_filenames, strict=True))
+)
+def test_lightleak(idlfile, compfile):
     IDL_map = Map(idlfile)
     input_map = Map(compfile)
-
     ll_removed_map_xrtpy = remove_lightleak(input_map)
-
-    if input_map.data.shape == (2048, 2048):
-        # Because of rebinning for full resolution images, the match is worse
-        # between IDL created images and XRTpy ones. IDL's method of rebinning
-        # is different from that used by sunpy.
-        assert allclose(ll_removed_map_xrtpy.data, IDL_map.data, atol=0.75)
-    else:
-        assert allclose(ll_removed_map_xrtpy.data, IDL_map.data, atol=1e-5)
+    # Because of rebinning for full resolution images, the match is worse
+    # between IDL created images and XRTpy ones. IDL's method of rebinning
+    # is different from that used by sunpy.
+    atol = 0.75 if input_map.data.shape == (2048, 2048) else 1e-5
+    assert np.allclose(ll_removed_map_xrtpy.data, IDL_map.data, atol=atol)

xrtpy/response/channel.py

@@ -432,6 +432,8 @@ def ccd_gain_left(self) -> u.electron / u.DN:
     @u.quantity_input
     def ccd_gain_right(self) -> u.electron / u.DN:
         """Gain when reading the right port of the CCD."""
+        # NOTE: Value for the right gain in the instrument data files is incorrect.
+        # See https://github.com/HinodeXRT/xrtpy/pull/76
         return u.Quantity(57.5, u.electron / u.DN)
 
     @property

xrtpy/response/effective_area.py

@@ -10,6 +10,7 @@
 
 import astropy.time
 import numpy as np
+import scipy.interpolate
 import scipy.io
 import sunpy.io.special
 import sunpy.time
@@ -522,29 +523,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):
@@ -586,12 +592,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

@@ -535,7 +535,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/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_channel.py

@@ -848,41 +848,6 @@ def test_channel_name2(channel_name):
     assert name == IDL_mirror_name_AUTO
 
 
-@pytest.mark.parametrize("channel_name", channel_names)
-def test_channel_wavelength(channel_name):
-    channel_filter = Channel(channel_name)
-
-    wavelength_length = int(channel_filter.number_of_wavelengths)
-    wavelength = channel_filter.wavelength[:wavelength_length]
-
-    idl_array_length = int(
-        v6_genx_s[_channel_name_to_index_mapping[channel_name]]["LENGTH"]
-    )
-    idl_wavelength_auto = (
-        v6_genx_s[_channel_name_to_index_mapping[channel_name]]["WAVE"][
-            :idl_array_length
-        ]
-        * u.angstrom
-    )
-
-    assert u.allclose(idl_wavelength_auto, wavelength)
-
-    idl_mirror_wavelength_manu = [
-        9.00000,
-        9.10000,
-        9.20000,
-        9.30000,
-        9.40000,
-        9.50000,
-        9.60000,
-        9.70000,
-        9.80000,
-        9.90000,
-    ] * u.angstrom
-
-    assert u.allclose(idl_mirror_wavelength_manu, wavelength[80:90])
-
-
 @pytest.mark.parametrize("channel_name", channel_names)
 def test_channel_transmission(channel_name):
     channel_filter = Channel(channel_name)