diff --git a/scripts/mass_via_integral/gaussian.py b/scripts/mass_via_integral/gaussian.py index 4424956..aadb8c0 100644 --- a/scripts/mass_via_integral/gaussian.py +++ b/scripts/mass_via_integral/gaussian.py @@ -29,6 +29,49 @@ def deflection_func(u, y, x, npow, axis_ratio, sigma): ) +""" +__Deflections via Integral__ + +Computes deflection angles by numerically integrating the Gaussian deflection +function. Uses the profile object's coordinate transform methods. +""" + + +def deflections_2d_via_integral_from(mp, grid): + transformed = mp.transformed_to_reference_frame_grid_from(grid) + axis_ratio = mp.axis_ratio() + + def calculate_deflection_component(npow, index): + deflection_grid = np.array(axis_ratio * transformed.array[:, index]) + + for i in range(transformed.shape[0]): + deflection_grid[i] *= ( + mp.intensity + * mp.mass_to_light_ratio + * quad( + deflection_func, + a=0.0, + b=1.0, + args=( + transformed.array[i, 0], + transformed.array[i, 1], + npow, + axis_ratio, + mp.sigma / np.sqrt(axis_ratio), + ), + )[0] + ) + + return deflection_grid + + deflection_y = calculate_deflection_component(1.0, 0) + deflection_x = calculate_deflection_component(0.0, 1) + + return mp.rotated_grid_from_reference_frame_from( + np.multiply(1.0, np.vstack((deflection_y, deflection_x)).T) + ) + + """ __Config 1__ """ @@ -40,13 +83,13 @@ def deflection_func(u, y, x, npow, axis_ratio, sigma): sigma=3.0, mass_to_light_ratio=1.0, ) -deflections = mp.deflections_2d_via_analytic_from( - grid=ag.Grid2DIrregular([[1.0, 0.0]]) +integral = deflections_2d_via_integral_from( + mp, 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) +np.testing.assert_allclose(integral, analytic.array, rtol=1e-3) print("Config 1: PASSED") """ @@ -60,13 +103,13 @@ def deflection_func(u, y, x, npow, axis_ratio, sigma): sigma=5.0, mass_to_light_ratio=1.0, ) -deflections = mp.deflections_2d_via_analytic_from( - grid=ag.Grid2DIrregular([[0.5, 0.2]]) +integral = deflections_2d_via_integral_from( + mp, 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) +np.testing.assert_allclose(integral, analytic.array, rtol=1e-3) print("Config 2: PASSED") """ @@ -80,13 +123,13 @@ def deflection_func(u, y, x, npow, axis_ratio, sigma): sigma=5.0, mass_to_light_ratio=2.0, ) -deflections = mp.deflections_2d_via_analytic_from( - grid=ag.Grid2DIrregular([[0.5, 0.2]]) +integral = deflections_2d_via_integral_from( + mp, 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) +np.testing.assert_allclose(integral, analytic.array, rtol=1e-3) print("Mass to light 2: PASSED") """ @@ -100,13 +143,13 @@ def deflection_func(u, y, x, npow, axis_ratio, sigma): sigma=5.0, mass_to_light_ratio=1.0, ) -deflections = mp.deflections_2d_via_analytic_from( - grid=ag.Grid2DIrregular([[0.5, 0.2]]) +integral = deflections_2d_via_integral_from( + mp, 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) +np.testing.assert_allclose(integral, 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 index 68b76de..a893274 100644 --- a/scripts/mass_via_integral/gnfw.py +++ b/scripts/mass_via_integral/gnfw.py @@ -11,6 +11,286 @@ import autogalaxy as ag +epsrel = 1.49e-6 + + +""" +__Tabulate Integral__ + +Helper to set up the tabulation grid for the gNFW integral. +""" + + +def tabulate_integral(mp, grid, tabulate_bins): + eta_min = 1.0e-4 + eta_max = 1.05 * np.max(mp.elliptical_radii_grid_from(grid=grid)) + + minimum_log_eta = np.log10(eta_min) + maximum_log_eta = np.log10(eta_max) + bin_size = (maximum_log_eta - minimum_log_eta) / (tabulate_bins - 1) + + return eta_min, eta_max, minimum_log_eta, maximum_log_eta, bin_size + + +""" +__Surface Density Integrand (Elliptical gNFW)__ +""" + + +def surface_density_integrand(x, kappa_radius, scale_radius, inner_slope): + return ( + (3 - inner_slope) + * (x + kappa_radius / scale_radius) ** (inner_slope - 4) + * (1 - np.sqrt(1 - x * x)) + ) + + +""" +__gNFW Deflection Function (Elliptical, tabulated)__ +""" + + +def gnfw_deflection_func( + u, + y, + x, + npow, + axis_ratio, + minimum_log_eta, + maximum_log_eta, + tabulate_bins, + surface_density_integral, +): + _eta_u = np.sqrt((u * ((x**2) + (y**2 / (1 - (1 - axis_ratio**2) * u))))) + bin_size = (maximum_log_eta - minimum_log_eta) / (tabulate_bins - 1) + i = 1 + int((np.log10(_eta_u) - minimum_log_eta) / bin_size) + r1 = 10.0 ** (minimum_log_eta + (i - 1) * bin_size) + r2 = r1 * 10.0**bin_size + kap = surface_density_integral[i] + ( + surface_density_integral[i + 1] - surface_density_integral[i] + ) * (_eta_u - r1) / (r2 - r1) + return kap / (1.0 - (1.0 - axis_ratio**2) * u) ** (npow + 0.5) + + +""" +__gNFW Potential Function (Elliptical, tabulated)__ +""" + + +def deflection_integrand_for_potential(x, kappa_radius, scale_radius, inner_slope): + return (x + kappa_radius / scale_radius) ** (inner_slope - 3) * ( + (1 - np.sqrt(1 - x**2)) / x + ) + + +def gnfw_potential_func( + u, + y, + x, + axis_ratio, + minimum_log_eta, + maximum_log_eta, + tabulate_bins, + potential_integral, +): + _eta_u = np.sqrt((u * ((x**2) + (y**2 / (1 - (1 - axis_ratio**2) * u))))) + bin_size = (maximum_log_eta - minimum_log_eta) / (tabulate_bins - 1) + i = 1 + int((np.log10(_eta_u) - minimum_log_eta) / bin_size) + r1 = 10.0 ** (minimum_log_eta + (i - 1) * bin_size) + r2 = r1 * 10.0**bin_size + angle = potential_integral[i] + ( + potential_integral[i + 1] - potential_integral[i] + ) * (_eta_u - r1) / (r2 - r1) + return _eta_u * (angle / u) / (1.0 - (1.0 - axis_ratio**2) * u) ** 0.5 + + +""" +__gNFWSph Deflection Integrand and Function__ +""" + + +def deflection_integrand_sph(y, eta, inner_slope): + return (y + eta) ** (inner_slope - 3) * ((1 - np.sqrt(1 - y**2)) / y) + + +def deflection_func_sph(mp, eta): + integral_y_2 = quad( + deflection_integrand_sph, + a=0.0, + b=1.0, + args=(eta, mp.inner_slope), + epsrel=1.49e-6, + )[0] + return eta ** (2 - mp.inner_slope) * ( + (1.0 / (3 - mp.inner_slope)) + * special.hyp2f1( + 3 - mp.inner_slope, 3 - mp.inner_slope, 4 - mp.inner_slope, -eta + ) + + integral_y_2 + ) + + +""" +__Deflections via Integral (Elliptical gNFW)__ +""" + + +def deflections_2d_via_integral_from(mp, grid, tabulate_bins=1000): + transformed = mp.transformed_to_reference_frame_grid_from(grid) + axis_ratio = mp.axis_ratio() + + ( + eta_min, + eta_max, + minimum_log_eta, + maximum_log_eta, + bin_size, + ) = tabulate_integral(mp, transformed, tabulate_bins) + + surface_density_integral_arr = np.zeros((tabulate_bins,)) + + for i in range(tabulate_bins): + eta = 10.0 ** (minimum_log_eta + (i - 1) * bin_size) + + integral = quad( + surface_density_integrand, + a=0.0, + b=1.0, + args=(eta, mp.scale_radius, mp.inner_slope), + epsrel=epsrel, + )[0] + + surface_density_integral_arr[i] = ( + (eta / mp.scale_radius) ** (1 - mp.inner_slope) + ) * (((1 + eta / mp.scale_radius) ** (mp.inner_slope - 3)) + integral) + + def calculate_deflection_component(npow, yx_index): + deflection_grid = np.zeros(transformed.shape[0]) + + for i in range(transformed.shape[0]): + deflection_grid[i] = ( + 2.0 + * mp.kappa_s + * axis_ratio + * transformed.array[i, yx_index] + * quad( + gnfw_deflection_func, + a=0.0, + b=1.0, + args=( + transformed.array[i, 0], + transformed.array[i, 1], + npow, + axis_ratio, + minimum_log_eta, + maximum_log_eta, + tabulate_bins, + surface_density_integral_arr, + ), + epsrel=epsrel, + )[0] + ) + + return deflection_grid + + deflection_y = calculate_deflection_component(npow=1.0, yx_index=0) + deflection_x = calculate_deflection_component(npow=0.0, yx_index=1) + + return mp.rotated_grid_from_reference_frame_from( + np.multiply(1.0, np.vstack((deflection_y, deflection_x)).T), + ) + + +""" +__Deflections via Integral (Spherical gNFWSph)__ +""" + + +def deflections_2d_via_integral_sph_from(mp, grid): + transformed = mp.transformed_to_reference_frame_grid_from(grid) + + eta = np.multiply( + 1.0 / mp.scale_radius, + mp.radial_grid_from(transformed, is_transformed=True).array, + ) + + deflection_grid = np.zeros(transformed.shape[0]) + + for i in range(transformed.shape[0]): + deflection_grid[i] = np.multiply( + 4.0 * mp.kappa_s * mp.scale_radius, deflection_func_sph(mp, eta[i]) + ) + + return mp._cartesian_grid_via_radial_from( + grid=transformed, radius=deflection_grid, xp=np + ) + + +""" +__Potential via Integral (Elliptical gNFW)__ +""" + + +def potential_2d_via_integral_from(mp, grid, tabulate_bins=1000): + transformed = mp.transformed_to_reference_frame_grid_from(grid) + axis_ratio = mp.axis_ratio() + + ( + eta_min, + eta_max, + minimum_log_eta, + maximum_log_eta, + bin_size, + ) = tabulate_integral(mp, transformed, tabulate_bins) + + deflection_integral = np.zeros((tabulate_bins,)) + + for i in range(tabulate_bins): + eta = 10.0 ** (minimum_log_eta + (i - 1) * bin_size) + + integral = quad( + deflection_integrand_for_potential, + a=0.0, + b=1.0, + args=(eta, mp.scale_radius, mp.inner_slope), + epsrel=epsrel, + )[0] + + deflection_integral[i] = ( + (eta / mp.scale_radius) ** (2 - mp.inner_slope) + ) * ( + (1.0 / (3 - mp.inner_slope)) + * special.hyp2f1( + 3 - mp.inner_slope, + 3 - mp.inner_slope, + 4 - mp.inner_slope, + -(eta / mp.scale_radius), + ) + + integral + ) + + potential_grid = np.zeros(transformed.shape[0]) + + for i in range(transformed.shape[0]): + potential_grid[i] = (2.0 * mp.kappa_s * axis_ratio) * quad( + gnfw_potential_func, + a=0.0, + b=1.0, + args=( + transformed.array[i, 0], + transformed.array[i, 1], + axis_ratio, + minimum_log_eta, + maximum_log_eta, + tabulate_bins, + deflection_integral, + ), + epsrel=epsrel, + )[0] + + return potential_grid + + """ __gNFWSph Config 1__ """ @@ -18,11 +298,11 @@ 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_mge_from( - grid=ag.Grid2DIrregular([[0.1875, 0.1625]]) +integral = deflections_2d_via_integral_sph_from( + mp, 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) +mge = mp.deflections_2d_via_mge_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) +np.testing.assert_allclose(integral, mge.array, rtol=1e-3) print("gNFWSph config 1: PASSED") """ @@ -32,11 +312,11 @@ 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_mge_from( - grid=ag.Grid2DIrregular([[0.1875, 0.1625]]) +integral = deflections_2d_via_integral_sph_from( + mp, 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) +mge = mp.deflections_2d_via_mge_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) +np.testing.assert_allclose(integral, mge.array, rtol=1e-3) print("gNFWSph config 2: PASSED") """ @@ -50,11 +330,11 @@ inner_slope=0.5, scale_radius=8.0, ) -deflections = mp.deflections_2d_via_mge_from( - grid=ag.Grid2DIrregular([[0.1875, 0.1625]]) +integral = deflections_2d_via_integral_from( + mp, 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) +mge = mp.deflections_2d_via_mge_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) +np.testing.assert_allclose(integral, mge.array, rtol=1e-3) print("gNFW ell config 1: PASSED") """ @@ -68,11 +348,11 @@ inner_slope=1.5, scale_radius=4.0, ) -deflections = mp.deflections_2d_via_mge_from( - grid=ag.Grid2DIrregular([[0.1875, 0.1625]]) +integral = deflections_2d_via_integral_from( + mp, 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) +mge = mp.deflections_2d_via_mge_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) +np.testing.assert_allclose(integral, mge.array, rtol=1e-3) print("gNFW ell config 2: PASSED") """ @@ -82,8 +362,10 @@ 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) +integral_potential = potential_2d_via_integral_from( + mp, grid=ag.Grid2DIrregular([[0.1625, 0.1875]]) +) +np.testing.assert_allclose(integral_potential, 0.00920, rtol=1e-3) print("Potential inner_slope=0.5: PASSED") """ @@ -93,8 +375,10 @@ 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) +integral_potential = potential_2d_via_integral_from( + mp, grid=ag.Grid2DIrregular([[0.1625, 0.1875]]) +) +np.testing.assert_allclose(integral_potential, 0.17448, rtol=1e-3) print("Potential inner_slope=1.5: PASSED") """ @@ -108,8 +392,10 @@ 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) +integral_potential = potential_2d_via_integral_from( + mp, grid=ag.Grid2DIrregular([[2.0, 2.0]]) +) +np.testing.assert_allclose(integral_potential, 2.4718, rtol=1e-4) print("Potential gNFW ell: PASSED") """ diff --git a/scripts/mass_via_integral/gnfw_virial_mass_conc.py b/scripts/mass_via_integral/gnfw_virial_mass_conc.py index a7cb738..06c553f 100644 --- a/scripts/mass_via_integral/gnfw_virial_mass_conc.py +++ b/scripts/mass_via_integral/gnfw_virial_mass_conc.py @@ -7,9 +7,64 @@ source code. """ import numpy as np +from scipy import special +from scipy.integrate import quad import autogalaxy as ag +""" +__gNFWSph Deflection Integrand and Function__ + +Same spherical integral as gNFWSph — gNFWVirialMassConcSph inherits from it. +""" + + +def deflection_integrand_sph(y, eta, inner_slope): + return (y + eta) ** (inner_slope - 3) * ((1 - np.sqrt(1 - y**2)) / y) + + +def deflection_func_sph(mp, eta): + integral_y_2 = quad( + deflection_integrand_sph, + a=0.0, + b=1.0, + args=(eta, mp.inner_slope), + epsrel=1.49e-6, + )[0] + return eta ** (2 - mp.inner_slope) * ( + (1.0 / (3 - mp.inner_slope)) + * special.hyp2f1( + 3 - mp.inner_slope, 3 - mp.inner_slope, 4 - mp.inner_slope, -eta + ) + + integral_y_2 + ) + + +""" +__Deflections via Integral (Spherical)__ +""" + + +def deflections_2d_via_integral_sph_from(mp, grid): + transformed = mp.transformed_to_reference_frame_grid_from(grid) + + eta = np.multiply( + 1.0 / mp.scale_radius, + mp.radial_grid_from(transformed, is_transformed=True).array, + ) + + deflection_grid = np.zeros(transformed.shape[0]) + + for i in range(transformed.shape[0]): + deflection_grid[i] = np.multiply( + 4.0 * mp.kappa_s * mp.scale_radius, deflection_func_sph(mp, eta[i]) + ) + + return mp._cartesian_grid_via_radial_from( + grid=transformed, radius=deflection_grid, xp=np + ) + + """ __Config 1__ """ @@ -23,11 +78,11 @@ redshift_source=1.0, inner_slope=1.0, ) -deflections = mp.deflections_2d_via_mge_from( - grid=ag.Grid2DIrregular([[0.1875, 0.1625]]) +integral = deflections_2d_via_integral_sph_from( + mp, 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) +mge = mp.deflections_2d_via_mge_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) +np.testing.assert_allclose(integral, mge.array, 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 index e8016eb..8e22e4f 100644 --- a/scripts/mass_via_integral/nfw.py +++ b/scripts/mass_via_integral/nfw.py @@ -67,16 +67,86 @@ def nfw_potential_func(u, y, x, axis_ratio, kappa_s, scale_radius): ) +""" +__Deflections via Integral__ +""" + + +def deflections_2d_via_integral_from(mp, grid): + transformed = mp.transformed_to_reference_frame_grid_from(grid) + axis_ratio = mp.axis_ratio() + + def calculate_deflection_component(npow, index): + deflection_grid = np.array(axis_ratio * transformed.array[:, index]) + + for i in range(transformed.shape[0]): + deflection_grid[i] *= ( + mp.kappa_s + * quad( + nfw_deflection_func, + a=0.0, + b=1.0, + args=( + transformed.array[i, 0], + transformed.array[i, 1], + npow, + axis_ratio, + mp.scale_radius, + ), + )[0] + ) + + return deflection_grid + + deflection_y = calculate_deflection_component(1.0, 0) + deflection_x = calculate_deflection_component(0.0, 1) + + return mp.rotated_grid_from_reference_frame_from( + np.multiply(1.0, np.vstack((deflection_y, deflection_x)).T) + ) + + +""" +__Potential via Integral__ +""" + + +def potential_2d_via_integral_from(mp, grid): + transformed = mp.transformed_to_reference_frame_grid_from(grid) + axis_ratio = mp.axis_ratio() + + potential_grid = np.zeros(transformed.shape[0]) + + for i in range(transformed.shape[0]): + potential_grid[i] = quad( + nfw_potential_func, + a=0.0, + b=1.0, + args=( + transformed.array[i, 0], + transformed.array[i, 1], + axis_ratio, + mp.kappa_s, + mp.scale_radius, + ), + epsrel=1.49e-5, + )[0] + + return potential_grid + + """ __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_analytic_from( +integral = deflections_2d_via_integral_from( + nfw, grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +) +analytic = nfw.deflections_2d_via_analytic_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) +np.testing.assert_allclose(integral, analytic.array, rtol=1e-3) print("NFWSph config 1: PASSED") """ @@ -84,11 +154,13 @@ def nfw_potential_func(u, y, x, axis_ratio, kappa_s, scale_radius): """ nfw = ag.mp.NFWSph(centre=(0.3, 0.2), kappa_s=2.5, scale_radius=4.0) -deflections = nfw.deflections_2d_via_analytic_from( +integral = deflections_2d_via_integral_from( + nfw, grid=ag.Grid2DIrregular([[0.1875, 0.1625]]) +) +analytic = nfw.deflections_2d_via_analytic_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) +np.testing.assert_allclose(integral, analytic.array, rtol=1e-3) print("NFWSph config 2: PASSED") """ @@ -101,11 +173,13 @@ def nfw_potential_func(u, y, x, axis_ratio, kappa_s, scale_radius): kappa_s=1.0, scale_radius=1.0, ) -deflections = nfw.deflections_2d_via_analytic_from( +integral = deflections_2d_via_integral_from( + nfw, grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +) +analytic = nfw.deflections_2d_via_analytic_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) +np.testing.assert_allclose(integral, analytic.array, rtol=1e-3) print("NFW ell config 1: PASSED") """ @@ -118,11 +192,13 @@ def nfw_potential_func(u, y, x, axis_ratio, kappa_s, scale_radius): kappa_s=2.5, scale_radius=4.0, ) -deflections = nfw.deflections_2d_via_analytic_from( +integral = deflections_2d_via_integral_from( + nfw, grid=ag.Grid2DIrregular([(0.1625, 0.1625)]) +) +analytic = nfw.deflections_2d_via_analytic_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) +np.testing.assert_allclose(integral, analytic.array, 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 index cb135e5..71ca1af 100644 --- a/scripts/mass_via_integral/sersic.py +++ b/scripts/mass_via_integral/sersic.py @@ -27,6 +27,49 @@ def deflection_func( ) / ((1 - (1 - axis_ratio**2) * u) ** (npow + 0.5)) +""" +__Deflections via Integral__ +""" + + +def deflections_2d_via_integral_from(mp, grid): + transformed = mp.transformed_to_reference_frame_grid_from(grid) + axis_ratio = mp.axis_ratio() + sersic_constant = mp.sersic_constant + + def calculate_deflection_component(npow, index): + deflection_grid = axis_ratio * transformed.array[:, index] + + for i in range(transformed.shape[0]): + deflection_grid[i] = deflection_grid[i] * ( + mp.intensity + * mp.mass_to_light_ratio + * quad( + deflection_func, + a=0.0, + b=1.0, + args=( + transformed.array[i, 0], + transformed.array[i, 1], + npow, + axis_ratio, + mp.sersic_index, + mp.effective_radius, + sersic_constant, + ), + )[0] + ) + + return deflection_grid + + deflection_y = calculate_deflection_component(1.0, 0) + deflection_x = calculate_deflection_component(0.0, 1) + + return mp.rotated_grid_from_reference_frame_from( + np.multiply(1.0, np.vstack((deflection_y, deflection_x)).T) + ) + + """ __Config 1__ """ @@ -39,11 +82,11 @@ def deflection_func( sersic_index=2.0, mass_to_light_ratio=1.0, ) -deflections = mp.deflections_2d_via_cse_from( - grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +integral = deflections_2d_via_integral_from( + mp, 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) +cse = mp.deflections_2d_via_cse_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) +np.testing.assert_allclose(integral, cse.array, rtol=1e-3) print("Config 1: PASSED") """ @@ -58,11 +101,11 @@ def deflection_func( sersic_index=3.0, mass_to_light_ratio=1.0, ) -deflections = mp.deflections_2d_via_cse_from( - grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +integral = deflections_2d_via_integral_from( + mp, 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) +cse = mp.deflections_2d_via_cse_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) +np.testing.assert_allclose(integral, cse.array, 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 index c8f567c..8238416 100644 --- a/scripts/mass_via_integral/sersic_gradient.py +++ b/scripts/mass_via_integral/sersic_gradient.py @@ -39,6 +39,49 @@ def deflection_func( ) +""" +__Deflections via Integral__ +""" + + +def deflections_2d_via_integral_from(mp, grid): + transformed = mp.transformed_to_reference_frame_grid_from(grid) + axis_ratio = mp.axis_ratio() + sersic_constant = mp.sersic_constant + + def calculate_deflection_component(npow, index): + deflection_grid = np.array(axis_ratio * transformed.array[:, index]) + + for i in range(transformed.shape[0]): + deflection_grid[i] *= ( + mp.intensity + * mp.mass_to_light_ratio + * quad( + deflection_func, + a=0.0, + b=1.0, + args=( + transformed.array[i, 0], + transformed.array[i, 1], + npow, + axis_ratio, + mp.sersic_index, + mp.effective_radius, + mp.mass_to_light_gradient, + sersic_constant, + ), + )[0] + ) + return deflection_grid + + deflection_y = calculate_deflection_component(1.0, 0) + deflection_x = calculate_deflection_component(0.0, 1) + + return mp.rotated_grid_from_reference_frame_from( + np.multiply(1.0, np.vstack((deflection_y, deflection_x)).T) + ) + + """ __Gradient 1__ """ @@ -52,11 +95,11 @@ def deflection_func( mass_to_light_ratio=1.0, mass_to_light_gradient=1.0, ) -deflections = mp.deflections_2d_via_cse_from( - grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +integral = deflections_2d_via_integral_from( + mp, 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) +cse = mp.deflections_2d_via_cse_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) +np.testing.assert_allclose(integral, cse.array, rtol=1e-3) print("Gradient 1: PASSED") """ @@ -72,11 +115,11 @@ def deflection_func( mass_to_light_ratio=1.0, mass_to_light_gradient=-1.0, ) -deflections = mp.deflections_2d_via_cse_from( - grid=ag.Grid2DIrregular([[0.1625, 0.1625]]) +integral = deflections_2d_via_integral_from( + mp, 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) +cse = mp.deflections_2d_via_cse_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) +np.testing.assert_allclose(integral, cse.array, rtol=1e-3) print("Gradient -1: PASSED") print("\nAll SersicGradient integral tests passed.")