From 6789cfe8078a752e8e4cfeb96a21a085a679ec5d Mon Sep 17 00:00:00 2001
From: Adrian Price-Whelan <583379+adrn@users.noreply.github.com>
Date: Tue, 21 Oct 2025 20:29:47 -0400
Subject: [PATCH 1/6] implement new tests to validate density-potential
 consistency

---
 .../potential/potential/potential_helpers.py  | 83 ++++++++++++++++++-
 tests/potential/potential/test_all_builtin.py | 46 +++++-----
 2 files changed, 102 insertions(+), 27 deletions(-)

diff --git a/tests/potential/potential/potential_helpers.py b/tests/potential/potential/potential_helpers.py
index 2c5539d16..b58651deb 100644
--- a/tests/potential/potential/potential_helpers.py
+++ b/tests/potential/potential/potential_helpers.py
@@ -6,11 +6,12 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pytest
+from astropy.constants import G
 
 from gala._compat_utils import SCIPY_LT_1_15
 from gala._optional_deps import HAS_SYMPY
 from gala.dynamics import PhaseSpacePosition
-from gala.potential import Hamiltonian, StaticFrame
+from gala.potential import Hamiltonian, SphericalSymmetry, StaticFrame
 from gala.potential.potential.io import load
 from gala.units import DimensionlessUnitSystem, UnitSystem
 
@@ -65,6 +66,9 @@ class PotentialTestBase:
     check_zero_at_infinity = True
     rotation = False
 
+    skip_hessian = False
+    skip_hessian_density = False
+
     def setup_method(self):
         # set up hamiltonian
         if self.frame is None:
@@ -159,6 +163,9 @@ def test_gradient(self):
             )
 
     def test_hessian(self):
+        if self.skip_hessian:
+            pytest.skip("Hessian not implemented for this potential")
+
         for arr, shp in zip(self.w0s, self._hess_return_shapes):
             g = self.potential.hessian(arr[: self.ndim])
             assert g.shape == shp
@@ -322,6 +329,80 @@ def compute_energy(xyz):
         ).T
         np.testing.assert_allclose(grad, num_grad, rtol=self.tol, atol=1e-8)
 
+    def test_numerical_density_vs_density(self):
+        """
+        Compare a numerically estimated Laplacian (trace of Hessian) times 4*pi*G
+        to the implemented density function via Poisson's equation:
+            Laplacian(Phi) = 4*pi*G * rho
+        """
+        if not isinstance(self.potential._symmetry, SphericalSymmetry):
+            pytest.skip("only for spherical potentials")
+
+        dx = 1e-3 * np.sqrt(np.sum(self.w0[: self.w0.size // 2] ** 2))
+        max_x = np.sqrt(np.sum([x**2 for x in self.w0[: self.w0.size // 2]]))
+
+        # use a slightly smaller grid to limit cost of nested derivatives
+        grid = np.linspace(-max_x, max_x, 6)
+        grid = grid[grid != 0.0]
+        grids = [grid for i in range(self.w0.size // 2)]
+        xyz = np.ascontiguousarray(
+            np.vstack(list(map(np.ravel, np.meshgrid(*grids)))).T
+        )
+
+        def compute_energy(xyz):
+            xyz = np.ascontiguousarray(np.atleast_2d(xyz))
+            return np.squeeze(self.potential._energy(xyz, t=np.array([0.0])))
+
+        def numeric_laplacian_at(pt):
+            # sum of second partial derivatives via nested numerical differentiation
+            s = 0.0
+            ndim = self.w0.size // 2
+            for dim in range(ndim):
+
+                def first_deriv(x):
+                    return partial_derivative(
+                        compute_energy,
+                        x,
+                        dim_ix=dim,  # noqa: B023
+                        dx=dx,
+                    )
+
+                s += partial_derivative(first_deriv, pt, dim_ix=dim, dx=dx)
+            return s
+
+        num_lap = np.zeros(xyz.shape[0])
+        for i in range(xyz.shape[0]):
+            num_lap[i] = numeric_laplacian_at(xyz[i])
+
+        num_rho = num_lap / (4.0 * np.pi * self.potential.G)
+
+        try:
+            dens = np.squeeze(
+                self.potential._density(np.ascontiguousarray(xyz), t=np.array([0.0]))
+            )
+        except Exception:
+            pytest.skip("density not implemented for this potential")
+
+        np.testing.assert_allclose(dens, num_rho, rtol=self.tol, atol=1e-6)
+
+    def test_hessian_density_consistency(self):
+        """
+        Check that the trace of the Hessian matches the density via Poisson's equation
+        """
+        if self.skip_hessian or self.skip_hessian_density:
+            pytest.skip("Hessian not implemented for this potential")
+
+        _G = G if not isinstance(self.potential.units, DimensionlessUnitSystem) else 1.0
+
+        for arr in self.w0s:
+            hess = self.potential.hessian(arr[: self.ndim])
+            lap = np.sum(np.diagonal(hess, axis1=0, axis2=1), axis=-1)
+
+            dens = self.potential.density(arr[: self.ndim])
+            rho_from_hess = lap / (4.0 * np.pi * _G)
+            print("YO", dens, rho_from_hess, np.diagonal(hess, axis1=0, axis2=1))
+            assert u.allclose(dens, rho_from_hess, rtol=self.tol)
+
     def test_orbit_integration(self, t1=0.0, t2=1000.0, nsteps=10000):
         """
         Make sure we can integrate an orbit in this potential
diff --git a/tests/potential/potential/test_all_builtin.py b/tests/potential/potential/test_all_builtin.py
index 3491955c1..9cc7177cf 100644
--- a/tests/potential/potential/test_all_builtin.py
+++ b/tests/potential/potential/test_all_builtin.py
@@ -34,6 +34,7 @@ class TestHarmonicOscillator1D(PotentialTestBase):
     sympy_density = False
     check_finite_at_origin = False
     check_zero_at_infinity = False
+    skip_hessian_density = True
 
     def test_plot(self):
         # Skip for now because contour plotting assumes 3D
@@ -46,6 +47,7 @@ class TestHarmonicOscillator2D(PotentialTestBase):
     sympy_density = False
     check_finite_at_origin = False
     check_zero_at_infinity = False
+    skip_hessian_density = True
 
     def test_plot(self):
         # Skip for now because contour plotting assumes 3D
@@ -64,6 +66,7 @@ def test_against_sympy(self):
 class TestNull(PotentialTestBase):
     potential = p.NullPotential()
     w0 = [1.0, 0.0, 0.0, 0.0, 2 * np.pi, 0.0]
+    skip_hessian_density = True
 
     def test_mass_enclosed(self):
         for arr, shp in zip(self.w0s, self._valu_return_shapes):
@@ -110,6 +113,7 @@ class TestHenonHeiles(PotentialTestBase):
     sympy_density = False
     check_finite_at_origin = False
     check_zero_at_infinity = False
+    skip_hessian_density = True
 
     @pytest.mark.skip(reason="Not relevant")
     def test_plot(self):
@@ -119,14 +123,15 @@ def test_plot(self):
 class TestKepler(PotentialTestBase):
     potential = p.KeplerPotential(units=solarsystem, m=1.0)
     w0 = [1.0, 0.0, 0.0, 0.0, 2 * np.pi, 0.0]
-    # show_plots = True
     check_finite_at_origin = False
+    skip_hessian_density = True
 
 
 class TestKeplerUnitInput(PotentialTestBase):
     potential = p.KeplerPotential(units=solarsystem, m=(1 * u.Msun).to(u.Mjup))
     w0 = [1.0, 0.0, 0.0, 0.0, 2 * np.pi, 0.0]
     check_finite_at_origin = False
+    skip_hessian_density = True
 
 
 class TestIsochrone(PotentialTestBase):
@@ -228,6 +233,7 @@ class TestKuzmin(PotentialTestBase):
     sympy_hessian = False
     sympy_density = False
     rotation = True
+    skip_hessian = True  # TODO: implement
 
 
 class TestStone(PotentialTestBase):
@@ -244,7 +250,7 @@ class TestPowerLawCutoff(PotentialTestBase):
 
     def setup_method(self):
         self.potential = p.PowerLawCutoffPotential(
-            units=galactic, m=1e10, r_c=1.0, alpha=1.8
+            units=galactic, m=1e10, r_c=10.0, alpha=1.8
         )
         super().setup_method()
 
@@ -259,6 +265,7 @@ class TestFlattenedNFW(PotentialTestBase):
     w0 = [19.0, 2.7, -6.9, 0.0352238, -0.03579493, 0.075]
     sympy_density = False  # not defined
     rotation = True
+    skip_hessian_density = True  # no density defined
 
     def test_against_spherical(self):
         """
@@ -276,6 +283,7 @@ class TestTriaxialNFW(PotentialTestBase):
     w0 = [19.0, 2.7, -6.9, 0.0352238, -0.03579493, 0.075]
     sympy_density = False  # not defined
     rotation = True
+    skip_hessian_density = True  # no density defined
 
 
 class TestSphericalNFWFromCircVel(PotentialTestBase):
@@ -342,6 +350,7 @@ class TestNFW(PotentialTestBase):
     )
     w0 = [19.0, 2.7, -0.9, 0.00352238, -0.15134, 0.0075]
     sympy_density = False  # like triaxial case
+    skip_hessian_density = True  # no density defined
 
     def test_compare(self):
         sph = p.NFWPotential(m=6e11 * u.Msun, r_s=20 * u.kpc, units=galactic)
@@ -416,6 +425,7 @@ class TestLeeSutoTriaxialNFW(PotentialTestBase):
     )
     w0 = [19.0, 2.7, -6.9, 0.0352238, -0.03579493, 0.075]
     rotation = True
+    skip_hessian = True  # TODO: implement
 
     @pytest.mark.skip(reason="to_sympy() not implemented yet")
     def test_against_sympy(self):
@@ -457,10 +467,7 @@ class TestLongMuraliBarRotate(PotentialTestBase):
     vc = potential.circular_velocity([19.0, 0, 0] * u.kpc).decompose(galactic).value[0]
     w0 = [19.0, 0.2, -0.9, 0.0, vc, 0.0]
 
-    def test_hessian(self):
-        # TODO: when hessian for rotated potentials implemented, remove this
-        with pytest.raises(NotImplementedError):
-            self.potential.hessian([1.0, 2.0, 3.0])
+    skip_hessian = True  # TODO: implement
 
     @pytest.mark.skip(reason="Not implemented for rotated potentials")
     def test_against_sympy(self):
@@ -479,10 +486,7 @@ class TestLongMuraliBarRotationScipy(PotentialTestBase):
     vc = potential.circular_velocity([19.0, 0, 0] * u.kpc).decompose(galactic).value[0]
     w0 = [19.0, 0.2, -0.9, 0.0, vc, 0.0]
 
-    def test_hessian(self):
-        # TODO: when hessian for rotated potentials implemented, remove this
-        with pytest.raises(NotImplementedError):
-            self.potential.hessian([1.0, 2.0, 3.0])
+    skip_hessian = True  # TODO: implement
 
     @pytest.mark.skip(reason="Not implemented for rotated potentials")
     def test_against_sympy(self):
@@ -531,6 +535,8 @@ class TestKepler3Body(CompositePotentialTestBase):
     frame = ConstantRotatingFrame(Omega=Omega)
     w0 = [0.5, 0, 0, 0.0, 1.05800316, 0.0]
 
+    skip_hessian_density = True
+
 
 @pytest.mark.skipif(not GSL_ENABLED, reason="requires GSL to run this test")
 class TestMultipoleInner(CompositePotentialTestBase):
@@ -541,10 +547,7 @@ class TestMultipoleInner(CompositePotentialTestBase):
     vc = potential.circular_velocity([19.0, 0, 0] * u.kpc).decompose(galactic).value[0]
     w0 = [19.0, 0.2, -0.9, 0.0, vc, 0.0]
     check_zero_at_infinity = False
-
-    @pytest.mark.skip(reason="Not implemented for multipole potentials")
-    def test_hessian(self):
-        pass
+    skip_hessian = True  # TODO: implement
 
     @pytest.mark.skip(reason="Not implemented for multipole potentials")
     def test_against_sympy(self):
@@ -560,10 +563,7 @@ class TestMultipoleOuter(CompositePotentialTestBase):
     vc = potential.circular_velocity([19.0, 0, 0] * u.kpc).decompose(galactic).value[0]
     w0 = [19.0, 0.2, -0.9, 0.0, vc, 0.0]
     check_finite_at_origin = False
-
-    @pytest.mark.skip(reason="Not implemented for multipole potentials")
-    def test_hessian(self):
-        pass
+    skip_hessian = True  # TODO: implement
 
     @pytest.mark.skip(reason="Not implemented for multipole potentials")
     def test_against_sympy(self):
@@ -573,6 +573,7 @@ def test_against_sympy(self):
 @pytest.mark.skipif(not GSL_ENABLED, reason="requires GSL to run this test")
 class TestCylspline(PotentialTestBase):
     check_finite_at_origin = True
+    skip_hessian = True  # TODO: implement
 
     def setup_method(self):
         self.potential = p.CylSplinePotential.from_file(
@@ -586,10 +587,6 @@ def setup_method(self):
     def test_density(self):
         pass
 
-    @pytest.mark.skip(reason="Not implemented for cylspline potentials")
-    def test_hessian(self):
-        pass
-
     @pytest.mark.skip(reason="Not implemented for cylspline potentials")
     def test_against_sympy(self):
         pass
@@ -612,15 +609,12 @@ class TestBurkert(PotentialTestBase):
     w0 = [1.0, 0.0, 0.0, 0.0, 0.1, 0.1]
 
     check_finite_at_origin = False
+    skip_hessian = True  # TODO: implement
 
     @pytest.mark.skip(reason="Not implemented for Burkert potentials")
     def test_against_sympy(self):
         pass
 
-    @pytest.mark.skip(reason="Hessian not implemented for Burkert potential")
-    def test_hessian(self):
-        pass
-
     def test_from_r0(self):
         # Test against values from Zhu+2023
         pot = p.BurkertPotential.from_r0(r0=11.87 * u.kpc, units=galactic)

From 053881bd9653e33558f6acbd3bbe7494e2669f81 Mon Sep 17 00:00:00 2001
From: Adrian Price-Whelan <583379+adrn@users.noreply.github.com>
Date: Wed, 22 Oct 2025 19:03:38 -0400
Subject: [PATCH 2/6] add new potential_helpers header with common functions

---
 src/gala/potential/potential/builtin/potential_helpers.h | 9 +++++++++
 1 file changed, 9 insertions(+)
 create mode 100644 src/gala/potential/potential/builtin/potential_helpers.h

diff --git a/src/gala/potential/potential/builtin/potential_helpers.h b/src/gala/potential/potential/builtin/potential_helpers.h
new file mode 100644
index 000000000..1e35cb67a
--- /dev/null
+++ b/src/gala/potential/potential/builtin/potential_helpers.h
@@ -0,0 +1,9 @@
+#include <math.h>
+
+static inline double norm2(const double *q) {
+    return sqrt(q[0]*q[0] + q[1]*q[1]);
+}
+
+static inline double norm3(const double *q) {
+    return sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+}

From 3a78c83cc2276bc44cd33c34f60030769016c457 Mon Sep 17 00:00:00 2001
From: Adrian Price-Whelan <583379+adrn@users.noreply.github.com>
Date: Wed, 22 Oct 2025 20:04:45 -0400
Subject: [PATCH 3/6] add helper functions and use constant double

---
 .../potential/builtin/builtin_potentials.cpp  | 393 ++++++++----------
 .../potential/builtin/potential_helpers.h     |  30 +-
 2 files changed, 197 insertions(+), 226 deletions(-)

diff --git a/src/gala/potential/potential/builtin/builtin_potentials.cpp b/src/gala/potential/potential/builtin/builtin_potentials.cpp
index 3204e4260..6c3d71abe 100644
--- a/src/gala/potential/potential/builtin/builtin_potentials.cpp
+++ b/src/gala/potential/potential/builtin/builtin_potentials.cpp
@@ -3,6 +3,7 @@
 #include <stdio.h>
 #include "extra_compile_macros.h"
 #include "src/vectorization.h"
+#include "potential_helpers.h"
 
 #if USE_GSL == 1
 #include <gsl/gsl_sf_gamma.h>
@@ -59,9 +60,7 @@ double kepler_value(double t, double *pars, double *q, int n_dim, void *state) {
             - G (Gravitational constant)
             - m (mass scale)
     */
-    double R;
-    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    return -pars[0] * pars[1] / R;
+    return -pars[0] * pars[1] / norm3(q);
 }
 
 void kepler_gradient_single(double t, double *__restrict__ pars, double6ptr q, int n_dim, double6ptr grad, void *__restrict__ state) {
@@ -69,9 +68,7 @@ void kepler_gradient_single(double t, double *__restrict__ pars, double6ptr q, i
             - G (Gravitational constant)
             - m (mass scale)
     */
-    double R, fac;
-    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    fac = pars[0] * pars[1] / (R*R*R);
+    const double fac = pars[0] * pars[1] / pow(norm3(q), 3);
 
     grad[0] = grad[0] + fac*q[0];
     grad[1] = grad[1] + fac*q[1];
@@ -83,10 +80,7 @@ double kepler_density(double t, double *pars, double *q, int n_dim, void *state)
             - G (Gravitational constant)
             - m (mass scale)
     */
-    double r2;
-    r2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
-
-    if (r2 == 0.) {
+    if (norm3_sq(q) == 0.) {
         return INFINITY;
     } else {
         return 0.;
@@ -136,9 +130,8 @@ double isochrone_value(double t, double *pars, double *q, int n_dim, void *state
             - m (mass scale)
             - b (core scale)
     */
-    double R2;
-    R2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
-    return -pars[0] * pars[1] / (sqrt(R2 + pars[2]*pars[2]) + pars[2]);
+    const double r2 = norm3_sq(q);
+    return -pars[0] * pars[1] / (sqrt(r2 + pars[2]*pars[2]) + pars[2]);
 }
 
 void isochrone_gradient_single(double t, double *__restrict__ pars, double6ptr q, int n_dim, double6ptr grad, void *__restrict__ state) {
@@ -147,10 +140,9 @@ void isochrone_gradient_single(double t, double *__restrict__ pars, double6ptr q
             - m (mass scale)
             - b (core scale)
     */
-    double sqrt_r2_b2, fac, denom;
-    sqrt_r2_b2 = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + pars[2]*pars[2]);
-    denom = sqrt_r2_b2 * (sqrt_r2_b2 + pars[2])*(sqrt_r2_b2 + pars[2]);
-    fac = pars[0] * pars[1] / denom;
+    const double sqrt_r2_b2 = sqrt(norm3_sq(q) + pars[2]*pars[2]);
+    const double denom = sqrt_r2_b2 * (sqrt_r2_b2 + pars[2])*(sqrt_r2_b2 + pars[2]);
+    const double fac = pars[0] * pars[1] / denom;
 
     grad[0] = grad[0] + fac*q[0];
     grad[1] = grad[1] + fac*q[1];
@@ -163,10 +155,9 @@ double isochrone_density(double t, double *pars, double *q, int n_dim, void *sta
             - m (mass scale)
             - b (core scale)
     */
-    double r2, a, b;
-    b = pars[2];
-    r2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
-    a = sqrt(b*b + r2);
+    const double b = pars[2];
+    const double r2 = norm3_sq(q);
+    const double a = sqrt(b*b + r2);
 
     return pars[1] * (3*(b+a)*a*a - r2*(b+3*a)) / (4*M_PI*pow(b+a,3)*a*a*a);
 }
@@ -223,9 +214,8 @@ double hernquist_value(double t, double *pars, double *q, int n_dim, void *state
             - m (mass scale)
             - c (length scale)
     */
-    double R;
-    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    return -pars[0] * pars[1] / (R + pars[2]);
+    const double r = norm3(q);
+    return -pars[0] * pars[1] / (r + pars[2]);
 }
 
 void hernquist_gradient_single(double t, double *__restrict__ pars, double6ptr q, int n_dim, double6ptr grad, void *__restrict__ state) {
@@ -234,9 +224,8 @@ void hernquist_gradient_single(double t, double *__restrict__ pars, double6ptr q
             - m (mass scale)
             - c (length scale)
     */
-    double R, fac;
-    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    fac = pars[0] * pars[1] / ((R + pars[2]) * (R + pars[2]) * R);
+    const double r = norm3(q);
+    const double fac = pars[0] * pars[1] / ((r + pars[2]) * (r + pars[2]) * r);
 
     grad[0] = grad[0] + fac*q[0];
     grad[1] = grad[1] + fac*q[1];
@@ -249,9 +238,8 @@ double hernquist_density(double t, double *pars, double *q, int n_dim, void *sta
             - m (mass scale)
             - c (length scale)
     */
-    double r, rho0;
-    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    rho0 = pars[1]/(2*M_PI*pars[2]*pars[2]*pars[2]);
+    const double r = norm3(q);
+    const double rho0 = pars[1]/(2*M_PI*pars[2]*pars[2]*pars[2]);
     return rho0 / ((r/pars[2]) * pow(1+r/pars[2],3));
 }
 
@@ -307,8 +295,8 @@ double plummer_value(double t, double *pars, double *q, int n_dim, void *state)
             - m (mass scale)
             - b (length scale)
     */
-    double R2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
-    return -pars[0]*pars[1] / sqrt(R2 + pars[2]*pars[2]);
+    const double r2 = norm3_sq(q);
+    return -pars[0]*pars[1] / sqrt(r2 + pars[2]*pars[2]);
 }
 
 void plummer_gradient_single(double t, double *__restrict__ pars, double6ptr q, int n_dim, double6ptr grad, void *__restrict__ state) {
@@ -317,9 +305,8 @@ void plummer_gradient_single(double t, double *__restrict__ pars, double6ptr q,
             - m (mass scale)
             - b (length scale)
     */
-    double R2b, fac;
-    R2b = q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + pars[2]*pars[2];
-    fac = pars[0] * pars[1] / sqrt(R2b) / R2b;
+    const double R2b = norm3_sq(q) + pars[2]*pars[2];
+    const double fac = pars[0] * pars[1] / sqrt(R2b) / R2b;
 
     grad[0] = grad[0] + fac*q[0];
     grad[1] = grad[1] + fac*q[1];
@@ -332,7 +319,7 @@ double plummer_density(double t, double *pars, double *q, int n_dim, void *state
             - m (mass scale)
             - b (length scale)
     */
-    double r2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
+    const double r2 = norm3_sq(q);
     return 3*pars[1] / (4*M_PI*pars[2]*pars[2]*pars[2]) * pow(1 + r2/(pars[2]*pars[2]), -2.5);
 }
 
@@ -381,9 +368,8 @@ double jaffe_value(double t, double *pars, double *q, int n_dim, void *state) {
             - m (mass scale)
             - c (length scale)
     */
-    double R;
-    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    return -pars[0] * pars[1] / pars[2] * log(1 + pars[2]/R);
+    const double r = norm3(q);
+    return -pars[0] * pars[1] / pars[2] * log(1 + pars[2] / r);
 }
 
 void jaffe_gradient_single(double t, double *__restrict__ pars, double6ptr q, int n_dim, double6ptr grad, void *__restrict__ state){
@@ -392,13 +378,12 @@ void jaffe_gradient_single(double t, double *__restrict__ pars, double6ptr q, in
             - m (mass scale)
             - c (length scale)
     */
-    double R, fac;
-    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    fac = pars[0] * pars[1] / pars[2] * (pars[2] / (R * (pars[2] + R)));
+    const double r = norm3(q);
+    const double fac = pars[0] * pars[1] / pars[2] * (pars[2] / (r * (pars[2] + r))) / r;
 
-    grad[0] = grad[0] + fac*q[0]/R;
-    grad[1] = grad[1] + fac*q[1]/R;
-    grad[2] = grad[2] + fac*q[2]/R;
+    grad[0] = grad[0] + fac * q[0];
+    grad[1] = grad[1] + fac * q[1];
+    grad[2] = grad[2] + fac * q[2];
 }
 
 double jaffe_density(double t, double *pars, double *q, int n_dim, void *state) {
@@ -407,9 +392,8 @@ double jaffe_density(double t, double *pars, double *q, int n_dim, void *state)
             - m (mass scale)
             - c (length scale)
     */
-    double r, rho0;
-    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    rho0 = pars[1] / (4*M_PI*pars[2]*pars[2]*pars[2]);
+    const double r = norm3(q);
+    const double rho0 = pars[1] / (4*M_PI*pars[2]*pars[2]*pars[2]);
     return rho0 / (pow(r/pars[2],2) * pow(1+r/pars[2],2));
 }
 
@@ -511,25 +495,22 @@ double powerlawcutoff_value(double t, double *pars, double *q, int n_dim, void *
             2 - a (power-law index)
             3 - c (cutoff radius)
     */
-    double G = pars[0];
-    double m = pars[1];
-    double alpha = pars[2];
-    double r_c = pars[3];
-    double x = q[0];
-    double y = q[1];
-    double z = q[2];
-    double r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    const double G = pars[0];
+    const double m = pars[1];
+    const double alpha = pars[2];
+    const double r_c = pars[3];
+    const double r = norm3(q);
 
     if (r == 0.) {
         return -INFINITY;
     } else {
-        double tmp_0 = alpha / 2.0;
-        double tmp_1 = -tmp_0;
-        double tmp_2 = tmp_1 + 1.5;
-        double tmp_3 = pow(x, 2) + pow(y, 2) + pow(z, 2);
-        double tmp_4 = tmp_3/pow(r_c, 2);
-        double tmp_5 = G*m;
-        double tmp_6 = tmp_5*safe_gamma_inc(tmp_2, tmp_4)/(sqrt(tmp_3)*tgamma(tmp_1 + 2.5));
+        const double tmp_0 = alpha / 2.0;
+        const double tmp_1 = -tmp_0;
+        const double tmp_2 = tmp_1 + 1.5;
+        const double tmp_3 = r * r;
+        const double tmp_4 = tmp_3 / pow(r_c, 2);
+        const double tmp_5 = G*m;
+        const double tmp_6 = tmp_5*safe_gamma_inc(tmp_2, tmp_4)/(sqrt(tmp_3)*tgamma(tmp_1 + 2.5));
 
         // Original potential
         double phi_r = tmp_0*tmp_6 - 3.0/2.0*tmp_6 + tmp_5*safe_gamma_inc(tmp_1 + 1, tmp_4)/(r_c*tgamma(tmp_2));
@@ -551,9 +532,8 @@ double powerlawcutoff_density(double t, double *pars, double *q, int n_dim, void
             2 - a (power-law index)
             3 - c (cutoff radius)
     */
-    double r, A;
-    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    A = pars[1] / (2*M_PI) * pow(pars[3], pars[2] - 3) / gsl_sf_gamma(0.5 * (3 - pars[2]));
+    const double r = norm3(q);
+    const double A = pars[1] / (2*M_PI) * pow(pars[3], pars[2] - 3) / gsl_sf_gamma(0.5 * (3 - pars[2]));
     return A * pow(r, -pars[2]) * exp(-r*r / (pars[3]*pars[3]));
 }
 
@@ -564,14 +544,13 @@ void powerlawcutoff_gradient_single(double t, double *__restrict__ pars, double6
             2 - a (power-law index)
             3 - c (cutoff radius)
     */
-    double r, dPhi_dr;
-    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    dPhi_dr = (pars[0] * pars[1] / (r*r) *
+    const double r = norm3(q);
+    const double dPhi_dr = (pars[0] * pars[1] / (r*r * r) *
         gsl_sf_gamma_inc_P(0.5 * (3-pars[2]), r*r/(pars[3]*pars[3]))); // / gsl_sf_gamma(0.5 * (3-pars[2])));
 
-    grad[0] = grad[0] + dPhi_dr * q[0]/r;
-    grad[1] = grad[1] + dPhi_dr * q[1]/r;
-    grad[2] = grad[2] + dPhi_dr * q[2]/r;
+    grad[0] = grad[0] + dPhi_dr * q[0];
+    grad[1] = grad[1] + dPhi_dr * q[1];
+    grad[2] = grad[2] + dPhi_dr * q[2];
 }
 
 void powerlawcutoff_hessian(double t, double *pars, double *q, int n_dim, double *hess, void *state) {
@@ -687,13 +666,11 @@ double stone_value(double t, double *pars, double *q, int n_dim, void *state) {
             - r_c (core radius)
             - r_h (halo radius)
     */
-    double r, u_c, u_h, fac;
-
-    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    u_c = r / pars[2];
-    u_h = r / pars[3];
+    const double r = norm3(q);
+    const double u_c = r / pars[2];
+    const double u_h = r / pars[3];
 
-    fac = 2*pars[0]*pars[1] / M_PI / (pars[3] - pars[2]);
+    const double fac = 2*pars[0]*pars[1] / M_PI / (pars[3] - pars[2]);
 
     if (r == 0) {
         return -fac * 0.5 * log(pars[3]*pars[3] / (pars[2] * pars[2]));
@@ -713,18 +690,17 @@ void stone_gradient_single(double t, double *__restrict__ pars, double6ptr q, in
             - r_c (core radius)
             - r_h (halo radius)
     */
-    double r, u_c, u_h, fac, dphi_dr;
 
-    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    u_c = r / pars[2];
-    u_h = r / pars[3];
+    const double r = norm3(q);
+    const double u_c = r / pars[2];
+    const double u_h = r / pars[3];
 
-    fac = 2*pars[0]*pars[1] / (M_PI*r*r) / (pars[2] - pars[3]);  // order flipped from value
-    dphi_dr = fac * (pars[2]*atan(u_c) - pars[3]*atan(u_h));
+    const double fac = 2*pars[0]*pars[1] / (M_PI*r*r * r) / (pars[2] - pars[3]);  // order flipped from value
+    const double dphi_dr = fac * (pars[2]*atan(u_c) - pars[3]*atan(u_h));
 
-    grad[0] = grad[0] + dphi_dr*q[0]/r;
-    grad[1] = grad[1] + dphi_dr*q[1]/r;
-    grad[2] = grad[2] + dphi_dr*q[2]/r;
+    grad[0] = grad[0] + dphi_dr*q[0];
+    grad[1] = grad[1] + dphi_dr*q[1];
+    grad[2] = grad[2] + dphi_dr*q[2];
 }
 
 double stone_density(double t, double *pars, double *q, int n_dim, void *state) {
@@ -734,12 +710,10 @@ double stone_density(double t, double *pars, double *q, int n_dim, void *state)
             - r_c (core radius)
             - r_h (halo radius)
     */
-    double r, u_c, u_t, rho;
-    rho = pars[1] * (pars[2] + pars[3]) / (2*M_PI*M_PI*pars[2]*pars[2]*pars[3]*pars[3]);
-
-    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    u_c = r / pars[2];
-    u_t = r / pars[3];
+    const double r = norm3(q);
+    const double rho = pars[1] * (pars[2] + pars[3]) / (2*M_PI*M_PI*pars[2]*pars[2]*pars[3]*pars[3]);
+    const double u_c = r / pars[2];
+    const double u_t = r / pars[3];
 
     return rho / ((1 + u_c*u_c)*(1 + u_t*u_t));
 }
@@ -848,10 +822,9 @@ double sphericalnfw_value(double t, double *pars, double *q, int n_dim, void *st
             - m (mass scale)
             - r_s (scale radius)
     */
-    double u, v_h2;
     // v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
-    v_h2 = -pars[0] * pars[1] / pars[2];
-    u = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]) / pars[2];
+    const double v_h2 = -pars[0] * pars[1] / pars[2];
+    const double u = norm3(q) / pars[2];
     if (u == 0) {
         return v_h2;
     } else {
@@ -865,12 +838,11 @@ void sphericalnfw_gradient_single(double t, double *__restrict__ pars, double6pt
             - m (mass scale)
             - r_s (scale radius)
     */
-    double fac, u, v_h2;
     // v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
-    v_h2 = pars[0] * pars[1] / pars[2];
+    const double v_h2 = pars[0] * pars[1] / pars[2];
 
-    u = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]) / pars[2];
-    fac = v_h2 / (u*u*u) / (pars[2]*pars[2]) * (log(1+u) - u/(1+u));
+    const double u = norm3(q) / pars[2];
+    const double fac = v_h2 / (u*u*u) / (pars[2]*pars[2]) * (log(1+u) - u/(1+u));
 
     grad[0] = grad[0] + fac*q[0];
     grad[1] = grad[1] + fac*q[1];
@@ -884,11 +856,10 @@ double sphericalnfw_density(double t, double *pars, double *q, int n_dim, void *
             - r_s (scale radius)
     */
     // double v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
-    double v_h2 = pars[0] * pars[1] / pars[2];
-    double r, rho0;
-    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    const double v_h2 = pars[0] * pars[1] / pars[2];
+    const double r = norm3(q);
 
-    rho0 = v_h2 / (4*M_PI*pars[0]*pars[2]*pars[2]);
+    const double rho0 = v_h2 / (4*M_PI*pars[0]*pars[2]*pars[2]);
     return rho0 / ((r/pars[2]) * pow(1+r/pars[2],2));
 }
 
@@ -960,10 +931,9 @@ double flattenednfw_value(double t, double *pars, double *q, int n_dim, void *st
             - b (ignore)
             - c (z flattening)
     */
-    double u, v_h2;
     // v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
-    v_h2 = -pars[0] * pars[1] / pars[2];
-    u = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]/(pars[5]*pars[5])) / pars[2];
+    const double v_h2 = -pars[0] * pars[1] / pars[2];
+    const double u = norm3_flat_z(q, pars[5]) / pars[2];
     if (u == 0) {
         return v_h2;
     } else {
@@ -980,12 +950,11 @@ void flattenednfw_gradient_single(double t, double *__restrict__ pars, double6pt
             - b (ignore)
             - c (z flattening)
     */
-    double fac, u, v_h2;
     // v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
-    v_h2 = pars[0] * pars[1] / pars[2];
-    u = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]/(pars[5]*pars[5])) / pars[2];
+    const double v_h2 = pars[0] * pars[1] / pars[2];
+    const double u = norm3_flat_z(q, pars[5]) / pars[2];
 
-    fac = v_h2 / (u*u*u) / (pars[2]*pars[2]) * (log(1+u) - u/(1+u));
+    const double fac = v_h2 / (u*u*u) / (pars[2]*pars[2]) * (log(1+u) - u/(1+u));
 
     grad[0] = grad[0] + fac*q[0];
     grad[1] = grad[1] + fac*q[1];
@@ -1182,8 +1151,7 @@ double satoh_value(double t, double *pars, double *q, int n_dim, void *state) {
             - a (length scale 1) TODO
             - b (length scale 2) TODO
     */
-    double S2;
-    S2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + pars[2]*(pars[2] + 2*sqrt(q[2]*q[2] + pars[3]*pars[3]));
+    const double S2 = norm3_sq(q) + pars[2]*(pars[2] + 2*sqrt(q[2]*q[2] + pars[3]*pars[3]));
     return -pars[0] * pars[1] / sqrt(S2);
 }
 
@@ -1195,8 +1163,8 @@ void satoh_gradient_single(double t, double *__restrict__ pars, double6ptr q, in
             - b (length scale 2) TODO
     */
 
-    double S2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + pars[2]*(pars[2] + 2*sqrt(q[2]*q[2] + pars[3]*pars[3]));
-    double dPhi_dS = pars[0] * pars[1] / S2;
+    const double S2 = norm3_sq(q) + pars[2]*(pars[2] + 2*sqrt(q[2]*q[2] + pars[3]*pars[3]));
+    const double dPhi_dS = pars[0] * pars[1] / S2;
 
     grad[0] = grad[0] + dPhi_dS*q[0]/sqrt(S2);
     grad[1] = grad[1] + dPhi_dS*q[1]/sqrt(S2);
@@ -1323,8 +1291,7 @@ double miyamotonagai_value(double t, double *pars, double *q, int n_dim, void *s
             - a (length scale 1) TODO
             - b (length scale 2) TODO
     */
-    double zd;
-    zd = (pars[2] + sqrt(q[2]*q[2] + pars[3]*pars[3]));
+    const double zd = (pars[2] + sqrt(q[2]*q[2] + pars[3]*pars[3]));
     return -pars[0] * pars[1] / sqrt(q[0]*q[0] + q[1]*q[1] + zd*zd);
 }
 
@@ -1335,11 +1302,9 @@ void miyamotonagai_gradient_single(double t, double *__restrict__ pars, double6p
             - a (length scale 1) TODO
             - b (length scale 2) TODO
     */
-    double sqrtz, zd, fac;
-
-    sqrtz = sqrt(q[2]*q[2] + pars[3]*pars[3]);
-    zd = pars[2] + sqrtz;
-    fac = pars[0]*pars[1] * pow(q[0]*q[0] + q[1]*q[1] + zd*zd, -1.5);
+    const double sqrtz = sqrt(q[2]*q[2] + pars[3]*pars[3]);
+    const double zd = pars[2] + sqrtz;
+    const double fac = pars[0]*pars[1] * pow(q[0]*q[0] + q[1]*q[1] + zd*zd, -1.5);
 
     grad[0] = grad[0] + fac*q[0];
     grad[1] = grad[1] + fac*q[1];
@@ -1354,17 +1319,16 @@ double miyamotonagai_density(double t, double *pars, double *q, int n_dim, void
             - b (length scale 2) TODO
     */
 
-    double M, a, b;
-    M = pars[1];
-    a = pars[2];
-    b = pars[3];
+    const double M = pars[1];
+    const double a = pars[2];
+    const double b = pars[3];
 
-    double R2 = q[0]*q[0] + q[1]*q[1];
-    double sqrt_zb = sqrt(q[2]*q[2] + b*b);
-    double numer = (b*b*M / (4*M_PI)) * (a*R2 + (a + 3*sqrt_zb)*(a + sqrt_zb)*(a + sqrt_zb));
-    double denom = pow(R2 + (a + sqrt_zb)*(a + sqrt_zb), 2.5) * sqrt_zb*sqrt_zb*sqrt_zb;
+    const double R2 = q[0]*q[0] + q[1]*q[1];
+    const double sqrt_zb = sqrt(q[2]*q[2] + b*b);
+    const double numer = (b*b*M / (4*M_PI)) * (a*R2 + (a + 3*sqrt_zb)*(a + sqrt_zb)*(a + sqrt_zb));
+    const double denom = pow(R2 + (a + sqrt_zb)*(a + sqrt_zb), 2.5) * sqrt_zb*sqrt_zb*sqrt_zb;
 
-    return numer/denom;
+    return numer / denom;
 }
 
 void miyamotonagai_hessian(double t, double *pars, double *q, int n_dim,
@@ -1602,11 +1566,9 @@ double logarithmic_value(double t, double *pars, double *q, int n_dim, void *sta
             - q2
             - q3
     */
-    double x, y, z;
-
-    x = q[0]*cos(pars[6]) + q[1]*sin(pars[6]);
-    y = -q[0]*sin(pars[6]) + q[1]*cos(pars[6]);
-    z = q[2];
+    const double x = q[0]*cos(pars[6]) + q[1]*sin(pars[6]);
+    const double y = -q[0]*sin(pars[6]) + q[1]*cos(pars[6]);
+    const double z = q[2];
 
     return 0.5*pars[1]*pars[1] * log(pars[2]*pars[2] + // scale radius
                                      x*x/(pars[3]*pars[3]) +
@@ -1647,16 +1609,15 @@ void logarithmic_gradient_single(double t, double *__restrict__ pars, double6ptr
             - q2
             - q3
     */
-    double x, y, z, ax, ay, az, fac;
 
-    x = q[0]*cos(pars[6]) + q[1]*sin(pars[6]);
-    y = -q[0]*sin(pars[6]) + q[1]*cos(pars[6]);
-    z = q[2];
+    const double x = q[0]*cos(pars[6]) + q[1]*sin(pars[6]);
+    const double y = -q[0]*sin(pars[6]) + q[1]*cos(pars[6]);
+    const double z = q[2];
 
-    fac = pars[1]*pars[1] / (pars[2]*pars[2] + x*x/(pars[3]*pars[3]) + y*y/(pars[4]*pars[4]) + z*z/(pars[5]*pars[5]));
-    ax = fac*x/(pars[3]*pars[3]);
-    ay = fac*y/(pars[4]*pars[4]);
-    az = fac*z/(pars[5]*pars[5]);
+    const double fac = pars[1]*pars[1] / (pars[2]*pars[2] + x*x/(pars[3]*pars[3]) + y*y/(pars[4]*pars[4]) + z*z/(pars[5]*pars[5]));
+    const double ax = fac*x/(pars[3]*pars[3]);
+    const double ay = fac*y/(pars[4]*pars[4]);
+    const double az = fac*z/(pars[5]*pars[5]);
 
     grad[0] = grad[0] + (ax*cos(pars[6]) - ay*sin(pars[6]));
     grad[1] = grad[1] + (ax*sin(pars[6]) + ay*cos(pars[6]));
@@ -1728,20 +1689,16 @@ double longmuralibar_value(double t, double *pars, double *q, int n_dim, void *s
         - c
         - alpha
     */
-    double x, y, z;
-    double a, b, c;
-    double Tm, Tp;
-
-    x = q[0]*cos(pars[5]) + q[1]*sin(pars[5]);
-    y = -q[0]*sin(pars[5]) + q[1]*cos(pars[5]);
-    z = q[2];
+    const double x = q[0]*cos(pars[5]) + q[1]*sin(pars[5]);
+    const double y = -q[0]*sin(pars[5]) + q[1]*cos(pars[5]);
+    const double z = q[2];
 
-    a = pars[2];
-    b = pars[3];
-    c = pars[4];
+    const double a = pars[2];
+    const double b = pars[3];
+    const double c = pars[4];
 
-    Tm = sqrt((a-x)*(a-x) + y*y + pow(b + sqrt(c*c + z*z),2));
-    Tp = sqrt((a+x)*(a+x) + y*y + pow(b + sqrt(c*c + z*z),2));
+    const double Tm = sqrt((a-x)*(a-x) + y*y + pow(b + sqrt(c*c + z*z),2));
+    const double Tp = sqrt((a+x)*(a+x) + y*y + pow(b + sqrt(c*c + z*z),2));
 
     return pars[0]*pars[1]/(2*a) * log((x - a + Tm) / (x + a + Tp));
 }
@@ -1757,30 +1714,25 @@ void longmuralibar_gradient_single(double t, double *__restrict__ pars, double6p
         - c
         - alpha
     */
-    double x, y, z;
-    double a, b, c;
-    double Tm, Tp, fac1, fac2, fac3, bcz;
-    double gx, gy, gz;
-
-    x = q[0]*cos(pars[5]) + q[1]*sin(pars[5]);
-    y = -q[0]*sin(pars[5]) + q[1]*cos(pars[5]);
-    z = q[2];
+    const double x = q[0]*cos(pars[5]) + q[1]*sin(pars[5]);
+    const double y = -q[0]*sin(pars[5]) + q[1]*cos(pars[5]);
+    const double z = q[2];
 
-    a = pars[2];
-    b = pars[3];
-    c = pars[4];
+    const double a = pars[2];
+    const double b = pars[3];
+    const double c = pars[4];
 
-    bcz = b + sqrt(c*c + z*z);
-    Tm = sqrt((a-x)*(a-x) + y*y + bcz*bcz);
-    Tp = sqrt((a+x)*(a+x) + y*y + bcz*bcz);
+    const double bcz = b + sqrt(c*c + z*z);
+    const double Tm = sqrt((a-x)*(a-x) + y*y + bcz*bcz);
+    const double Tp = sqrt((a+x)*(a+x) + y*y + bcz*bcz);
 
-    fac1 = pars[0]*pars[1] / (2*Tm*Tp);
-    fac2 = 1 / (y*y + bcz*bcz);
-    fac3 = Tp + Tm - (4*x*x)/(Tp+Tm);
+    const double fac1 = pars[0]*pars[1] / (2*Tm*Tp);
+    const double fac2 = 1 / (y*y + bcz*bcz);
+    const double fac3 = Tp + Tm - (4*x*x)/(Tp+Tm);
 
-    gx = 4 * fac1 * x / (Tp + Tm);
-    gy = fac1 * y * fac2 * fac3;
-    gz = fac1 * z * fac2 * fac3 * bcz / sqrt(z*z + c*c);
+    const double gx = 4 * fac1 * x / (Tp + Tm);
+    const double gy = fac1 * y * fac2 * fac3;
+    const double gz = fac1 * z * fac2 * fac3 * bcz / sqrt(z*z + c*c);
 
     grad[0] = grad[0] + (gx*cos(pars[5]) - gy*sin(pars[5]));
     grad[1] = grad[1] + (gx*sin(pars[5]) + gy*cos(pars[5]));
@@ -2022,7 +1974,7 @@ double spherical_spline_density_value(double t, double *pars, double *q, int n_d
             1 to 1+n_knots-1 - r_knots (radial knot locations)
             n_knots to 2*n_knots-1 - density_values (density at each knot)
     */
-    double r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    const double r = norm3(q);
     spherical_spline_state *spl_state = (spherical_spline_state *)state;
 
     // Check bounds
@@ -2031,17 +1983,15 @@ double spherical_spline_density_value(double t, double *pars, double *q, int n_d
     }
 
     // Calculate enclosed mass M(r) = 4π ∫[0 to r] ρ(r') r'² dr'
-    double r_min = spl_state->r_knots[0];
-    double integral_mass = gsl_spline_eval_integ(spl_state->rho_r2_spline, r_min, r, spl_state->rho_r2_acc);
-    double M_r = 4.0 * M_PI * integral_mass;
+    const double r_min = spl_state->r_knots[0];
+    const double integral_mass = gsl_spline_eval_integ(spl_state->rho_r2_spline, r_min, r, spl_state->rho_r2_acc);
+    const double M_r = 4.0 * M_PI * integral_mass;
 
     // Calculate potential from density
     // For spherical symmetry: Φ(r) = -G M(r) / r - 4πG ∫[r to ∞] ρ(r') r' dr'
-    double r_max = spl_state->r_knots[spl_state->n_knots-1];
-    double integral_outer = gsl_spline_eval_integ(spl_state->rho_r_spline, r, r_max, spl_state->rho_r_acc);
-    double potential = -pars[0] * M_r / r - 4.0 * M_PI * pars[0] * integral_outer;
-
-    return potential;
+    const double r_max = spl_state->r_knots[spl_state->n_knots-1];
+    const double integral_outer = gsl_spline_eval_integ(spl_state->rho_r_spline, r, r_max, spl_state->rho_r_acc);
+    return -pars[0] * M_r / r - 4.0 * M_PI * pars[0] * integral_outer;
 }
 
 void spherical_spline_density_gradient_single(double t, double *__restrict__ pars, double6ptr q, int n_dim, double6ptr grad, void *__restrict__ state) {
@@ -2050,7 +2000,7 @@ void spherical_spline_density_gradient_single(double t, double *__restrict__ par
             1 to 1+n_knots-1 - r_knots (radial knot locations)
             n_knots to 2*n_knots-1 - density_values (density at each knot)
     */
-    double r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    const double r = norm3(q);
     if (r == 0.0) return;
 
     spherical_spline_state *spl_state = (spherical_spline_state *)state;
@@ -2062,9 +2012,9 @@ void spherical_spline_density_gradient_single(double t, double *__restrict__ par
 
     // Calculate enclosed mass M(r) = 4π ∫[0 to r] ρ(r') r'² dr'
     // Use the pre-computed ρ(r) * r² spline
-    double r_min = spl_state->r_knots[0];
-    double integral = gsl_spline_eval_integ(spl_state->rho_r2_spline, r_min, r, spl_state->rho_r2_acc);
-    double M_r = 4.0 * M_PI * integral;
+    const double r_min = spl_state->r_knots[0];
+    const double integral = gsl_spline_eval_integ(spl_state->rho_r2_spline, r_min, r, spl_state->rho_r2_acc);
+    const double M_r = 4.0 * M_PI * integral;
 
     // Gradient: dΦ/dr = GM(r)/r²
     double dPhi_dr = pars[0] * M_r / (r * r);
@@ -2081,7 +2031,7 @@ double spherical_spline_density_density(double t, double *pars, double *q, int n
             1 to 1+n_knots-1 - r_knots (radial knot locations)
             n_knots to 2*n_knots-1 - density_values (density at each knot)
     */
-    double r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    const double r = norm3(q);
     spherical_spline_state *spl_state = (spherical_spline_state *)state;
 
     // Check bounds
@@ -2102,7 +2052,7 @@ double spherical_spline_mass_value(double t, double *pars, double *q, int n_dim,
             1 to 1+n_knots-1 - r_knots (radial knot locations)
             n_knots to 2*n_knots-1 - mass_values (mass enclosed at each knot)
     */
-    double r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    const double r = norm3(q);
     spherical_spline_state *spl_state = (spherical_spline_state *)state;
 
     double M_r;
@@ -2122,12 +2072,13 @@ double spherical_spline_mass_value(double t, double *pars, double *q, int n_dim,
     // Calculate potential: Φ(r) = -G ∫[r to ∞] M(r')/r'² dr'
     // For finite extent with maximum radius r_max, we assume M(r') = M(r_max) for r' > r_max
     // So: Φ(r) = -G ∫[r to r_max] M(r')/r'² dr' - G M(r_max) / r_max
-    double r_max = spl_state->r_knots[spl_state->n_knots-1];
-    double M_max = spl_state->values[spl_state->n_knots-1];
+    const double r_max = spl_state->r_knots[spl_state->n_knots-1];
+    const double M_max = spl_state->values[spl_state->n_knots-1];
 
     // Use numerical integration from r to r_max
+    // TODO: allow number of integration points to be a parameter
     int n_integration_points = 1000;
-    double dr = (r_max - r) / n_integration_points;
+    const double dr = (r_max - r) / n_integration_points;
     double potential = 0.0;
 
     for (int i = 0; i < n_integration_points; i++) {
@@ -2148,7 +2099,7 @@ void spherical_spline_mass_gradient_single(double t, double *__restrict__ pars,
             1 to 1+n_knots-1 - r_knots (radial knot locations)
             n_knots to 2*n_knots-1 - mass_values (mass enclosed at each knot)
     */
-    double r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    const double r = norm3(q);
     if (r == 0.0) return;
 
     spherical_spline_state *spl_state = (spherical_spline_state *)state;
@@ -2164,12 +2115,12 @@ void spherical_spline_mass_gradient_single(double t, double *__restrict__ pars,
     }
 
     // Gradient: dΦ/dr = GM(r)/r²
-    double dPhi_dr = pars[0] * M_r / (r * r);
+    const double dPhi_dr = pars[0] * M_r / (r * r * r);
 
     // Convert to Cartesian gradients
-    grad[0] += dPhi_dr * q[0] / r;
-    grad[1] += dPhi_dr * q[1] / r;
-    grad[2] += dPhi_dr * q[2] / r;
+    grad[0] += dPhi_dr * q[0];
+    grad[1] += dPhi_dr * q[1];
+    grad[2] += dPhi_dr * q[2];
 }
 
 double spherical_spline_mass_density(double t, double *pars, double *q, int n_dim, void *state) {
@@ -2178,7 +2129,7 @@ double spherical_spline_mass_density(double t, double *pars, double *q, int n_di
             1 to 1+n_knots-1 - r_knots (radial knot locations)
             n_knots to 2*n_knots-1 - mass_values (mass enclosed at each knot)
     */
-    double r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    const double r = norm3(q);
     if (r == 0.0) return 0.0;
 
     spherical_spline_state *spl_state = (spherical_spline_state *)state;
@@ -2189,7 +2140,7 @@ double spherical_spline_mass_density(double t, double *pars, double *q, int n_di
     }
 
     // Calculate density using: ρ(r) = (1/4πr²) dM/dr
-    double dM_dr = gsl_spline_eval_deriv(spl_state->spline, r, spl_state->acc);
+    const double dM_dr = gsl_spline_eval_deriv(spl_state->spline, r, spl_state->acc);
     return dM_dr / (4.0 * M_PI * r * r);
 }
 
@@ -2202,13 +2153,13 @@ double spherical_spline_potential_value(double t, double *pars, double *q, int n
             1 to 1+n_knots-1 - r_knots (radial knot locations)
             n_knots to 2*n_knots-1 - potential_values (potential at each knot)
     */
-    double r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    const double r = norm3(q);
     spherical_spline_state *spl_state = (spherical_spline_state *)state;
 
     // Check bounds - extrapolate beyond grid
     if (r < spl_state->r_knots[0]) {
         // Linear extrapolation to smaller radii
-        double slope = (spl_state->values[1] - spl_state->values[0]) / (spl_state->r_knots[1] - spl_state->r_knots[0]);
+        const double slope = (spl_state->values[1] - spl_state->values[0]) / (spl_state->r_knots[1] - spl_state->r_knots[0]);
         return spl_state->values[0] + slope * (r - spl_state->r_knots[0]);
     }
     if (r > spl_state->r_knots[spl_state->n_knots-1]) {
@@ -2226,7 +2177,7 @@ void spherical_spline_potential_gradient_single(double t, double *__restrict__ p
             1 to 1+n_knots-1 - r_knots (radial knot locations)
             n_knots to 2*n_knots-1 - potential_values (potential at each knot)
     */
-    double r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    const double r = norm3(q);
     if (r == 0.0) return;
 
     spherical_spline_state *spl_state = (spherical_spline_state *)state;
@@ -2260,7 +2211,7 @@ double spherical_spline_potential_density(double t, double *pars, double *q, int
             1 to 1+n_knots-1 - r_knots (radial knot locations)
             n_knots to 2*n_knots-1 - potential_values (potential at each knot)
     */
-    double r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    const double r = norm3(q);
     if (r == 0.0) return 0.0;
 
     spherical_spline_state *spl_state = (spherical_spline_state *)state;
@@ -2272,8 +2223,8 @@ double spherical_spline_potential_density(double t, double *pars, double *q, int
 
     // Calculate density using Poisson equation: ∇²Φ = 4πGρ
     // For spherical symmetry: ρ = (1/4πG) [d²Φ/dr² + (2/r) dΦ/dr]
-    double dPhi_dr = gsl_spline_eval_deriv(spl_state->spline, r, spl_state->acc);
-    double d2Phi_dr2 = gsl_spline_eval_deriv2(spl_state->spline, r, spl_state->acc);
+    const double dPhi_dr = gsl_spline_eval_deriv(spl_state->spline, r, spl_state->acc);
+    const double d2Phi_dr2 = gsl_spline_eval_deriv2(spl_state->spline, r, spl_state->acc);
 
     return (d2Phi_dr2 + 2.0 * dPhi_dr / r) / (4.0 * M_PI * pars[0]);
 }
@@ -2290,9 +2241,8 @@ double burkert_value(double t, double *pars, double *q, int n_dim, void *state)
             - rho (mass scale)
             - r0
     */
-    double R, x;
-    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    x = R / pars[2];
+    const double r = norm3(q);
+    const double x = r / pars[2];
 
     // pi G rho r0^2 (pi - 2(1 - r0/r)arctan(r/r0) + 2(1 - r0/r)log(1 + r/r0) - (1 - r0/r)log(1 + (r/r0)^2))
     return -M_PI * pars[0] * pars[1] * pars[2] * pars[2] * (M_PI - 2 * (1 + 1 / x) * atan(x) + 2 * (1 + 1/x) * log(1 + x) - (1 - 1/x) * log(1 + x * x) );
@@ -2305,15 +2255,14 @@ void burkert_gradient_single(double t, double *__restrict__ pars, double6ptr q,
             - rho (mass scale)
             - r0
     */
-    double R, x, dphi_dr;
-    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    x = R / pars[2];
+    const double r = norm3(q);
+    const double x = r / pars[2];
 
-    dphi_dr = -M_PI * pars[0] * pars[1] * pars[2] / (x * x) * (2 * atan(x) - 2 * log(1 + x) - log(1 + x * x));
+    const double dphi_dr = -M_PI * pars[0] * pars[1] * pars[2] / (x * x) * (2 * atan(x) - 2 * log(1 + x) - log(1 + x * x));
 
-    grad[0] = grad[0] + dphi_dr*q[0]/R;
-    grad[1] = grad[1] + dphi_dr*q[1]/R;
-    grad[2] = grad[2] + dphi_dr*q[2]/R;
+    grad[0] = grad[0] + dphi_dr*q[0]/r;
+    grad[1] = grad[1] + dphi_dr*q[1]/r;
+    grad[2] = grad[2] + dphi_dr*q[2]/r;
 }
 
 
@@ -2323,13 +2272,9 @@ double burkert_density(double t, double *pars, double *q, int n_dim, void *state
             - rho (mass scale)
             - r0
     */
-    double r, x, rho;
-
-    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
-    x = r / pars[2];
-    rho = pars[1] / ((1 + x) * (1 + x * x));
-
-    return rho;
+    const double r = norm3(q);
+    const double x = r / pars[2];
+    return pars[1] / ((1 + x) * (1 + x * x));
 }
 
 DEFINE_VECTORIZED_GRADIENT(burkert)
diff --git a/src/gala/potential/potential/builtin/potential_helpers.h b/src/gala/potential/potential/builtin/potential_helpers.h
index 1e35cb67a..ac42f3923 100644
--- a/src/gala/potential/potential/builtin/potential_helpers.h
+++ b/src/gala/potential/potential/builtin/potential_helpers.h
@@ -1,9 +1,35 @@
 #include <math.h>
+#include "src/vectorization.h"
+
+static inline double norm2_sq(const double *q) {
+    return q[0]*q[0] + q[1]*q[1];
+}
 
 static inline double norm2(const double *q) {
-    return sqrt(q[0]*q[0] + q[1]*q[1]);
+    return sqrt(norm2_sq(q));
+}
+
+static inline double norm3_sq(const double *q) {
+    return q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
 }
 
 static inline double norm3(const double *q) {
-    return sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
+    return sqrt(norm3_sq(q));
+}
+
+static inline double norm3_flat_z(const double *q, const double qz) {
+    return sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]/(qz*qz));
+}
+
+// Helper functions for computing norms with double6ptr
+static inline double norm3_sq(const double6ptr& q) {
+    return (*q.x) * (*q.x) + (*q.y) * (*q.y) + (*q.z) * (*q.z);
+}
+
+static inline double norm3(const double6ptr& q) {
+    return sqrt(norm3_sq(q));
+}
+
+static inline double norm3_flat_z(const double6ptr& q, const double qz) {
+    return sqrt((*q.x) * (*q.x) + (*q.y) * (*q.y) + (*q.z) * (*q.z) / (qz * qz));
 }

From a86a5a19c0da71476922999c7f095ce27e504826 Mon Sep 17 00:00:00 2001
From: Adrian Price-Whelan <583379+adrn@users.noreply.github.com>
Date: Wed, 22 Oct 2025 19:04:00 -0400
Subject: [PATCH 4/6] i dunno - math doesn't seem to be working today

---
 .../potential/builtin/builtin_potentials.cpp  | 225 ++++++++++++++++++
 .../potential/builtin/builtin_potentials.h    |   8 +
 src/gala/potential/potential/builtin/core.py  |  53 +++++
 .../potential/potential/builtin/cybuiltin.pxd |   8 +
 .../potential/potential/builtin/cybuiltin.pyx |  23 ++
 tests/potential/potential/test_all_builtin.py |  58 +++++
 6 files changed, 375 insertions(+)

diff --git a/src/gala/potential/potential/builtin/builtin_potentials.cpp b/src/gala/potential/potential/builtin/builtin_potentials.cpp
index 6c3d71abe..8dec06b5e 100644
--- a/src/gala/potential/potential/builtin/builtin_potentials.cpp
+++ b/src/gala/potential/potential/builtin/builtin_potentials.cpp
@@ -2277,6 +2277,229 @@ double burkert_density(double t, double *pars, double *q, int n_dim, void *state
     return pars[1] / ((1 + x) * (1 + x * x));
 }
 
+#if USE_GSL == 1
+
+// --- Helper functions for Einasto profiles ---
+
+static inline double _s_of_r(double r, double r_m2, double alpha) {
+    // s(r) = (2/α) (r/r_-2 )^α - 1
+    if (r <= 0.0) return 0.0;
+    return (2.0 / alpha) * pow(r / r_m2, alpha);
+}
+
+static inline double gamma_beta(double beta, double x1, double x2) {
+    return gsl_sf_gamma_inc(beta, x1) - gsl_sf_gamma_inc(beta, x2);
+}
+
+static inline double gamma_tilde_beta(
+    double beta, double x1, double x2, double alpha, double r_s
+) {
+    const double scale = alpha * pow(r_s, alpha) / 2.0;
+    return pow(scale, beta) * gamma_beta(beta, x1, x2);
+}
+
+static inline double Gamma_beta(
+    double beta, double x
+) {
+    // Upper incomplete gamma:
+    return gsl_sf_gamma(beta) - gsl_sf_gamma_inc(beta, x);
+}
+
+static inline double Gamma_tilde_beta(
+    double beta, double x, double alpha, double r_s
+) {
+    const double scale = pow(alpha * pow(r_s, alpha) / 2.0, beta);
+    return scale * Gamma_beta(beta, x);
+}
+
+static inline double _einasto_mass_enclosed(double rho_m2, double r_m2, double alpha, double r) {
+    if (r <= 0.0) return 0.0;
+
+    const double sr  = _s_of_r(r, r_m2, alpha);
+    const double a3 = 3.0 / alpha;
+    const double h = r_m2 / pow(2.0 / alpha, 1.0 / alpha);
+    const double M = 4.0*M_PI * rho_m2 * pow(h, 3) / alpha * gsl_sf_gamma(a3);
+
+    return M * (1 - gsl_sf_gamma_inc(a3, sr) / gsl_sf_gamma(a3));
+
+}
+
+/* ---------------------------------------------------------------------------
+    Einasto profile
+    See: https://arxiv.org/abs/2004.10817
+*/
+double einasto_value(double t, double *pars, double *q, int n_dim, void *state) {
+    /*  pars:
+            - G (Gravitational constant)
+            - rho_-2 - scale density (density when log slope is -2)
+            - r_-2 - radius at which logarithmic slope of the density profile is equal to -2
+            - alpha - shape parameter
+    */
+    const double G = pars[0];
+    const double rho_m2 = pars[1];
+    const double r_m2 = pars[2];
+    const double alpha = pars[3];
+
+    const double r = norm3(q);
+
+    const double sr  = _s_of_r(r, r_m2, alpha);
+    const double a2 = 2.0 / alpha;
+    const double a3 = 3.0 / alpha;
+
+    const double h = r_m2 / pow(a2, 1.0 / alpha);
+    const double M = 4.0*M_PI * rho_m2 * pow(h, 3) / alpha * gsl_sf_gamma(a3);
+
+    if (r == 0.0) {
+        // Analytic limit at r = 0
+        return - G * M / h * gsl_sf_gamma(a2) / gsl_sf_gamma(a3);
+    }
+
+    const double factor = G * M / h / sr;
+    const double term1 = gsl_sf_gamma_inc(a3, pow(sr, alpha)) / gsl_sf_gamma(a3);
+    const double term2 = sr * gsl_sf_gamma_inc(a2, pow(sr, alpha)) / gsl_sf_gamma(a3);
+
+    return - factor * (1 - term1 + term2);
+
+}
+
+void einasto_gradient_single(double t, double *__restrict__ pars, double6ptr q, int n_dim, double6ptr grad, void *__restrict__ state) {
+    /*  pars:
+            - G (Gravitational constant)
+            - rho_-2 - scale density (density when log slope is -2)
+            - r_-2 - radius at which logarithmic slope of the density profile is equal to -2
+            - alpha - shape parameter
+    */
+    const double G = pars[0];
+    const double rho_m2 = pars[1];
+    const double r_m2 = pars[2];
+    const double alpha = pars[3];
+
+    const double r = norm3(&q[0]);
+    if (r == 0.0) {
+        return;
+    }
+
+    const double Menc = _einasto_mass_enclosed(rho_m2, r_m2, alpha, r);
+    const double dphi_dr = -(G * Menc) / (r * r * r);
+
+    grad[0] = grad[0] + dphi_dr * q[0];
+    grad[1] = grad[1] + dphi_dr * q[1];
+    grad[2] = grad[2] + dphi_dr * q[2];
+}
+
+
+double einasto_density(double t, double *pars, double *q, int n_dim, void *state) {
+    /*  pars:
+            - G (Gravitational constant)
+            - rho_-2 - scale density (density when log slope is -2)
+            - r_-2 - radius at which logarithmic slope of the density profile is equal to -2
+            - alpha - shape parameter
+    */
+    const double r = norm3(q);
+    const double alpha = pars[3];
+    return pars[1] * exp( - _s_of_r(r, pars[2], alpha) + (2.0 / alpha));
+}
+
+
+/* ---------------------------------------------------------------------------
+    Core Einasto potential
+    See: https://arxiv.org/abs/2004.10817
+*/
+
+static inline double _cEinasto_mass_enclosed(
+    double rho_s, double r_s, double alpha, double r_c, double r
+) {
+    // if (r <= 0.0) return 0.0;
+
+    // const double s0 = (2.0 / alpha) * pow(r_c / r_s, alpha);
+    // const double sr = (2.0 / alpha) * pow((r + r_c) / r_s, alpha);
+
+    // const double a1 = 1.0 / alpha;
+    // const double a2 = 2.0 / alpha;
+    // const double a3 = 3.0 / alpha;
+
+    // const double term1 =
+    // const double term2 = r_c * r_c * gamma_tilde_beta(a1, s0, sr, alpha, r_s);
+    // const double term3 = -2 * r_c * gamma_tilde_beta(a2, s0, sr, alpha, r_s);
+
+    // const double factor = (4.0 * M_PI * rho_s * exp(a2)) / alpha;
+    // return factor * ( term1 + term2 + term3 );
+    return 0.0;
+}
+
+double cEinasto_value(double t, double *pars, double *q, int n_dim, void *state) {
+    (void)t; (void)n_dim; (void)state;
+
+    const double G = pars[0];
+    const double rho_s = pars[1];
+    const double r_s = pars[2];
+    const double alpha = pars[3];
+    const double r_c = pars[4];
+
+    const double r = norm3(q);
+
+    const double a1 = 1.0 / alpha;
+    const double a2 = 2.0 / alpha;
+    const double a3 = 3.0 / alpha;
+
+    const double factor = (4.0 * M_PI * G * rho_s * exp(2.0/alpha)) / alpha;
+
+    const double s0 = _s_of_r(r_c, r_s, alpha);
+
+    if (r == 0.0) {
+        // Analytic limit at r = 0
+        const double term4 = Gamma_tilde_beta(a2, s0, alpha, r_s);
+        const double term5 = -r_c * Gamma_tilde_beta(a1, s0, alpha, r_s);
+        return factor * (term4 + term5);
+    }
+
+    const double sr = _s_of_r(r + r_c, r_s, alpha);
+
+    const double term1 = gamma_tilde_beta(a3, s0, sr, alpha, r_s);
+    const double term2 = r_c * r_c * gamma_tilde_beta(a1, s0, sr, alpha, r_s);
+    const double term3 = -2 * r_c * gamma_tilde_beta(a2, s0, sr, alpha, r_s);
+    const double term4 = Gamma_tilde_beta(a2, sr, alpha, r_s);
+    const double term5 = -r_c * Gamma_tilde_beta(a1, sr, alpha, r_s);
+
+    return factor * (1/r * (term1 + term2 + term3) + term4 + term5);
+}
+
+void cEinasto_gradient_single(
+    double t, double *__restrict__ pars, double6ptr q, int n_dim, double6ptr grad,
+    void *__restrict__ state
+) {
+    const double G = pars[0];
+    const double rho_s = pars[1];
+    const double r_s = pars[2];
+    const double alpha = pars[3];
+    const double r_c = pars[4];
+
+    const double r = norm3(&q[0]);
+    if (r == 0.0) {
+        return;
+    }
+
+    const double Menc = _cEinasto_mass_enclosed(rho_s, r_s, alpha, r_c, r);
+    const double dphi_dr = (G * Menc) / (r * r * r);
+
+    grad[0] = grad[0] + dphi_dr * q[0];
+    grad[1] = grad[1] + dphi_dr * q[1];
+    grad[2] = grad[2] + dphi_dr * q[2];
+}
+
+double cEinasto_density(double t, double *pars, double *q, int n_dim, void *state) {
+    const double rho_s = pars[1];
+    const double r_s = pars[2];
+    const double alpha = pars[3];
+    const double r_c = pars[4];
+
+    const double r = norm3(q);
+
+    return rho_s * exp( - 2.0 / alpha * (pow((r + r_c) / r_s, alpha) - 1) );
+}
+
+#endif
+
 DEFINE_VECTORIZED_GRADIENT(burkert)
 DEFINE_VECTORIZED_GRADIENT(flattenednfw)
 DEFINE_VECTORIZED_GRADIENT(henon_heiles)
@@ -2299,6 +2522,8 @@ DEFINE_VECTORIZED_GRADIENT(stone)
 DEFINE_VECTORIZED_GRADIENT(triaxialnfw)
 
 #if USE_GSL == 1
+DEFINE_VECTORIZED_GRADIENT(cEinasto)
+DEFINE_VECTORIZED_GRADIENT(einasto)
 DEFINE_VECTORIZED_GRADIENT(powerlawcutoff)
 DEFINE_VECTORIZED_GRADIENT(spherical_spline_density)
 DEFINE_VECTORIZED_GRADIENT(spherical_spline_mass)
diff --git a/src/gala/potential/potential/builtin/builtin_potentials.h b/src/gala/potential/potential/builtin/builtin_potentials.h
index 4d4a72342..d81bf65c4 100644
--- a/src/gala/potential/potential/builtin/builtin_potentials.h
+++ b/src/gala/potential/potential/builtin/builtin_potentials.h
@@ -145,6 +145,14 @@ extern double burkert_value(double t, double *pars, double *q, int n_dim, void *
 extern void burkert_gradient(double t, double *pars, double *q, int n_dim, size_t N, double *grad, void *state);
 extern double burkert_density(double t, double *pars, double *q, int n_dim, void *state);
 
+extern double einasto_value(double t, double *pars, double *q, int n_dim, void *state);
+extern void einasto_gradient(double t, double *pars, double *q, int n_dim, size_t N, double *grad, void *state);
+extern double einasto_density(double t, double *pars, double *q, int n_dim, void *state);
+
+extern double cEinasto_value(double t, double *pars, double *q, int n_dim, void *state);
+extern void cEinasto_gradient(double t, double *pars, double *q, int n_dim, size_t N, double *grad, void *state);
+extern double cEinasto_density(double t, double *pars, double *q, int n_dim, void *state);
+
 // Spherical spline interpolated potentials
 extern double spherical_spline_density_value(double t, double *pars, double *q, int n_dim, void *state);
 extern void spherical_spline_density_gradient(double t, double *pars, double *q, int n_dim, size_t N, double *grad, void *state);
diff --git a/src/gala/potential/potential/builtin/core.py b/src/gala/potential/potential/builtin/core.py
index 8ae610a57..956a9506f 100644
--- a/src/gala/potential/potential/builtin/core.py
+++ b/src/gala/potential/potential/builtin/core.py
@@ -18,7 +18,9 @@
 from gala.potential.common import PotentialParameter
 from gala.potential.potential.builtin.cybuiltin import (
     BurkertWrapper,
+    CoreEinastoWrapper,
     CylSplineWrapper,
+    EinastoWrapper,
     FlattenedNFWWrapper,
     HenonHeilesWrapper,
     HernquistWrapper,
@@ -53,8 +55,10 @@
 
 __all__ = [
     "BurkertPotential",
+    "CoreEinastoPotential",
     "CylSplinePotential",
     "EXPPotential",
+    "EinastoPotential",
     "HenonHeilesPotential",
     "HernquistPotential",
     "IsochronePotential",
@@ -449,6 +453,55 @@ def from_r0(cls, r0, units=None):
         return cls(rho=rho_d0, r0=r0, units=units)
 
 
+@format_doc(common_doc=_potential_docstring)
+class EinastoPotential(CPotentialBase):
+    r"""The Einasto model.
+
+    Parameters
+    ----------
+    rho_m2 : :class:`~astropy.units.Quantity`, numeric [mass density]
+        Density at the radius where the logarithmic slope of the density profile is -2.
+    r_m2 : :class:`~astropy.units.Quantity`, numeric [length]
+        Radius where the logarithmic slope of the density profile is -2.
+    alpha :
+        Shape parameter.
+    {common_doc}
+    """
+
+    rho_m2 = PotentialParameter("rho_m2", physical_type="mass density")
+    r_m2 = PotentialParameter("r_m2", physical_type="length")
+    alpha = PotentialParameter("alpha", physical_type="dimensionless")
+
+    Wrapper = EinastoWrapper
+    _symmetry = SphericalSymmetry()
+
+
+@format_doc(common_doc=_potential_docstring)
+class CoreEinastoPotential(CPotentialBase):
+    r"""The core Einasto model.
+
+    Parameters
+    ----------
+    rho_s : :class:`~astropy.units.Quantity`, numeric [mass density]
+        Scale density (analogous to rho_-2 in the standard Einasto profile).
+    r_s : :class:`~astropy.units.Quantity`, numeric [length]
+        Scale radius (analogous to r_-2 in the standard Einasto profile).
+    alpha :
+        Shape parameter.
+    r_c : :class:`~astropy.units.Quantity`, numeric [length]
+        The core radius.
+    {common_doc}
+    """
+
+    rho_s = PotentialParameter("rho_s", physical_type="mass density")
+    r_s = PotentialParameter("r_s", physical_type="length")
+    alpha = PotentialParameter("alpha", physical_type="dimensionless")
+    r_c = PotentialParameter("r_c", physical_type="length")
+
+    Wrapper = CoreEinastoWrapper
+    _symmetry = SphericalSymmetry()
+
+
 # ============================================================================
 # Flattened, axisymmetric models
 #
diff --git a/src/gala/potential/potential/builtin/cybuiltin.pxd b/src/gala/potential/potential/builtin/cybuiltin.pxd
index a6953fade..bdbf93fce 100644
--- a/src/gala/potential/potential/builtin/cybuiltin.pxd
+++ b/src/gala/potential/potential/builtin/cybuiltin.pxd
@@ -100,3 +100,11 @@ cdef extern from "potential/potential/builtin/builtin_potentials.h":
     double burkert_value(double t, double *pars, double *q, int n_dim, void *state) nogil
     void burkert_gradient(double t, double *pars, double *q, int n_dim, size_t N, double *grad, void *state) nogil
     double burkert_density(double t, double *pars, double *q, int n_dim, void *state) nogil
+
+    double einasto_value(double t, double *pars, double *q, int n_dim, void *state) nogil
+    void einasto_gradient(double t, double *pars, double *q, int n_dim, size_t N, double *grad, void *state) nogil
+    double einasto_density(double t, double *pars, double *q, int n_dim, void *state) nogil
+
+    double cEinasto_value(double t, double *pars, double *q, int n_dim, void *state) nogil
+    void cEinasto_gradient(double t, double *pars, double *q, int n_dim, size_t N, double *grad, void *state) nogil
+    double cEinasto_density(double t, double *pars, double *q, int n_dim, void *state) nogil
diff --git a/src/gala/potential/potential/builtin/cybuiltin.pyx b/src/gala/potential/potential/builtin/cybuiltin.pyx
index 64b4a7f9a..a5be783de 100644
--- a/src/gala/potential/potential/builtin/cybuiltin.pyx
+++ b/src/gala/potential/potential/builtin/cybuiltin.pyx
@@ -223,6 +223,7 @@ cdef class PowerLawCutoffWrapper(CPotentialWrapper):
             self.cpotential.gradient[0] = <gradientfunc>(powerlawcutoff_gradient)
             self.cpotential.hessian[0] = <hessianfunc>(powerlawcutoff_hessian)
 
+
 cdef class BurkertWrapper(CPotentialWrapper):
 
     def __init__(self, G, parameters, q0, R):
@@ -234,6 +235,28 @@ cdef class BurkertWrapper(CPotentialWrapper):
         self.cpotential.gradient[0] = <gradientfunc>(burkert_gradient)
 
 
+cdef class EinastoWrapper(CPotentialWrapper):
+
+    def __init__(self, G, parameters, q0, R):
+        self.init([G] + list(parameters),
+                  np.ascontiguousarray(q0),
+                  np.ascontiguousarray(R))
+        self.cpotential.value[0] = <energyfunc>(einasto_value)
+        self.cpotential.density[0] = <densityfunc>(einasto_density)
+        self.cpotential.gradient[0] = <gradientfunc>(einasto_gradient)
+
+
+cdef class CoreEinastoWrapper(CPotentialWrapper):
+
+    def __init__(self, G, parameters, q0, R):
+        self.init([G] + list(parameters),
+                  np.ascontiguousarray(q0),
+                  np.ascontiguousarray(R))
+        self.cpotential.value[0] = <energyfunc>(cEinasto_value)
+        self.cpotential.density[0] = <densityfunc>(cEinasto_density)
+        self.cpotential.gradient[0] = <gradientfunc>(cEinasto_gradient)
+
+
 # ============================================================================
 # Flattened, axisymmetric models
 #
diff --git a/tests/potential/potential/test_all_builtin.py b/tests/potential/potential/test_all_builtin.py
index 9cc7177cf..c24d3de2e 100644
--- a/tests/potential/potential/test_all_builtin.py
+++ b/tests/potential/potential/test_all_builtin.py
@@ -624,3 +624,61 @@ def test_from_r0(self):
 
         # Check a 1% tolerance on inferred density against published values
         assert abs(rho - rho_check) / rho_check < 0.01
+
+
+@pytest.mark.skipif(not GSL_ENABLED, reason="requires GSL to run this test")
+class TestEinasto(PotentialTestBase):
+    potential = p.EinastoPotential(
+        units=galactic, rho_m2=5e-25 * u.g / u.cm**3, r_m2=12 * u.kpc, alpha=0.16
+    )
+    w0 = [1.0, 0.0, 0.0, 0.0, 0.1, 0.1]
+
+    check_finite_at_origin = True
+    skip_hessian = True  # TODO: implement
+
+    @pytest.mark.skip(reason="Not implemented for Burkert potentials")
+    def test_against_sympy(self):
+        pass
+
+    # def test_from_r0(self):
+    #     # TODO: test alternate constructor
+
+    #     # Test against values from Zhu+2023
+    #     pot = p.BurkertPotential.from_r0(r0=11.87 * u.kpc, units=galactic)
+
+    #     rho = pot.parameters["rho"].to(u.g / u.cm**3)
+    #     rho_check = 5.93e-25 * u.g / u.cm**3
+
+    #     # Check a 1% tolerance on inferred density against published values
+    #     assert abs(rho - rho_check) / rho_check < 0.01
+
+
+@pytest.mark.skipif(not GSL_ENABLED, reason="requires GSL to run this test")
+class TestCoreEinasto(PotentialTestBase):
+    potential = p.CoreEinastoPotential(
+        units=galactic,
+        rho_s=5e-25 * u.g / u.cm**3,
+        r_s=12 * u.kpc,
+        alpha=0.16,
+        r_c=1 * u.kpc,
+    )
+    w0 = [1.0, 0.0, 0.0, 0.0, 0.1, 0.1]
+    skip_hessian = True  # TODO: implement
+
+    check_finite_at_origin = True
+
+    @pytest.mark.skip(reason="Not implemented for Burkert potentials")
+    def test_against_sympy(self):
+        pass
+
+    # def test_from_r0(self):
+    #     # TODO: test alternate constructor
+
+    #     # Test against values from Zhu+2023
+    #     pot = p.BurkertPotential.from_r0(r0=11.87 * u.kpc, units=galactic)
+
+    #     rho = pot.parameters["rho"].to(u.g / u.cm**3)
+    #     rho_check = 5.93e-25 * u.g / u.cm**3
+
+    #     # Check a 1% tolerance on inferred density against published values
+    #     assert abs(rho - rho_check) / rho_check < 0.01

From 9e9bdef2ce0d9115582c95b04c9a4685ad843982 Mon Sep 17 00:00:00 2001
From: Adrian Price-Whelan <583379+adrn@users.noreply.github.com>
Date: Wed, 22 Oct 2025 20:57:15 -0400
Subject: [PATCH 5/6] some fixes in einasto definitions

---
 .../potential/builtin/builtin_potentials.cpp  | 105 ++++++++++--------
 1 file changed, 61 insertions(+), 44 deletions(-)

diff --git a/src/gala/potential/potential/builtin/builtin_potentials.cpp b/src/gala/potential/potential/builtin/builtin_potentials.cpp
index 8dec06b5e..3571a8103 100644
--- a/src/gala/potential/potential/builtin/builtin_potentials.cpp
+++ b/src/gala/potential/potential/builtin/builtin_potentials.cpp
@@ -2282,13 +2282,13 @@ double burkert_density(double t, double *pars, double *q, int n_dim, void *state
 // --- Helper functions for Einasto profiles ---
 
 static inline double _s_of_r(double r, double r_m2, double alpha) {
-    // s(r) = (2/α) (r/r_-2 )^α - 1
+    // s(r) = (2/α) (r/r_-2 )^α
     if (r <= 0.0) return 0.0;
     return (2.0 / alpha) * pow(r / r_m2, alpha);
 }
 
 static inline double gamma_beta(double beta, double x1, double x2) {
-    return gsl_sf_gamma_inc(beta, x1) - gsl_sf_gamma_inc(beta, x2);
+    return gsl_sf_gamma_inc(beta, x2) - gsl_sf_gamma_inc(beta, x1);
 }
 
 static inline double gamma_tilde_beta(
@@ -2315,13 +2315,12 @@ static inline double Gamma_tilde_beta(
 static inline double _einasto_mass_enclosed(double rho_m2, double r_m2, double alpha, double r) {
     if (r <= 0.0) return 0.0;
 
-    const double sr  = _s_of_r(r, r_m2, alpha);
+    const double s = _s_of_r(r, r_m2, alpha);
     const double a3 = 3.0 / alpha;
-    const double h = r_m2 / pow(2.0 / alpha, 1.0 / alpha);
-    const double M = 4.0*M_PI * rho_m2 * pow(h, 3) / alpha * gsl_sf_gamma(a3);
-
-    return M * (1 - gsl_sf_gamma_inc(a3, sr) / gsl_sf_gamma(a3));
 
+    const double Mtot = (4.0 * M_PI * rho_m2 * exp(2.0/alpha) / alpha)
+                        * pow(r_m2, 3.0) * pow(alpha/2.0, a3) * gsl_sf_gamma(a3);
+    return Mtot * gsl_sf_gamma_inc_P(a3, s);
 }
 
 /* ---------------------------------------------------------------------------
@@ -2346,19 +2345,17 @@ double einasto_value(double t, double *pars, double *q, int n_dim, void *state)
     const double a2 = 2.0 / alpha;
     const double a3 = 3.0 / alpha;
 
-    const double h = r_m2 / pow(a2, 1.0 / alpha);
-    const double M = 4.0*M_PI * rho_m2 * pow(h, 3) / alpha * gsl_sf_gamma(a3);
+    const double pref  = (4.0 * M_PI * G * rho_m2 * exp(2.0/alpha)) / alpha;
+    const double scale = alpha * pow(r_m2, alpha) / 2.0;
 
     if (r == 0.0) {
-        // Analytic limit at r = 0
-        return - G * M / h * gsl_sf_gamma(a2) / gsl_sf_gamma(a3);
+        const double term1 = (pow(r_m2, 3.0) / 3.0) * pow(2.0/alpha, a3);
+        const double term2 = pow(scale, a2) * gsl_sf_gamma(a2);
+        return - pref * (term1 + term2);
     }
-
-    const double factor = G * M / h / sr;
-    const double term1 = gsl_sf_gamma_inc(a3, pow(sr, alpha)) / gsl_sf_gamma(a3);
-    const double term2 = sr * gsl_sf_gamma_inc(a2, pow(sr, alpha)) / gsl_sf_gamma(a3);
-
-    return - factor * (1 - term1 + term2);
+    const double term1 = pow(scale, a3) * gsl_sf_gamma(a3) * gsl_sf_gamma_inc_P(a3, sr);
+    const double term2 = pow(scale, a2) * gsl_sf_gamma(a2) * gsl_sf_gamma_inc_Q(a2, sr);
+    return - pref * ( (term1 / r) + term2 );
 
 }
 
@@ -2374,13 +2371,13 @@ void einasto_gradient_single(double t, double *__restrict__ pars, double6ptr q,
     const double r_m2 = pars[2];
     const double alpha = pars[3];
 
-    const double r = norm3(&q[0]);
+    const double r = norm3(q);
     if (r == 0.0) {
         return;
     }
 
     const double Menc = _einasto_mass_enclosed(rho_m2, r_m2, alpha, r);
-    const double dphi_dr = -(G * Menc) / (r * r * r);
+    const double dphi_dr = (G * Menc) / (r * r * r);
 
     grad[0] = grad[0] + dphi_dr * q[0];
     grad[1] = grad[1] + dphi_dr * q[1];
@@ -2409,27 +2406,34 @@ double einasto_density(double t, double *pars, double *q, int n_dim, void *state
 static inline double _cEinasto_mass_enclosed(
     double rho_s, double r_s, double alpha, double r_c, double r
 ) {
-    // if (r <= 0.0) return 0.0;
+    if (r <= 0.0) return 0.0;
+
+    const double s0 = _s_of_r(r_c, r_s, alpha);
+    const double sr = _s_of_r(r + r_c, r_s, alpha);
+
+    const double a1 = 1.0 / alpha;
+    const double a2 = 2.0 / alpha;
+    const double a3 = 3.0 / alpha;
 
-    // const double s0 = (2.0 / alpha) * pow(r_c / r_s, alpha);
-    // const double sr = (2.0 / alpha) * pow((r + r_c) / r_s, alpha);
+    const double scale = alpha * pow(r_s, alpha) / 2.0;
 
-    // const double a1 = 1.0 / alpha;
-    // const double a2 = 2.0 / alpha;
-    // const double a3 = 3.0 / alpha;
+    const double seg1 = pow(scale, a1) * (gsl_sf_gamma_inc(a1, sr) - gsl_sf_gamma_inc(a1, s0));
+    const double seg2 = pow(scale, a2) * (gsl_sf_gamma_inc(a2, sr) - gsl_sf_gamma_inc(a2, s0));
+    const double seg3 = pow(scale, a3) * (gsl_sf_gamma_inc(a3, sr) - gsl_sf_gamma_inc(a3, s0));
 
-    // const double term1 =
-    // const double term2 = r_c * r_c * gamma_tilde_beta(a1, s0, sr, alpha, r_s);
-    // const double term3 = -2 * r_c * gamma_tilde_beta(a2, s0, sr, alpha, r_s);
+    const double pref = (4.0 * M_PI * rho_s * exp(2.0/alpha)) / alpha;
 
-    // const double factor = (4.0 * M_PI * rho_s * exp(a2)) / alpha;
-    // return factor * ( term1 + term2 + term3 );
-    return 0.0;
+    return pref * ( seg3 + (r_c*r_c)*seg1 - 2.0*r_c*seg2 );
 }
 
 double cEinasto_value(double t, double *pars, double *q, int n_dim, void *state) {
-    (void)t; (void)n_dim; (void)state;
-
+    /*  pars:
+            - G (Gravitational constant)
+            - rho_s - scale density
+            - r_s - scale radius
+            - alpha - shape parameter
+            - r_c - core radius
+    */
     const double G = pars[0];
     const double rho_s = pars[1];
     const double r_s = pars[2];
@@ -2442,26 +2446,39 @@ double cEinasto_value(double t, double *pars, double *q, int n_dim, void *state)
     const double a2 = 2.0 / alpha;
     const double a3 = 3.0 / alpha;
 
-    const double factor = (4.0 * M_PI * G * rho_s * exp(2.0/alpha)) / alpha;
+    const double pref  = (4.0 * M_PI * G * rho_s * exp(2.0/alpha)) / alpha;
+    const double scale = alpha * pow(r_s, alpha) / 2.0;
 
     const double s0 = _s_of_r(r_c, r_s, alpha);
 
     if (r == 0.0) {
-        // Analytic limit at r = 0
-        const double term4 = Gamma_tilde_beta(a2, s0, alpha, r_s);
-        const double term5 = -r_c * Gamma_tilde_beta(a1, s0, alpha, r_s);
-        return factor * (term4 + term5);
+        const double dsdr0 = 2.0 * pow(r_c, alpha - 1.0) / pow(r_s, alpha);
+        const double e_m_s0 = exp(-s0);
+
+        // (1/r)*[γ̃_{3/α} + r_c^2 γ̃_{1/α} - 2 r_c γ̃_{2/α}]  -->  ds/dr|0 * e^{-s0} * sum
+        const double lim_over_r =
+            dsdr0 * e_m_s0 *
+            ( pow(scale, a3) * pow(s0, a3 - 1.0)
+            + (r_c*r_c) * pow(scale, a1) * pow(s0, a1 - 1.0)
+            - 2.0 * r_c * pow(scale, a2) * pow(s0, a2 - 1.0) );
+
+        // Upper tilded terms at r=0: Γ̃_b(s0) = scale^b * Γ(b) * Q(b, s0)
+        const double up2 = pow(scale, a2) * gsl_sf_gamma(a2) * gsl_sf_gamma_inc_Q(a2, s0);
+        const double up1 = pow(scale, a1) * gsl_sf_gamma(a1) * gsl_sf_gamma_inc_Q(a1, s0);
+
+        return -pref * ( lim_over_r + (up2 - r_c * up1) );
     }
 
     const double sr = _s_of_r(r + r_c, r_s, alpha);
 
-    const double term1 = gamma_tilde_beta(a3, s0, sr, alpha, r_s);
-    const double term2 = r_c * r_c * gamma_tilde_beta(a1, s0, sr, alpha, r_s);
-    const double term3 = -2 * r_c * gamma_tilde_beta(a2, s0, sr, alpha, r_s);
-    const double term4 = Gamma_tilde_beta(a2, sr, alpha, r_s);
-    const double term5 = -r_c * Gamma_tilde_beta(a1, sr, alpha, r_s);
+    const double seg1 = pow(scale, a1) * (gsl_sf_gamma_inc(a1, sr) - gsl_sf_gamma_inc(a1, s0));
+    const double seg2 = pow(scale, a2) * (gsl_sf_gamma_inc(a2, sr) - gsl_sf_gamma_inc(a2, s0));
+    const double seg3 = pow(scale, a3) * (gsl_sf_gamma_inc(a3, sr) - gsl_sf_gamma_inc(a3, s0));
+
+    const double up2 = pow(scale, a2) * gsl_sf_gamma(a2) * gsl_sf_gamma_inc_Q(a2, sr);
+    const double up1 = pow(scale, a1) * gsl_sf_gamma(a1) * gsl_sf_gamma_inc_Q(a1, sr);
 
-    return factor * (1/r * (term1 + term2 + term3) + term4 + term5);
+    return -pref * ( ((seg3 + (r_c*r_c)*seg1 - 2.0*r_c*seg2) / r) + (up2 - r_c * up1) );
 }
 
 void cEinasto_gradient_single(

From 9eed01e981774821aaeb3f2fbffab03f6258940d Mon Sep 17 00:00:00 2001
From: Adrian Price-Whelan <583379+adrn@users.noreply.github.com>
Date: Thu, 23 Oct 2025 08:12:33 -0400
Subject: [PATCH 6/6] minor cEinasto tweaks

---
 src/gala/potential/potential/builtin/builtin_potentials.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/gala/potential/potential/builtin/builtin_potentials.cpp b/src/gala/potential/potential/builtin/builtin_potentials.cpp
index 3571a8103..799ec7d94 100644
--- a/src/gala/potential/potential/builtin/builtin_potentials.cpp
+++ b/src/gala/potential/potential/builtin/builtin_potentials.cpp
@@ -2491,7 +2491,7 @@ void cEinasto_gradient_single(
     const double alpha = pars[3];
     const double r_c = pars[4];
 
-    const double r = norm3(&q[0]);
+    const double r = norm3(q);
     if (r == 0.0) {
         return;
     }
@@ -2512,7 +2512,7 @@ double cEinasto_density(double t, double *pars, double *q, int n_dim, void *stat
 
     const double r = norm3(q);
 
-    return rho_s * exp( - 2.0 / alpha * (pow((r + r_c) / r_s, alpha) - 1) );
+    return rho_s * exp( - _s_of_r(r + r_c, r_s, alpha) + (2.0 / alpha) );
 }
 
 #endif