Add Einasto potential #760

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.

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

jobovy merged 4 commits into jobovy:main from johnwez1:potential/einasto

Nov 3, 2025

.pre-commit-config.yaml

-Original file line number
+Diff line change
@@ Expand Up / @@ -33,7 +33,7 @@ repos: @@
       rev: "v2.4.1"
       hooks:
         - id: codespell
-          args: ["-L", "thisE,thise,mye,tE,te,hist,ro,sav,ccompiler,aas,floatIn,dOmin,indx",
+          args: ["-L", "thisE,thise,mye,tE,te,hist,ro,sav,ccompiler,aas,floatIn,dOmin,indx,Hese",
           "-x","doc/source/_static/try-galpy.js"]
     - repo: https://github.com/astral-sh/ruff-pre-commit
       rev: v0.13.3
@@ Expand Down @@

doc/source/reference/potential.rst

-Original file line number
+Diff line change
@@ Expand Up / @@ -171,6 +171,7 @@ ellipsoidal triaxial potentials below. @@
        potentialdoublepowerspher.rst
        potentialcoredehnen.rst
        potentialdehnen.rst
+       potentialeinasto.rst
        potentialhernquist.rst
        potentialhomogsphere.rst
        potentialinterpsphere.rst
@@ Expand Down @@

doc/source/reference/potentialeinasto.rst

-Original file line number
+Diff line change
@@ -0,0 +1,5 @@
+    Einasto potential
+    ======================
+    .. autoclass:: galpy.potential.EinastoPotential
+       :members: __init__

galpy/orbit/integrateFullOrbit.py

-Original file line number
+Diff line change
@@ Expand Up @@
                         )
             elif isinstance(p, potential.NullPotential):
                 pot_type.append(40)
+            elif isinstance(p, potential.EinastoPotential):
+                pot_type.append(41)
+                pot_args.extend([p._amp, p.h, p.n])
                 # No arguments, zero forces
             ############################## WRAPPERS ###############################
             elif isinstance(p, potential.DehnenSmoothWrapperPotential):
@@ Expand Down @@

galpy/orbit/integratePlanarOrbit.py

-Original file line number
+Diff line change
@@ Expand Up / @@ -466,6 +466,11 @@ def _parse_pot(pot): @@
                 p._Pot, potential.NullPotential
             ):
                 pot_type.append(40)
+            elif isinstance(p, planarPotentialFromRZPotential) and isinstance(
+                p._Pot, potential.EinastoPotential
+            ):
+                pot_type.append(41)
+                pot_args.extend([p._Pot._amp, p._Pot.h, p._Pot.n])
             ############################## WRAPPERS ###############################
             elif (
                 (
@@ Expand Down @@

galpy/orbit/orbit_c_ext/integrateFullOrbit.c

-Original file line number
+Diff line change
@@ Expand Up / @@ -530,6 +530,21 @@ void parse_leapFuncArgs_Full(int npot, @@
           potentialArgs->ntfuncs= 0;
           potentialArgs->requiresVelocity= false;
           break;
+        case 41: //EinastoPotential
+          potentialArgs->potentialEval= &SphericalPotentialEval;
+          potentialArgs->Rforce = &SphericalPotentialRforce;
+          potentialArgs->zforce = &SphericalPotentialzforce;
+          potentialArgs->phitorque= &ZeroForce;
+          potentialArgs->dens= &SphericalPotentialDens;
+          // Also assign functions specific to SphericalPotential
+          potentialArgs->revaluate= &EinastoPotentialrevaluate;
+          potentialArgs->rforce= &EinastoPotentialrforce;
+          potentialArgs->r2deriv= &EinastoPotentialr2deriv;
+          potentialArgs->rdens= &EinastoPotentialrdens;
+          potentialArgs->nargs = 3;
+          potentialArgs->ntfuncs= 0;
+          potentialArgs->requiresVelocity= false;
+          break;
     //////////////////////////////// WRAPPERS /////////////////////////////////////
         case -1: //DehnenSmoothWrapperPotential
           potentialArgs->potentialEval= &DehnenSmoothWrapperPotentialEval;
@@ Expand Down @@

galpy/orbit/orbit_c_ext/integratePlanarOrbit.c

-Original file line number
+Diff line change
@@ Expand Up @@
           potentialArgs->ntfuncs= 0;
           potentialArgs->requiresVelocity= false;
           break;
+        case 41: //EinastoPotential
+          potentialArgs->potentialEval= &SphericalPotentialEval;
+          potentialArgs->planarRforce = &SphericalPotentialPlanarRforce;
+          potentialArgs->planarphitorque= &ZeroPlanarForce;
+          potentialArgs->planarR2deriv= &SphericalPotentialPlanarR2deriv;
+          potentialArgs->planarphi2deriv= &ZeroPlanarForce;
+          potentialArgs->planarRphideriv= &ZeroPlanarForce;
+          // Also assign functions specific to SphericalPotential
+          potentialArgs->revaluate= &EinastoPotentialrevaluate;
+          potentialArgs->rforce= &EinastoPotentialrforce;
+          potentialArgs->r2deriv= &EinastoPotentialr2deriv;
+          potentialArgs->nargs = 3;
+          potentialArgs->ntfuncs= 0;
+          potentialArgs->requiresVelocity= false;
+          break;
     //////////////////////////////// WRAPPERS /////////////////////////////////////
         case -1: //DehnenSmoothWrapperPotential
           potentialArgs->potentialEval= &DehnenSmoothWrapperPotentialEval;
@@ Expand Down @@

galpy/potential/EinastoPotential.py

-Original file line number
+Diff line change
@@ -0,0 +1,156 @@
+    ###############################################################################
+    #   BurkertPotential.py: Potential with a Burkert density
+    ###############################################################################
+    import numpy
+    from scipy import special
+    from scipy.optimize import fsolve
+    from ..util import conversion
+    from .SphericalPotential import SphericalPotential
+    class EinastoPotential(SphericalPotential):
+        """Potential with an Einasto [1]_ density. Class implements the following interchangeable conventions:
+        .. math::
+            \\rho(r) = \\mathrm{amp}\\,\\exp\\left(-d_n\\left[\\left(\\frac{r}{r_s}\\right)^\\frac{1}{n}-1\\right]\\right)
+        or
+        .. math::
+            \\rho(r) = \\mathrm{amp}\\,\\exp\\left(-2n\\left[\\left(\\frac{r}{r_{-2}}\\right)^\\frac{1}{n}-1\\right]\\right)
+        or
+        .. math::
+            \\rho(r) = \\mathrm{amp}\\,\\exp\\left(-\\left(\\frac{r}{h}\\right)^\\frac{1}{n}\\right)
+        With conventions taken from [2]_.
+        """
+        def __init__(
+            self, amp=1.0, h=2.0, n=1, rs=None, rm2=None, normalize=False, ro=None, vo=None
+        ):
+            """
+            Initialize a Einasto-density potential [1]_.
+            Parameters
+            ----------
+            amp : float or Quantity
+                Amplitude to be applied to the potential. Can be a Quantity with units of mass density or Gxmass density.
+            h : float or Quantity
+                Scale length.
+            rs : float or Quantity
+                Radius of the sphere that contains half of the total mass.
+            rm2 : float or Quantity
+                Radius at which rho(r) ∝ r^-2.
+            n : float
+                The Einasto index. A shape parameter defining the steepness of the power law
+            normalize : bool or float, optional
+                If True, normalize such that vc(1.,0.)=1., or, if given as a number, such that the force is this fraction of the force necessary to make vc(1.,0.)=1. Default is False.
+            ro : float or Quantity, optional
+                Distance scale for translation into internal units (default from configuration file).
+            vo : float or Quantity, optional
+                Velocity scale for translation into internal units (default from configuration file).
+            Notes
+            -----
+            - Either specify h or rs or rm2.
+            - 2025-09-12 - Written - John Weatherall
+            References
+            ----------
+            .. [1] Einasto (1965), Trudy Inst. Astroz. Alma-Ata, No. 17, 1 ADS: https://ui.adsabs.harvard.edu/abs/1965TrAlm...5...87E.
+            .. [2] Retana-Montenegro, E., Van Hese, E., Gentile, G., Baes, M., & Frutos-Alfaro, F. 2012, A&A, 540, A70 ADS: https://ui.adsabs.harvard.edu/abs/2012A&A...540A..70R.
+            """
+            SphericalPotential.__init__(self, amp=amp, ro=ro, vo=vo, amp_units="density")
+            if rs is not None:
+                rs = conversion.parse_length(rs, ro=self._ro, vo=self._vo)
+                # convert to h
+                dn = self._estimate_dn(n)
+                dn = self._calculate_dn(n, dn)
+                self.amp = amp * numpy.e**dn
+                h = rs / (dn**n)
+            elif rm2 is not None:
+                rm2 = conversion.parse_length(rm2, ro=self._ro, vo=self._vo)
+                # convert to h
+                self.amp = amp * numpy.e ** (2 * n)
+                h = rm2 / ((2 * n) ** n)
+            else:
+                h = conversion.parse_length(h, ro=self._ro, vo=self._vo)
+            self.h = h
+            self.n = n
+            self._scale = self.h
+            if normalize or (
+                isinstance(normalize, (int, float)) and not isinstance(normalize, bool)
+            ):  # pragma: no cover
+                self.normalize(normalize)
+            self.hasC = True
+            self.hasC_dxdv = True
+            self.hasC_dens = True
+            return None
+        def _revaluate(self, r, t=0.0):
+            """Potential as a function of r and time"""
+            s = r / self.h
+            gamma_3n = special.gamma(3 * self.n)
+            gamma_2n = special.gamma(2 * self.n)
+            gamma_upper_3n = special.gammaincc(3 * self.n, (s ** (1 / self.n)))
+            gamma_upper_2n = special.gammaincc(2 * self.n, (s ** (1 / self.n)))
+            # written to handle s = numpy.inf
+            out = -(4 * numpy.pi * (self.h**2) * self.n * gamma_3n) * (
+                (1 - gamma_upper_3n) / s + gamma_upper_2n * (gamma_2n / gamma_3n)
+            )
+            core = -(4 * numpy.pi * (self.h**2) * self.n) * special.gamma(2 * self.n)
+            if isinstance(r, (float, int)):
+                if r == 0:
+                    return core
+                else:
+                    return out
+            else:
+                out[r == 0] = core
+                return out
+        def _rforce(self, r, t=0.0):
+            s = r / self.h
+            gamma_3n = special.gamma(3 * self.n)
+            gamma_upper_3n = special.gammaincc(3 * self.n, (s ** (1 / self.n)))
+            return (
+                (4 * numpy.pi * self.h * self.n * gamma_3n) * (s**-2) * (gamma_upper_3n - 1)
+            )
+        def _r2deriv(self, r, t=0.0):
+            s = r / self.h
+            gamma_3n = special.gamma(3 * self.n)
+            gamma_upper_3n = special.gammaincc(3 * self.n, (s ** (1 / self.n)))
+            # (self.h**2)
+            return -(4 * numpy.pi * self.n * gamma_3n) * (
+                (-2 * (s**-3)) * (gamma_upper_3n - 1)
+                - ((1 / self.n) * (numpy.e ** -(s ** (1 / self.n))) / gamma_3n)
+            )
+        def _rdens(self, r, t=0.0):
+            return numpy.e ** -((r / self.h) ** (1 / self.n))
+        def _estimate_dn(self, n):
+            # see [2]
+            return (
+* n
+                - 1 / 3
+                + (8 / (1215 * n))
+                + (184 / (229635 * n**2))
+                + (1048 / (31000725 * n**3))
+                - (17557576 / (1242974068875 * n**4))
+            )
+        def _calculate_dn(self, n, est_dn):
+            # use numerical solver
+            def func(x):
+                gamma_3n = special.gamma(3 * n)
+                gamma_3n_upper = special.gammaincc(3 * n, x) * gamma_3n
+                return 2 * gamma_3n_upper - gamma_3n
+            return fsolve(func, est_dn)[0]

galpy/potential/__init__.py

-Original file line number
+Diff line change
@@ Expand Up / @@ -10,6 +10,7 @@ @@
         DehnenSmoothWrapperPotential,
         DiskSCFPotential,
         DoubleExponentialDiskPotential,
+        EinastoPotential,
         EllipticalDiskPotential,
         FDMDynamicalFrictionForce,
         FerrersPotential,
@@ Expand Down Expand Up / @@ -158,6 +159,7 @@ @@
     )
     DehnenSphericalPotential = TwoPowerSphericalPotential.DehnenSphericalPotential
     DehnenCoreSphericalPotential = TwoPowerSphericalPotential.DehnenCoreSphericalPotential
+    EinastoPotential = EinastoPotential.EinastoPotential
     NFWPotential = TwoPowerSphericalPotential.NFWPotential
     JaffePotential = TwoPowerSphericalPotential.JaffePotential
     HernquistPotential = TwoPowerSphericalPotential.HernquistPotential
@@ Expand Down @@

galpy/potential/potential_c_ext/EinastoPotential.c

-Original file line number
+Diff line change
@@ -0,0 +1,83 @@
+    #include <math.h>
+    #include <gsl/gsl_sf_gamma.h>
+    #include <galpy_potentials.h>
+    #ifndef M_PI
+    #define M_PI 3.14159265358979323846
+    #endif
+    //EinastoPotential
+    // 3 arguments: amp, h, n
+    double EinastoPotentialrevaluate(double r,double t,
+      struct potentialArg * potentialArgs){
+      double * args= potentialArgs->args;
+      //Get args
+      double amp= *args;
+      double h= *(args+1);
+      double n= *(args+2);
+      double s = r / h;
+      double s_1n = pow(s,1.0/n);
+      return (
+          -(4.0 * M_PI * pow(h,2) * n * gsl_sf_gamma(3*n))
+          * pow(s,-1)
+          * (
+- gsl_sf_gamma_inc_Q(3*n, s_1n) + s * (gsl_sf_gamma_inc(2*n, s_1n)) / gsl_sf_gamma(3*n)
+          )
+      );
+    }
+    double EinastoPotentialrforce(double r,double t,
+      struct potentialArg * potentialArgs){
+      double * args= potentialArgs->args;
+      //Get args
+      double amp= *args;
+      double h= *(args+1);
+      double n= *(args+2);
+      double s = r / h;
+      double gamma_3n = gsl_sf_gamma(3 * n);
+      double gamma_upper_3n = gsl_sf_gamma_inc_Q(3 * n, pow(s, 1.0 / n));
+      double s_2 = pow(s, -2.0);
+      return (
+        (4.0 * M_PI * h * n * gamma_3n)
+        * s_2
+        * (gamma_upper_3n - 1.0)
+      );
+    }
+    double EinastoPotentialr2deriv(double r,double t,
+       struct potentialArg * potentialArgs){
+      double * args= potentialArgs->args;
+      //Get args
+      double amp= *args;
+      double h= *(args+1);
+      double n= *(args+2);
+      double s = r / h;
+      double s_1n = pow(s,1.0/n);
+      double gamma_3n = gsl_sf_gamma(3*n);
+      double gamma_upper_3n = gsl_sf_gamma_inc_Q(3*n, s_1n);
+      return (
+          - (4.0 * M_PI * n * gamma_3n)
+          * ((-2 * pow(s,-3)) * (gamma_upper_3n - 1)
+          - ((1/n) * exp(-s_1n)/gamma_3n))
+      );
+    }
+    double EinastoPotentialrdens(double r,
+    			    double t,
+    			    struct potentialArg * potentialArgs){
+      double * args= potentialArgs->args;
+      //Get args
+      double amp= *args;
+      double h= *(args+1);
+      double n= *(args+2);
+      //Calculate potential
+      return exp(-pow(r/h, 1.0 / n));
+    }

galpy/potential/potential_c_ext/galpy_potentials.h

-Original file line number
+Diff line change
@@ Expand Up @@
     double interpSphericalPotentialrforce(double,double,struct potentialArg *);
     double interpSphericalPotentialr2deriv(double,double,struct potentialArg *);
     double interpSphericalPotentialrdens(double,double,struct potentialArg *);
+    //EinastoPotential: uses SphericalPotential, only need revaluate, rforce, r2deriv
+    double EinastoPotentialrevaluate(double,double,struct potentialArg *);
+    double EinastoPotentialrforce(double,double,struct potentialArg *);
+    double EinastoPotentialr2deriv(double,double,struct potentialArg *);
+    double EinastoPotentialrdens(double,double,struct potentialArg *);
     //TriaxialGaussian: uses EllipsoidalPotential, only need psi, dens, densDeriv
     double TriaxialGaussianPotentialpsi(double,double *);
@@ Expand Down @@

tests/test_actionAngle.py

-Original file line number
+Diff line change
@@ Expand Up @@
             rforce=potential.HomogeneousSpherePotential(normalize=1.0, R=1.1),
             rgrid=numpy.linspace(0.0, 1.1, 201),
         )
+        ep = potential.EinastoPotential(normalize=1.0, h=2.2)
         msmlpwtdp = (
             mockSmoothedLogarithmicHaloPotentialwTimeDependentAmplitudeWrapperPotential()
         )
@@ Expand All @@
             homp,
             ihomp,
             msmlpwtdp,
+            ep,
         ]
         for pot in pots:
             aAS = actionAngleStaeckel(pot=pot, c=True, delta=0.01)
@@ Expand Down @@

tests/test_dynamfric.py

-Original file line number
+Diff line change
@@ Expand Up / @@ -281,6 +281,7 @@ def test_dynamfric_c(): @@
                 normalize=1.0,
                 a=3.5,
             ),
+            potential.EinastoPotential(normalize=1.0, h=2.2),
             MWPotential3021,
             McMillan17,  # SCF + DiskSCF
         ]
@@ Expand Down @@

tests/test_orbit.py

-Original file line number
+Diff line change
@@ Expand Up / @@ -1893,6 +1893,7 @@ def test_analytic_ecc_rperi_rap(): @@
         tol["default"] = -10.0
         tol["NFWPotential"] = -9.0  # these are more difficult
         tol["PlummerPotential"] = -9.0  # these are more difficult
+        tol["EinastoPotential"] = -9.0  # these are more difficult
         tol["DoubleExponentialDiskPotential"] = -6.0  # these are more difficult
         tol["RazorThinExponentialDiskPotential"] = -8.0  # these are more difficult
         tol["IsochronePotential"] = -6.0  # these are more difficult
@@ Expand Down @@

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add Einasto potential #760

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!

Uh oh!