diff --git a/scripts/mass_via_integral/__init__.py b/scripts/mass_via_integral/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/scripts/mass_via_integral/gaussian.py b/scripts/mass_via_integral/gaussian.py new file mode 100644 index 0000000..e010c8e --- /dev/null +++ b/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.") diff --git a/scripts/mass_via_integral/gnfw.py b/scripts/mass_via_integral/gnfw.py new file mode 100644 index 0000000..f51a7d4 --- /dev/null +++ b/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.") diff --git a/scripts/mass_via_integral/gnfw_virial_mass_conc.py b/scripts/mass_via_integral/gnfw_virial_mass_conc.py new file mode 100644 index 0000000..08c1a19 --- /dev/null +++ b/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.") diff --git a/scripts/mass_via_integral/nfw.py b/scripts/mass_via_integral/nfw.py new file mode 100644 index 0000000..1d22d91 --- /dev/null +++ b/scripts/mass_via_integral/nfw.py @@ -0,0 +1,128 @@ +""" +NFW Mass Profile — Deflections and Potential via Integral +========================================================= + +This script preserves the integral-based deflection and potential calculations +for the NFW mass profile, which were removed from the autogalaxy source code. +""" +import numpy as np +from scipy.integrate import quad +import autogalaxy as ag + + +""" +__NFW Deflection Function__ +""" + + +def nfw_deflection_func(u, y, x, npow, axis_ratio, scale_radius): + _eta_u = (1.0 / scale_radius) * np.sqrt( + (u * ((x**2) + (y**2 / (1 - (1 - axis_ratio**2) * u)))) + ) + if _eta_u > 1: + _eta_u_2 = (1.0 / np.sqrt(_eta_u**2 - 1)) * np.arctan( + np.sqrt(_eta_u**2 - 1) + ) + elif _eta_u < 1: + _eta_u_2 = (1.0 / np.sqrt(1 - _eta_u**2)) * np.arctanh( + np.sqrt(1 - _eta_u**2) + ) + else: + _eta_u_2 = 1 + return ( + 2.0 + * (1 - _eta_u_2) + / (_eta_u**2 - 1) + / ((1 - (1 - axis_ratio**2) * u) ** (npow + 0.5)) + ) + + +""" +__NFW Potential Function__ +""" + + +def nfw_potential_func(u, y, x, axis_ratio, kappa_s, scale_radius): + _eta_u = (1.0 / scale_radius) * np.sqrt( + (u * ((x**2) + (y**2 / (1 - (1 - axis_ratio**2) * u)))) + ) + if _eta_u > 1: + _eta_u_2 = (1.0 / np.sqrt(_eta_u**2 - 1)) * np.arctan( + np.sqrt(_eta_u**2 - 1) + ) + elif _eta_u < 1: + _eta_u_2 = (1.0 / np.sqrt(1 - _eta_u**2)) * np.arctanh( + np.sqrt(1 - _eta_u**2) + ) + else: + _eta_u_2 = 1 + return ( + 4.0 + * kappa_s + * scale_radius + * (axis_ratio / 2.0) + * (_eta_u / u) + * ((np.log(_eta_u / 2.0) + _eta_u_2) / _eta_u) + / ((1 - (1 - axis_ratio**2) * u) ** 0.5) + ) + + +""" +__NFWSph Config 1__ +""" + +nfw = ag.mp.NFWSph(centre=(0.0, 0.0), kappa_s=1.0, scale_radius=1.0) +deflections = nfw.deflections_2d_via_integral_from( + grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +) +np.testing.assert_allclose(deflections[0, 0], 0.56194, rtol=1e-3) +np.testing.assert_allclose(deflections[0, 1], 0.56194, rtol=1e-3) +print("NFWSph config 1: PASSED") + +""" +__NFWSph Config 2__ +""" + +nfw = ag.mp.NFWSph(centre=(0.3, 0.2), kappa_s=2.5, scale_radius=4.0) +deflections = nfw.deflections_2d_via_integral_from( + grid=ag.Grid2DIrregular([[0.1875, 0.1625]]) +) +np.testing.assert_allclose(deflections[0, 0], -2.08909, rtol=1e-3) +np.testing.assert_allclose(deflections[0, 1], -0.69636, rtol=1e-3) +print("NFWSph config 2: PASSED") + +""" +__NFW Elliptical Config 1__ +""" + +nfw = ag.mp.NFW( + centre=(0.0, 0.0), + ell_comps=(0.0, 0.0), + kappa_s=1.0, + scale_radius=1.0, +) +deflections = nfw.deflections_2d_via_integral_from( + grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +) +np.testing.assert_allclose(deflections[0, 0], 0.56194, rtol=1e-3) +np.testing.assert_allclose(deflections[0, 1], 0.56194, rtol=1e-3) +print("NFW ell config 1: PASSED") + +""" +__NFW Elliptical Config 2__ +""" + +nfw = ag.mp.NFW( + centre=(0.3, 0.2), + ell_comps=(0.03669, 0.172614), + kappa_s=2.5, + scale_radius=4.0, +) +deflections = nfw.deflections_2d_via_integral_from( + grid=ag.Grid2DIrregular([(0.1625, 0.1625)]) +) +np.testing.assert_allclose(deflections[0, 0], -2.59480, rtol=1e-3) +np.testing.assert_allclose(deflections[0, 1], -0.44204, rtol=1e-3) +print("NFW ell config 2: PASSED") + +print("\nAll NFW integral tests passed.") diff --git a/scripts/mass_via_integral/sersic.py b/scripts/mass_via_integral/sersic.py new file mode 100644 index 0000000..980f235 --- /dev/null +++ b/scripts/mass_via_integral/sersic.py @@ -0,0 +1,68 @@ +""" +Sersic Mass Profile — Deflections via Integral +=============================================== + +This script preserves the integral-based deflection calculation for the +Sersic mass profile, which was removed from the autogalaxy source code. +""" +import numpy as np +from scipy.integrate import quad +import autogalaxy as ag + + +""" +__Deflection Function__ +""" + + +def deflection_func( + u, y, x, npow, axis_ratio, sersic_index, effective_radius, sersic_constant +): + _eta_u = np.sqrt(axis_ratio) * np.sqrt( + (u * ((x**2) + (y**2 / (1 - (1 - axis_ratio**2) * u)))) + ) + return np.exp( + -sersic_constant + * (((_eta_u / effective_radius) ** (1.0 / sersic_index)) - 1) + ) / ((1 - (1 - axis_ratio**2) * u) ** (npow + 0.5)) + + +""" +__Config 1__ +""" + +mp = ag.mp.Sersic( + centre=(-0.4, -0.2), + ell_comps=(-0.07142, -0.085116), + intensity=5.0, + effective_radius=0.2, + sersic_index=2.0, + mass_to_light_ratio=1.0, +) +deflections = mp.deflections_2d_via_integral_from( + grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +) +np.testing.assert_allclose(deflections[0, 0], 1.1446, rtol=1e-3) +np.testing.assert_allclose(deflections[0, 1], 0.79374, rtol=1e-3) +print("Config 1: PASSED") + +""" +__Config 2__ +""" + +mp = ag.mp.Sersic( + centre=(-0.4, -0.2), + ell_comps=(-0.07142, -0.085116), + intensity=10.0, + effective_radius=0.2, + sersic_index=3.0, + mass_to_light_ratio=1.0, +) +deflections = mp.deflections_2d_via_integral_from( + grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +) +np.testing.assert_allclose(deflections[0, 0], 2.6134, rtol=1e-3) +np.testing.assert_allclose(deflections[0, 1], 1.80719, rtol=1e-3) +print("Config 2: PASSED") + +print("\nAll Sersic integral tests passed.") diff --git a/scripts/mass_via_integral/sersic_gradient.py b/scripts/mass_via_integral/sersic_gradient.py new file mode 100644 index 0000000..c7c0f0f --- /dev/null +++ b/scripts/mass_via_integral/sersic_gradient.py @@ -0,0 +1,82 @@ +""" +Sersic Gradient Mass Profile — Deflections via Integral +======================================================= + +This script preserves the integral-based deflection calculation for the +SersicGradient mass profile, which was removed from the autogalaxy source code. +""" +import numpy as np +from scipy.integrate import quad +import autogalaxy as ag + + +""" +__Deflection Function__ +""" + + +def deflection_func( + u, + y, + x, + npow, + axis_ratio, + sersic_index, + effective_radius, + mass_to_light_gradient, + sersic_constant, +): + _eta_u = np.sqrt(axis_ratio) * np.sqrt( + (u * ((x**2) + (y**2 / (1 - (1 - axis_ratio**2) * u)))) + ) + return ( + (((axis_ratio * _eta_u) / effective_radius) ** -mass_to_light_gradient) + * np.exp( + -sersic_constant + * (((_eta_u / effective_radius) ** (1.0 / sersic_index)) - 1) + ) + / ((1 - (1 - axis_ratio**2) * u) ** (npow + 0.5)) + ) + + +""" +__Gradient 1__ +""" + +mp = ag.mp.SersicGradient( + centre=(-0.4, -0.2), + ell_comps=(-0.07142, -0.085116), + intensity=5.0, + effective_radius=0.2, + sersic_index=2.0, + mass_to_light_ratio=1.0, + mass_to_light_gradient=1.0, +) +deflections = mp.deflections_2d_via_integral_from( + grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +) +np.testing.assert_allclose(deflections[0, 0], 3.60324873535244, rtol=1e-3) +np.testing.assert_allclose(deflections[0, 1], 2.3638898009652, rtol=1e-3) +print("Gradient 1: PASSED") + +""" +__Gradient -1__ +""" + +mp = ag.mp.SersicGradient( + centre=(-0.4, -0.2), + ell_comps=(-0.07142, -0.085116), + intensity=5.0, + effective_radius=0.2, + sersic_index=2.0, + mass_to_light_ratio=1.0, + mass_to_light_gradient=-1.0, +) +deflections = mp.deflections_2d_via_integral_from( + grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +) +np.testing.assert_allclose(deflections[0, 0], 0.97806399756448, rtol=1e-3) +np.testing.assert_allclose(deflections[0, 1], 0.725459334118341, rtol=1e-3) +print("Gradient -1: PASSED") + +print("\nAll SersicGradient integral tests passed.")