feat: add mass_via_integral reference scripts

scripts/mass_via_integral/gaussian.py

@@ -0,0 +1,112 @@
+"""
+Gaussian Mass Profile — Deflections via Integral
+=================================================
+
+This script preserves the integral-based deflection calculation for the
+Gaussian mass profile, which was removed from the autogalaxy source code.
+
+The values asserted here were originally verified by the autogalaxy unit tests.
+"""
+import numpy as np
+from scipy.integrate import quad
+import autogalaxy as ag
+
+
+"""
+__Deflection Function__
+
+The integrand used by scipy.integrate.quad to compute deflection angles
+for the Gaussian mass profile.
+"""
+
+
+def deflection_func(u, y, x, npow, axis_ratio, sigma):
+    _eta_u = np.sqrt(axis_ratio) * np.sqrt(
+        (u * ((x**2) + (y**2 / (1 - (1 - axis_ratio**2) * u))))
+    )
+    return np.exp(-0.5 * np.square(np.divide(_eta_u, sigma))) / (
+        (1 - (1 - axis_ratio**2) * u) ** (npow + 0.5)
+    )
+
+
+"""
+__Config 1__
+"""
+
+mp = ag.mp.Gaussian(
+    centre=(0.0, 0.0),
+    ell_comps=(0.0, 0.05263),
+    intensity=1.0,
+    sigma=3.0,
+    mass_to_light_ratio=1.0,
+)
+deflections = mp.deflections_2d_via_integral_from(
+    grid=ag.Grid2DIrregular([[1.0, 0.0]])
+)
+analytic = mp.deflections_2d_via_analytic_from(
+    grid=ag.Grid2DIrregular([[1.0, 0.0]])
+)
+np.testing.assert_allclose(deflections, analytic.array, rtol=1e-3)
+print("Config 1: PASSED")
+
+"""
+__Config 2__
+"""
+
+mp = ag.mp.Gaussian(
+    centre=(0.0, 0.0),
+    ell_comps=(0.0, 0.111111),
+    intensity=1.0,
+    sigma=5.0,
+    mass_to_light_ratio=1.0,
+)
+deflections = mp.deflections_2d_via_integral_from(
+    grid=ag.Grid2DIrregular([[0.5, 0.2]])
+)
+analytic = mp.deflections_2d_via_analytic_from(
+    grid=ag.Grid2DIrregular([[0.5, 0.2]])
+)
+np.testing.assert_allclose(deflections, analytic.array, rtol=1e-3)
+print("Config 2: PASSED")
+
+"""
+__Mass to Light Ratio 2__
+"""
+
+mp = ag.mp.Gaussian(
+    centre=(0.0, 0.0),
+    ell_comps=(0.0, 0.111111),
+    intensity=1.0,
+    sigma=5.0,
+    mass_to_light_ratio=2.0,
+)
+deflections = mp.deflections_2d_via_integral_from(
+    grid=ag.Grid2DIrregular([[0.5, 0.2]])
+)
+analytic = mp.deflections_2d_via_analytic_from(
+    grid=ag.Grid2DIrregular([[0.5, 0.2]])
+)
+np.testing.assert_allclose(deflections, analytic.array, rtol=1e-3)
+print("Mass to light 2: PASSED")
+
+"""
+__Intensity 2__
+"""
+
+mp = ag.mp.Gaussian(
+    centre=(0.0, 0.0),
+    ell_comps=(0.0, 0.111111),
+    intensity=2.0,
+    sigma=5.0,
+    mass_to_light_ratio=1.0,
+)
+deflections = mp.deflections_2d_via_integral_from(
+    grid=ag.Grid2DIrregular([[0.5, 0.2]])
+)
+analytic = mp.deflections_2d_via_analytic_from(
+    grid=ag.Grid2DIrregular([[0.5, 0.2]])
+)
+np.testing.assert_allclose(deflections, analytic.array, rtol=1e-3)
+print("Intensity 2: PASSED")
+
+print("\nAll Gaussian integral tests passed.")

scripts/mass_via_integral/gnfw.py

@@ -0,0 +1,126 @@
+"""
+gNFW Mass Profile — Deflections, Potential, and Convergence via Integral
+========================================================================
+
+This script preserves the integral-based calculations for the gNFW mass profile,
+which were removed from the autogalaxy source code.
+"""
+import numpy as np
+from scipy import special
+from scipy.integrate import quad
+import autogalaxy as ag
+
+
+"""
+__gNFWSph Config 1__
+"""
+
+mp = ag.mp.gNFWSph(
+    centre=(0.0, 0.0), kappa_s=1.0, inner_slope=0.5, scale_radius=8.0
+)
+deflections = mp.deflections_2d_via_integral_from(
+    grid=ag.Grid2DIrregular([[0.1875, 0.1625]])
+)
+np.testing.assert_allclose(deflections[0, 0], 0.43501, rtol=1e-3)
+np.testing.assert_allclose(deflections[0, 1], 0.37701, rtol=1e-3)
+print("gNFWSph config 1: PASSED")
+
+"""
+__gNFWSph Config 2__
+"""
+
+mp = ag.mp.gNFWSph(
+    centre=(0.3, 0.2), kappa_s=2.5, inner_slope=1.5, scale_radius=4.0
+)
+deflections = mp.deflections_2d_via_integral_from(
+    grid=ag.Grid2DIrregular([[0.1875, 0.1625]])
+)
+np.testing.assert_allclose(deflections[0, 0], -9.31254, rtol=1e-3)
+np.testing.assert_allclose(deflections[0, 1], -3.10418, rtol=1e-3)
+print("gNFWSph config 2: PASSED")
+
+"""
+__gNFW Elliptical Config 1__
+"""
+
+mp = ag.mp.gNFW(
+    centre=(0.0, 0.0),
+    kappa_s=1.0,
+    ell_comps=ag.convert.ell_comps_from(axis_ratio=0.3, angle=100.0),
+    inner_slope=0.5,
+    scale_radius=8.0,
+)
+deflections = mp.deflections_2d_via_integral_from(
+    grid=ag.Grid2DIrregular([[0.1875, 0.1625]])
+)
+np.testing.assert_allclose(deflections[0, 0], 0.26604, rtol=1e-3)
+np.testing.assert_allclose(deflections[0, 1], 0.58988, rtol=1e-3)
+print("gNFW ell config 1: PASSED")
+
+"""
+__gNFW Elliptical Config 2__
+"""
+
+mp = ag.mp.gNFW(
+    centre=(0.3, 0.2),
+    kappa_s=2.5,
+    ell_comps=ag.convert.ell_comps_from(axis_ratio=0.5, angle=100.0),
+    inner_slope=1.5,
+    scale_radius=4.0,
+)
+deflections = mp.deflections_2d_via_integral_from(
+    grid=ag.Grid2DIrregular([[0.1875, 0.1625]])
+)
+np.testing.assert_allclose(deflections[0, 0], -5.99032, rtol=1e-3)
+np.testing.assert_allclose(deflections[0, 1], -4.02541, rtol=1e-3)
+print("gNFW ell config 2: PASSED")
+
+"""
+__Potential — gNFWSph Inner Slope 0.5__
+"""
+
+mp = ag.mp.gNFWSph(
+    centre=(0.0, 0.0), kappa_s=1.0, inner_slope=0.5, scale_radius=8.0
+)
+potential = mp.potential_2d_from(grid=ag.Grid2DIrregular([[0.1625, 0.1875]]))
+np.testing.assert_allclose(potential, 0.00920, rtol=1e-3)
+print("Potential inner_slope=0.5: PASSED")
+
+"""
+__Potential — gNFWSph Inner Slope 1.5__
+"""
+
+mp = ag.mp.gNFWSph(
+    centre=(0.0, 0.0), kappa_s=1.0, inner_slope=1.5, scale_radius=8.0
+)
+potential = mp.potential_2d_from(grid=ag.Grid2DIrregular([[0.1625, 0.1875]]))
+np.testing.assert_allclose(potential, 0.17448, rtol=1e-3)
+print("Potential inner_slope=1.5: PASSED")
+
+"""
+__Potential — gNFW Elliptical__
+"""
+
+mp = ag.mp.gNFW(
+    centre=(1.0, 1.0),
+    kappa_s=5.0,
+    ell_comps=ag.convert.ell_comps_from(axis_ratio=0.5, angle=100.0),
+    inner_slope=1.0,
+    scale_radius=10.0,
+)
+potential = mp.potential_2d_from(grid=ag.Grid2DIrregular([[2.0, 2.0]]))
+np.testing.assert_allclose(potential, 2.4718, rtol=1e-4)
+print("Potential gNFW ell: PASSED")
+
+"""
+__Convergence — gNFWSph__
+"""
+
+mp = ag.mp.gNFWSph(
+    centre=(0.0, 0.0), kappa_s=1.0, inner_slope=1.5, scale_radius=1.0
+)
+convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[2.0, 0.0]]))
+np.testing.assert_allclose(convergence, 0.30840, rtol=1e-2)
+print("Convergence gNFWSph: PASSED")
+
+print("\nAll gNFW integral tests passed.")

scripts/mass_via_integral/gnfw_virial_mass_conc.py

@@ -0,0 +1,33 @@
+"""
+gNFW Virial Mass Concentration — Deflections via Integral
+=========================================================
+
+This script preserves the integral-based deflection calculation for the
+gNFWVirialMassConcSph mass profile, which was removed from the autogalaxy
+source code.
+"""
+import numpy as np
+import autogalaxy as ag
+
+
+"""
+__Config 1__
+"""
+
+mp = ag.mp.gNFWVirialMassConcSph(
+    centre=(0.0, 0.0),
+    log10m_vir=12.0,
+    c_2=10.0,
+    overdens=0.0,
+    redshift_object=0.5,
+    redshift_source=1.0,
+    inner_slope=1.0,
+)
+deflections = mp.deflections_2d_via_integral_from(
+    grid=ag.Grid2DIrregular([[0.1875, 0.1625]])
+)
+np.testing.assert_allclose(deflections[0, 0], 0.0466231, rtol=1e-3)
+np.testing.assert_allclose(deflections[0, 1], 0.04040671, rtol=1e-3)
+print("gNFWVirialMassConcSph config 1: PASSED")
+
+print("\nAll gNFWVirialMassConcSph integral tests passed.")