Skip to content

refactor: remove integral-based deflection/potential methods from mass profiles #324

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Jammy2211 merged 1 commit into from
Apr 6, 2026
Merged

refactor: remove integral-based deflection/potential methods from mass profiles #324

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 0 additions & 20 deletions autogalaxy/profiles/mass/dark/abstract.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,26 +67,6 @@ def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):

return self.convergence_func(grid_radius=grid_eta, xp=xp)

def tabulate_integral(self, grid, tabulate_bins, **kwargs):
"""Tabulate an integral over the convergence of deflection potential of a mass profile. This is used in \
the gNFW profile classes to speed up the integration procedure.

Parameters
----------
grid
The grid of (y,x) arc-second coordinates the potential / deflection_stacks are computed on.
tabulate_bins
The number of bins to tabulate the inner integral of this profile.
"""
eta_min = 1.0e-4
eta_max = 1.05 * np.max(self.elliptical_radii_grid_from(grid=grid, **kwargs))

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

def density_3d_func(self, r, xp=np):
x = r.array / self.scale_radius

Expand Down
267 changes: 1 addition & 266 deletions autogalaxy/profiles/mass/dark/gnfw.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,118 +28,6 @@ def deflections_2d_via_mge_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs)
)
return deflections_via_mge

@aa.grid_dec.to_vector_yx
@aa.grid_dec.transform
def deflections_2d_via_integral_from(
self, grid: aa.type.Grid2DLike, xp=np, tabulate_bins=1000, **kwargs
):
"""
Calculate the deflection angles at a given set of arc-second gridded coordinates.

Parameters
----------
grid
The grid of (y,x) arc-second coordinates the deflection angles are computed on.
tabulate_bins
The number of bins to tabulate the inner integral of this profile.

"""

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))
)

def calculate_deflection_component(npow, yx_index):

from scipy.integrate import quad

deflection_grid = np.zeros(grid.shape[0])

for i in range(grid.shape[0]):
deflection_grid[i] = (
2.0
* self.kappa_s
* self.axis_ratio(xp)
* grid[i, yx_index]
* quad(
self.deflection_func,
a=0.0,
b=1.0,
args=(
grid.array[i, 0],
grid.array[i, 1],
npow,
self.axis_ratio(xp),
minimum_log_eta,
maximum_log_eta,
tabulate_bins,
surface_density_integral,
),
epsrel=gNFW.epsrel,
)[0]
)

return deflection_grid

(
eta_min,
eta_max,
minimum_log_eta,
maximum_log_eta,
bin_size,
) = self.tabulate_integral(grid, tabulate_bins)

surface_density_integral = np.zeros((tabulate_bins,))

for i in range(tabulate_bins):
from scipy.integrate import quad

eta = 10.0 ** (minimum_log_eta + (i - 1) * bin_size)

integral = quad(
surface_density_integrand,
a=0.0,
b=1.0,
args=(eta, self.scale_radius, self.inner_slope),
epsrel=gNFW.epsrel,
)[0]

surface_density_integral[i] = (
(eta / self.scale_radius) ** (1 - self.inner_slope)
) * (((1 + eta / self.scale_radius) ** (self.inner_slope - 3)) + integral)

deflection_y = calculate_deflection_component(npow=1.0, yx_index=0)
deflection_x = calculate_deflection_component(npow=0.0, yx_index=1)

return self.rotated_grid_from_reference_frame_from(
np.multiply(1.0, np.vstack((deflection_y, deflection_x)).T),
)

@staticmethod
def 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)

def convergence_func(self, grid_radius: float, xp=np) -> float:

from scipy.integrate import quad
Expand Down Expand Up @@ -170,108 +58,6 @@ def integral_y(y, eta):

return grid_radius

@aa.over_sample
@aa.grid_dec.to_array
@aa.grid_dec.transform
def potential_2d_from(
self, grid: aa.type.Grid2DLike, xp=np, tabulate_bins=1000, **kwargs
):
"""
Calculate the potential at a given set of arc-second gridded coordinates.

Parameters
----------
grid
The grid of (y,x) arc-second coordinates the deflection angles are computed on.
tabulate_bins
The number of bins to tabulate the inner integral of this profile.

"""

from scipy import special
from scipy.integrate import quad

def deflection_integrand(x, kappa_radius, scale_radius, inner_slope):
return (x + kappa_radius / scale_radius) ** (inner_slope - 3) * (
(1 - np.sqrt(1 - x**2)) / x
)

(
eta_min,
eta_max,
minimum_log_eta,
maximum_log_eta,
bin_size,
) = self.tabulate_integral(grid, tabulate_bins)

potential_grid = np.zeros(grid.shape[0])

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,
a=0.0,
b=1.0,
args=(eta, self.scale_radius, self.inner_slope),
epsrel=gNFW.epsrel,
)[0]

deflection_integral[i] = (
(eta / self.scale_radius) ** (2 - self.inner_slope)
) * (
(1.0 / (3 - self.inner_slope))
* special.hyp2f1(
3 - self.inner_slope,
3 - self.inner_slope,
4 - self.inner_slope,
-(eta / self.scale_radius),
)
+ integral
)

for i in range(grid.shape[0]):
potential_grid[i] = (2.0 * self.kappa_s * self.axis_ratio(xp)) * quad(
self.potential_func,
a=0.0,
b=1.0,
args=(
grid.array[i, 0],
grid.array[i, 1],
self.axis_ratio(xp),
minimum_log_eta,
maximum_log_eta,
tabulate_bins,
deflection_integral,
),
epsrel=gNFW.epsrel,
)[0]

return potential_grid

@staticmethod
def 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


class gNFWSph(gNFW):
def __init__(
Expand Down Expand Up @@ -306,55 +92,4 @@ def __init__(
scale_radius=scale_radius,
)

@aa.grid_dec.to_vector_yx
@aa.grid_dec.transform
def deflections_2d_via_integral_from(
self, grid: aa.type.Grid2DLike, xp=np, **kwargs
):
"""
Calculate the deflection angles at a given set of arc-second gridded coordinates.

Parameters
----------
grid
The grid of (y,x) arc-second coordinates the deflection angles are computed on.
"""

eta = np.multiply(
1.0 / self.scale_radius, self.radial_grid_from(grid, **kwargs).array
)

deflection_grid = np.zeros(grid.shape[0])

for i in range(grid.shape[0]):
deflection_grid[i] = np.multiply(
4.0 * self.kappa_s * self.scale_radius, self.deflection_func_sph(eta[i])
)

return self._cartesian_grid_via_radial_from(
grid=grid, radius=deflection_grid, xp=xp
)

@staticmethod
def deflection_integrand(y, eta, inner_slope):
return (y + eta) ** (inner_slope - 3) * ((1 - np.sqrt(1 - y**2)) / y)

def deflection_func_sph(self, eta):

from scipy import special
from scipy.integrate import quad

integral_y_2 = quad(
self.deflection_integrand,
a=0.0,
b=1.0,
args=(eta, self.inner_slope),
epsrel=1.49e-6,
)[0]
return eta ** (2 - self.inner_slope) * (
(1.0 / (3 - self.inner_slope))
* special.hyp2f1(
3 - self.inner_slope, 3 - self.inner_slope, 4 - self.inner_slope, -eta
)
+ integral_y_2
)
pass
Loading
Loading