Just testing, pushing foward to start a new branch

xrtpy/response/tests/test_effective_area.py

@@ -113,6 +113,52 @@ def get_IDL_data_files():
     return files
 
 
+# #Working testing
+# @pytest.mark.parametrize("filename", get_IDL_data_files())
+# def test_effective_area_compare_idl(filename):
+#     print(f"\n\nTesting file: {filename}\n")
+
+#     # Read the filter name and observation date from the file
+#     with filename.open() as f:
+#         filter_name = f.readline().split()[1]
+#         filter_obs_date = " ".join(f.readline().split()[1:])
+#         print(f"Filter name: {filter_name}, Observation date: {filter_obs_date}")  # Debugging output
+
+#     # Correct non-standard date format
+#     filter_obs_date = filter_obs_date.replace("Sept", "Sep")
+
+#     # Load IDL data from the file
+#     IDL_data = np.loadtxt(filename, skiprows=3)
+#     IDL_wavelength = IDL_data[:, 0] * u.AA
+#     IDL_effective_area = IDL_data[:, 1] * u.cm**2
+#     #print(f"Loaded IDL data: {IDL_data.shape}\n")
+
+#     # Compute effective area using XRTpy
+#     instance = EffectiveAreaFundamental(filter_name, filter_obs_date)
+#     actual_effective_area = instance.effective_area()
+#     #print(f"Calculated effective area (first 5): {actual_effective_area[:5]}\n")
+
+#     # Interpolate IDL effective area values to align with XRTpy wavelengths
+#     IDL_effective_area_interp = np.interp(instance.wavelength, IDL_wavelength, IDL_effective_area)
+#     #print(f"Interpolated IDL effective area (first 5): {IDL_effective_area_interp[:5]}\n")
+
+#     # Compare XRTpy and IDL effective areas using rtol=1e-4
+#     assert u.allclose(
+#         actual_effective_area,
+#         IDL_effective_area_interp,
+#         rtol=1e-5,
+#     ), f"Effective areas differ for filter {filter_name} on {filter_obs_date}"
+
+
+################################################################################################
+################################################################################################
+################################################################################################
+
+import matplotlib.pyplot as plt
+from pathlib import Path
+
+OUTPUT_DIR = Path("effective_area_plots")
+OUTPUT_DIR.mkdir(exist_ok=True)  # Create the main output directory if it doesn't exist
 
 @pytest.mark.parametrize("filename", get_IDL_data_files())
 def test_effective_area_compare_idl(filename):
@@ -123,33 +169,177 @@ def test_effective_area_compare_idl(filename):
         filter_name = f.readline().split()[1]
         filter_obs_date = " ".join(f.readline().split()[1:])
         print(f"Filter name: {filter_name}, Observation date: {filter_obs_date}")  # Debugging output
-
+    
     # Correct non-standard date format
     filter_obs_date = filter_obs_date.replace("Sept", "Sep")
 
     # Load IDL data from the file
     IDL_data = np.loadtxt(filename, skiprows=3)
     IDL_wavelength = IDL_data[:, 0] * u.AA
     IDL_effective_area = IDL_data[:, 1] * u.cm**2
-    #print(f"Loaded IDL data: {IDL_data.shape}\n")
 
     # Compute effective area using XRTpy
     instance = EffectiveAreaFundamental(filter_name, filter_obs_date)
     actual_effective_area = instance.effective_area()
-    #print(f"Calculated effective area (first 5): {actual_effective_area[:5]}\n")
 
-    # Interpolate IDL effective area values to align with XRTpy wavelengths
-    IDL_effective_area_interp = np.interp(instance.wavelength, IDL_wavelength, IDL_effective_area)
-    #print(f"Interpolated IDL effective area (first 5): {IDL_effective_area_interp[:5]}\n")
+    # Interpolate XRTpy effective area onto the IDL wavelength grid
+    XRTpy_effective_area_interp = np.interp(
+        IDL_wavelength.value,  # Target grid (IDL wavelengths)
+        instance.wavelength.value,  # Source grid (XRTpy wavelengths)
+        actual_effective_area.value  # Data to interpolate (XRTpy effective area)
+    )
+
+    # Make both arrays dimensionless for comparison
+    IDL_effective_area_unitless = IDL_effective_area.value
+    XRTpy_effective_area_unitless = XRTpy_effective_area_interp
+
+    # Compare effective areas using relative tolerance (unitless comparison)
+    rtol = 1e-4
+    differences = np.abs(XRTpy_effective_area_unitless - IDL_effective_area_unitless)
+    max_diff = np.max(differences)
+    failed_indices = np.where(differences > rtol * np.abs(IDL_effective_area_unitless))[0]
+
+    if failed_indices.size > 0:
+        print(f"Test failed for filter {filter_name} on {filter_obs_date}")
+        print(f"Max difference: {max_diff}")
+
+        # Save the plot
+        save_effective_area_plot(
+            IDL_wavelength.value,
+            IDL_effective_area_unitless,
+            XRTpy_effective_area_unitless,
+            failed_indices,
+            filter_name,
+            filter_obs_date,
+        )
+
+        # Fail the test
+        assert False, f"Effective areas differ for filter {filter_name} on {filter_obs_date}"
+
+
+def save_effective_area_plot(wavelength, IDL_area, XRTpy_area, failed_indices, filter_name, obs_date):
+    """
+    Saves the effective area comparison plot to a subfolder for the given filter.
+    """
+    # Create a subdirectory for the filter
+    filter_dir = OUTPUT_DIR / filter_name
+    filter_dir.mkdir(parents=True, exist_ok=True)
+
+    # Define the output filename
+    output_file = filter_dir / f"{filter_name}_{obs_date.replace(':', '').replace(' ', '_')}.png"
+
+    # Plot the data
+    plt.figure(figsize=(10, 6))
+    plt.plot(wavelength, IDL_area, label="IDL Effective Area", color="blue", lw=2)
+    plt.plot(wavelength, XRTpy_area, label="XRTpy Effective Area", color="green", linestyle="--", lw=2)
+    plt.scatter(
+        wavelength[failed_indices],
+        IDL_area[failed_indices],
+        color="red",
+        label="Failed Points",
+        zorder=5,
+    )
+    plt.xlabel("Wavelength (Å)", fontsize=14)
+    plt.ylabel("Effective Area (cm²)", fontsize=14)
+    plt.title(f"Effective Area Comparison for {filter_name} on {obs_date}", fontsize=16)
+    plt.legend(fontsize=12)
+    plt.grid(True, linestyle="--", alpha=0.7)
+    plt.tight_layout()
+
+    # Save the plot
+    plt.savefig(output_file, dpi=300)
+    plt.close()
+    print(f"Saved plot to {output_file}")
+
+
+# import matplotlib.pyplot as plt
+
+# @pytest.mark.parametrize("filename", get_IDL_data_files())
+# def test_effective_area_compare_idl(filename):
+#     print(f"\n\nTesting file: {filename}\n")
+
+#     # Read the filter name and observation date from the file
+#     with filename.open() as f:
+#         filter_name = f.readline().split()[1]
+#         filter_obs_date = " ".join(f.readline().split()[1:])
+#         print(f"Filter name: {filter_name}, Observation date: {filter_obs_date}")  # Debugging output
+
+#     # Correct non-standard date format
+#     filter_obs_date = filter_obs_date.replace("Sept", "Sep")
+
+#     # Load IDL data from the file
+#     IDL_data = np.loadtxt(filename, skiprows=3)
+#     IDL_wavelength = IDL_data[:, 0] * u.AA
+#     IDL_effective_area = IDL_data[:, 1] * u.cm**2
+
+#     # Compute effective area using XRTpy
+#     instance = EffectiveAreaFundamental(filter_name, filter_obs_date)
+#     actual_effective_area = instance.effective_area()
 
-    # Compare XRTpy and IDL effective areas using rtol=1e-4
-    assert u.allclose(
-        actual_effective_area,
-        IDL_effective_area_interp,
-        rtol=1e-4,
-    ), f"Effective areas differ for filter {filter_name} on {filter_obs_date}"
+#     # Interpolate XRTpy effective area onto the IDL wavelength grid
+#     XRTpy_effective_area_interp = np.interp(
+#         IDL_wavelength.value,  # Target grid (IDL wavelengths)
+#         instance.wavelength.value,  # Source grid (XRTpy wavelengths)
+#         actual_effective_area.value  # Data to interpolate (XRTpy effective area)
+#     )
+
+#     # Make both arrays dimensionless for comparison
+#     IDL_effective_area_unitless = IDL_effective_area.value
+#     XRTpy_effective_area_unitless = XRTpy_effective_area_interp
+
+#     # Compare effective areas using relative tolerance (unitless comparison)
+#     rtol = 1e-4
+#     differences = np.abs(XRTpy_effective_area_unitless - IDL_effective_area_unitless)
+#     max_diff = np.max(differences)
+#     failed_indices = np.where(differences > rtol * np.abs(IDL_effective_area_unitless))[0]
+
+#     if failed_indices.size > 0:
+#         print(f"Test failed for filter {filter_name} on {filter_obs_date}")
+#         print(f"Max difference: {max_diff}")
+
+#         # Plot the results
+#         plot_effective_area_comparison(
+#             IDL_wavelength.value,
+#             IDL_effective_area_unitless,
+#             XRTpy_effective_area_unitless,
+#             failed_indices,
+#             filter_name,
+#             filter_obs_date,
+#         )
+
+#         # Fail the test
+#         assert False, f"Effective areas differ for filter {filter_name} on {filter_obs_date}"
+
+
+# def plot_effective_area_comparison(wavelength, IDL_area, XRTpy_area, failed_indices, filter_name, obs_date):
+#     """
+#     Plots the effective area comparison between XRTpy and IDL.
+#     Marks failed points in red.
+#     """
+#     plt.figure(figsize=(10, 6))
+#     plt.plot(wavelength, IDL_area, label="IDL Effective Area", color="blue", lw=2)
+#     plt.plot(wavelength, XRTpy_area, label="XRTpy Effective Area", color="green", linestyle="--", lw=2)
+#     plt.scatter(
+#         wavelength[failed_indices],
+#         IDL_area[failed_indices],
+#         color="red",
+#         label="Failed Points",
+#         zorder=5,
+#     )
+#     plt.xlabel("Wavelength (Å)", fontsize=14)
+#     plt.ylabel("Effective Area (cm²)", fontsize=14)
+#     plt.title(f"Effective Area Comparison for {filter_name} on {obs_date}", fontsize=16)
+#     plt.legend(fontsize=12)
+#     plt.grid(True, linestyle="--", alpha=0.7)
+#     plt.tight_layout()
+#     plt.show()
+
 
 
+########################################################################################################################
+########################################################################################################################
+########################################################################################################################
+# Original 
 ## 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.