Added Burkert cpotential

.gitignore

@@ -76,4 +76,6 @@ pip-wheel-metadata/
 notebooks
 plots
 scripts
-.hypothesis
+.hypothesis
+
+gala_env

CHANGES.rst

@@ -7,6 +7,8 @@ New Features
 - Added an option to specify a multiprocessing or parallel processing pool when
   computing basis function coefficients for the SCF potential from discrete particles.
 
+- Added the Burkert potential as a built-in cpotential.
+
 Bug fixes
 ---------
 

gala/potential/potential/builtin/builtin_potentials.c

@@ -1980,3 +1980,56 @@ void longmuralibar_hessian(double t, double *pars, double *q, int n_dim,
     hess[7] = hess[7] + tmp_88;
     hess[8] = hess[8] + tmp_38*tmp_76*tmp_94 - tmp_40*tmp_77*tmp_94 + tmp_52*(tmp_14*tmp_92*tmp_96 - tmp_22*tmp_45*tmp_97 - tmp_28*tmp_78*tmp_98 + tmp_35*tmp_48*tmp_97 + tmp_89*tmp_95 - tmp_90*tmp_96 + tmp_91 - tmp_92*tmp_95 - tmp_93 + tmp_50*tmp_98/tmp_17);
 }
+
+
+/* ---------------------------------------------------------------------------
+    Burkert potential
+    (from Mori and Burkert 2000: https://iopscience.iop.org/article/10.1086/309140/fulltext/50172.text.html)
+*/
+double burkert_value(double t, double *pars, double *q, int n_dim) {
+    /*  pars:
+            - G (Gravitational constant)
+            - 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];
+
+    // 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) );
+}
+
+
+void burkert_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
+    /*  pars:
+            - G (Gravitational constant)
+            - 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];
+
+    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;
+}
+
+
+double burkert_density(double t, double *pars, double *q, int n_dim) {
+    /*  pars:
+            - G (Gravitational constant)
+            - 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;
+}

gala/potential/potential/builtin/builtin_potentials.h

@@ -92,3 +92,7 @@ extern double longmuralibar_value(double t, double *pars, double *q, int n_dim);
 extern void longmuralibar_gradient(double t, double *pars, double *q, int n_dim, double *grad);
 extern double longmuralibar_density(double t, double *pars, double *q, int n_dim);
 extern void longmuralibar_hessian(double t, double *pars, double *q, int n_dim, double *hess);
+
+extern double burkert_value(double t, double *pars, double *q, int n_dim);
+extern void burkert_gradient(double t, double *pars, double *q, int n_dim, double *grad);
+extern double burkert_density(double t, double *pars, double *q, int n_dim);

gala/potential/potential/builtin/core.py

@@ -38,7 +38,8 @@
     LongMuraliBarWrapper,
     NullWrapper,
     MultipoleWrapper,
-    CylSplineWrapper
+    CylSplineWrapper,
+    BurkertWrapper
 )
 
 __all__ = [
@@ -60,8 +61,9 @@
     "LogarithmicPotential",
     "LongMuraliBarPotential",
     "MultipolePotential",
-    "CylSplinePotential"
-]
+    "CylSplinePotential",
+    "BurkertPotential"
+]   
 
 
 def __getattr__(name):
@@ -378,6 +380,26 @@ def to_sympy(cls, v, p):
         return expr, v, p
 
 
+@format_doc(common_doc=_potential_docstring)
+class BurkertPotential(CPotentialBase):
+    r"""
+    The Burkert potential that well-matches the rotation curve of dwarf galaxies.
+    See https://iopscience.iop.org/article/10.1086/309140/fulltext/50172.text.html
+
+    Parameters
+    ----------
+    rho : :class:`~astropy.units.Quantity`, numeric [mass density]
+        Central mass density.
+    r0 : :class:`~astropy.units.Quantity`, numeric [length]
+        The core radius.
+    {common_doc}
+    """
+    rho = PotentialParameter("rho", physical_type="mass density")
+    r0 = PotentialParameter("r0", physical_type="length")
+
+    Wrapper = BurkertWrapper
+
+
 # ============================================================================
 # Flattened, axisymmetric models
 #

gala/potential/potential/builtin/cybuiltin.pyx

@@ -121,6 +121,11 @@ cdef extern from "potential/potential/builtin/builtin_potentials.h":
     double longmuralibar_density(double t, double *pars, double *q, int n_dim) nogil
     void longmuralibar_hessian(double t, double *pars, double *q, int n_dim, double *hess) nogil
 
+    double burkert_value(double t, double *pars, double *q, int n_dim) nogil
+    void burkert_gradient(double t, double *pars, double *q, int n_dim, double *grad) nogil
+    double burkert_density(double t, double *pars, double *q, int n_dim) nogil
+
+
 cdef extern from "potential/potential/builtin/multipole.h":
     double mp_potential(double t, double *pars, double *q, int n_dim) nogil
     void mp_gradient(double t, double *pars, double *q, int n_dim, double *grad) nogil
@@ -175,6 +180,7 @@ __all__ = [
     'NullWrapper',
     'MultipoleWrapper',
     'CylSplineWrapper'
+    'BurkertWrapper'
 ]
 
 # ============================================================================
@@ -279,6 +285,16 @@ 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):
+        self.init([G] + list(parameters),
+                  np.ascontiguousarray(q0),
+                  np.ascontiguousarray(R))
+        self.cpotential.value[0] = <energyfunc>(burkert_value)
+        self.cpotential.density[0] = <densityfunc>(burkert_density)
+        self.cpotential.gradient[0] = <gradientfunc>(burkert_gradient)
+
 
 # ============================================================================
 # Flattened, axisymmetric models

gala/potential/potential/tests/test_all_builtin.py

@@ -538,3 +538,18 @@ def test_against_agama(self):
         assert u.allclose(gala_ene, agama_tbl['pot'][:, 0], rtol=1e-3)
         for i in range(3):
             assert u.allclose(gala_acc[i], agama_tbl['acc'][:, i], rtol=1e-2)
+
+
+class TestBurkert(PotentialTestBase):
+    potential = p.BurkertPotential(units=galactic, rho=5e-25 *u.g/u.cm**3, r0=12 * u.kpc)
+    w0 = [1., 0., 0., 0., 0.1, 0.1]
+
+    check_finite_at_origin = False
+
+    @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