succesfully producing output for the phi functions though mismatching…

… with the reference. TO be checked

agnpy/photo_meson/tests/test_kernels.py

@@ -0,0 +1,52 @@
+# test the kernels
+import numpy as np
+from pathlib import Path
+import pytest
+from agnpy.photo_meson.kernels import PhiKernel
+from agnpy.utils.validation_utils import (
+    make_comparison_plot,
+    extract_columns_sample_file,
+    check_deviation,
+    clean_and_make_dir,
+)
+
+agnpy_dir = Path(__file__).parent.parent.parent  # go to the agnpy root
+# where to read sampled files
+data_dir = agnpy_dir / "data"
+# where to save figures, clean-up before making the new
+figures_dir = clean_and_make_dir(agnpy_dir, "crosschecks/figures/photo_meson")
+
+
+class TestKernels:
+    """Class to test the phi functions representing the integration kernels"""
+
+    @pytest.mark.parametrize("particle", ["gamma", "positron", "muon_neutrino", "muon_antineutrino"])
+    @pytest.mark.parametrize("eta", ["1.5", "30"])
+    def test_phi_interpolation(self, particle, eta):
+        """Test the interpolations against the plots presented in the [KelnerAharonian2008]_ paper."""
+        # read the reference from file
+        x_ref, phi_ref = extract_columns_sample_file(
+            f"{data_dir}/photo_meson/kelner_aharonian_2008/phi_values/phi_{particle}_eta_{eta}.csv"
+        )
+
+        phi = PhiKernel(particle)
+        eta = float(eta)
+        phi_agnpy = phi(eta, x_ref)
+
+        # comparison plot
+        make_comparison_plot(
+            x=x_ref,
+            y_comp=x_ref * phi_agnpy,
+            y_ref=x_ref * phi_ref,
+            comp_label="agnpy",
+            ref_label="Kelner and Aharonian (2008)",
+            fig_title=r"$\phi$" + f" {particle.replace('_', ' ')}",
+            fig_path=f"{figures_dir}/phi_comparison_particle_{particle}_eta_{eta}.png",
+            plot_type="custom",
+            x_label=r"$x = E_{\gamma} / E_{\rm p}$",
+            y_label=r"$x \phi(\eta, x)$",
+            y_range=None,
+            comparison_range=None
+        )
+        # requires that the SED points deviate less than 25% from the figure
+        assert check_deviation(x_ref, phi_agnpy, phi_ref, 0.25)

agnpy/utils/validation_utils.py

@@ -48,10 +48,10 @@ def clean_and_make_dir(main_dir, sub_dir=None):
     return _dir
 
 
-def extract_columns_sample_file(sample_file, x_unit, y_unit=None):
+def extract_columns_sample_file(sample_file, x_unit=None, y_unit=None):
     """Return two arrays of quantities from a sample file."""
     sample_table = np.loadtxt(sample_file, delimiter=",", comments="#")
-    x = sample_table[:, 0] * u.Unit(x_unit)
+    x = sample_table[:, 0] if y_unit is None else sample_table[:, 0] * u.Unit(x_unit)
     y = sample_table[:, 1] if y_unit is None else sample_table[:, 1] * u.Unit(y_unit)
     return x, y
 
@@ -73,14 +73,16 @@ def check_deviation(x, y_comp, y_ref, rtol, x_range=None):
 
 
 def make_comparison_plot(
-    nu,
+    x,
     y_comp,
     y_ref,
     comp_label,
     ref_label,
     fig_title,
     fig_path,
-    plot_type,
+    plot_type="custom",
+    x_label="",
+    y_label="",
     y_range=None,
     comparison_range=None
 ):
@@ -125,11 +127,9 @@ def make_comparison_plot(
         x_label = TAU_X_LABEL
         y_label = TAU_Y_LABEL
         deviation_label = TAU_DEVIATION_LABEL
-    else:
+    elif plot_type == "custom": 
         # set a custom y label, keep the x-axis in frequency
-        x_label = SED_X_LABEL
-        y_label = plot_type
-        deviation_label = f"({plot_type} agnpy / {plot_type} ref.) - 1"
+        deviation_label = f"(agnpy / reference) - 1"
     # make the plot
     fig, ax = plt.subplots(
         2,
@@ -140,9 +140,9 @@ def make_comparison_plot(
 
     # plot the SEDs or TAUs in the upper panel
     # plot the reference sed with a continuous line and agnpy sed with a dashed one
-    ax[0].loglog(nu, y_ref, marker=".", ls="-", color="k", lw=1.5, label=ref_label)
+    ax[0].loglog(x, y_ref, marker=".", ls="-", color="k", lw=1.5, label=ref_label)
     ax[0].loglog(
-        nu, y_comp, marker=".", ls="--", color="crimson", lw=1.5, label=comp_label
+        x, y_comp, marker=".", ls="--", color="crimson", lw=1.5, label=comp_label
     )
     ax[0].set_ylabel(y_label)
     ax[0].set_title(fig_title)
@@ -163,7 +163,7 @@ def make_comparison_plot(
     ax[1].axhline(-0.3, ls=":", color="darkgray")
     ax[1].set_ylim([-0.5, 0.5])
     ax[1].semilogx(
-        nu,
+        x,
         deviation,
         marker=".",
         ls="--",