Repository: adrn/gala (language: python)

Issue:
Implement more compatibility with galpy (to/from methods)
#185 started this by allowing exporting Gala orbit objects to Galpy orbits, but other things we can do:

- [x] Add a `from_galpy_orbit()` classmethod
- [x] Convert potentials to galpy potentials


Discussion / hints:

Fix patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index b67f12fc..5463edef 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -7,6 +7,10 @@ New Features
 - ``UnitSystem`` objects can now be created with custom units passed in as
   Astropy ``Quantity`` objects.
 
+- Added functionality to convert Gala potential objects to Galpy potential
+  objects, or to create Gala potential objects from a pre-existing Galpy
+  potential.
+
 
 Bug fixes
 ---------
@@ -17,6 +21,12 @@ Bug fixes
 - Fixed a bug with ``GreatCircleICRSFrame.from_endpoints()`` that caused an
   error when the input coordinates had associated velocity data.
 
+- Fixed a bug with the ``JaffePotential`` density evaluation, which was too low
+  by a factor of two.
+
+- Implemented a density function for ``LogarithmicPotential``, which was
+  missing previously.
+
 
 API changes
 -----------
diff --git a/docs/interop.rst b/docs/interop.rst
new file mode 100644
index 00000000..c319a714
--- /dev/null
+++ b/docs/interop.rst
@@ -0,0 +1,96 @@
+.. _gala-interop:
+
+*********************************************
+Interoperability with Other Dynamics Packages
+*********************************************
+
+Gala provides interfaces with other common Galactic dynamics packages, which
+enables easily converting objects between these packages. Some examples are
+shown below. As always, if something does not work as expected or you would like
+more interoperability with any of these packages, please `open an issue
+<https://github.com/adrn/gala/issues/new>`_ on GitHub.
+
+Here are some imports we will use below in examples::
+
+    >>> import astropy.units as u
+    >>> import numpy as np
+    >>> import gala.dynamics as gd
+    >>> import gala.potential as gp
+    >>> from gala.units import galactic
+
+Galpy
+=====
+
+`Galpy <https://docs.galpy.org/en/>`_ is another popular Python package for
+Galactic dynamics with similar functionality to Gala. For example, Galpy
+supports creating gravitational potential objects and numerically integrating
+orbits (among other things).
+
+Gala provides an interface for converting representations of orbits from Gala to
+Galpy, or from Galpy to Gala. To convert a Gala :class:`~gala.dynamics.Orbit`
+object to a Galpy ``Orbit``, use the
+:meth:`~gala.dynamics.Orbit.to_galpy_orbit()` method:
+
+.. doctest-requires:: galpy
+
+    >>> w0 = gd.PhaseSpacePosition(pos=[10., 0, 0] * u.kpc,
+    ...                            vel=[0, 0, 200.] * u.km/u.s)
+    >>> mw = gp.Hamiltonian(gp.MilkyWayPotential())
+    >>> orbit = mw.integrate_orbit(w0, dt=1, n_steps=1000)
+    >>> orbit
+    <Orbit cartesian, dim=3, shape=(1001,)>
+    >>> galpy_orbit = orbit.to_galpy_orbit()
+    >>> galpy_orbit  # doctest: +SKIP
+    <galpy.orbit.Orbits.Orbit object at 0x7f99b3c8fe80>
+
+Similarly, a Galpy ``Orbit`` can be converted to a Gala
+:class:`~gala.dynamics.Orbit` using the
+:meth:`~gala.dynamics.Orbit.from_galpy_orbit()` classmethod:
+
+.. doctest-requires:: galpy
+
+    >>> import galpy.potential as galpy_p
+    >>> import galpy.orbit as galpy_o
+    >>> mp = galpy_p.MiyamotoNagaiPotential(a=0.5, b=0.0375, amp=1.,
+    ...                                     normalize=1.)
+    >>> galpy_orbit = galpy_o.Orbit([1., 0.1, 1.1, 0., 0.1, 1.])
+    >>> ts = np.linspace(0, 100, 10000)
+    >>> galpy_orbit.integrate(ts, mp, method='odeint')
+    >>> orbit = gd.Orbit.from_galpy_orbit(galpy_orbit)
+    >>> orbit
+    <Orbit cylindrical, dim=3, shape=(10000,)>
+
+Gala also provides tools for converting potential objects to `galpy` potential
+objects, or creating Gala potential objects from existing `galpy` potentials.
+To convert a Gala potential to a Galpy potential, use the
+:meth:`~gala.potential.PotentialBase.to_galpy_potential()` method on any Gala
+potential object. For example:
+
+.. doctest-requires:: galpy
+
+    >>> pot = gp.HernquistPotential(m=1e10*u.Msun, c=1.5*u.kpc, units=galactic)
+    >>> galpy_pot = pot.to_galpy_potential()
+    >>> galpy_pot  # doctest: +SKIP
+    <galpy.potential.TwoPowerSphericalPotential.HernquistPotential at 0x7faa00263b20>
+    >>> galpy_pot.Rforce(1., 0.)  # doctest: +FLOAT_CMP
+    -0.48737954713808573
+
+To convert from a Galpy potential to a Gala potential, use the
+:func:`~gala.potential.potential.galpy_to_gala_potential()` function. For
+example:
+
+.. doctest-requires:: galpy
+
+    >>> import galpy.potential as galpy_gp
+    >>> galpy_pot = galpy_gp.HernquistPotential(amp=1., a=0.5)
+    >>> pot = gp.galpy_to_gala_potential(galpy_pot)
+    >>> pot
+    <HernquistPotential: m=4.50e+10, c=4.00 (kpc,Myr,solMass,rad)>
+
+
+Agama
+=====
+
+Coming soon, but we could use your help! Please leave a note `in this issue
+<https://github.com/adrn/gala/issues/230>`_ if you would find interoperability
+with Agama useful.
diff --git a/docs/user_guide.rst b/docs/user_guide.rst
index 6bb017f8..7402dcf5 100644
--- a/docs/user_guide.rst
+++ b/docs/user_guide.rst
@@ -22,6 +22,7 @@ pieces of ``gala`` work together to meet more realistic research needs.
     dynamics/index
     units
     util
+    interop
 
 .. toctree::
     :hidden:
diff --git a/gala/dynamics/core.py b/gala/dynamics/core.py
index 72ba39f2..8570f343 100644
--- a/gala/dynamics/core.py
+++ b/gala/dynamics/core.py
@@ -44,7 +44,7 @@ class RegexRepresentationMapping(RepresentationMapping):
     pass
 
 
-class PhaseSpacePosition(object):
+class PhaseSpacePosition:
 
     representation_mappings = {
         r.CartesianRepresentation: [
@@ -857,57 +857,3 @@ def shape(self):
         accessed by doing, e.g., ``obj.x.shape``.
         """
         return self.pos.shape
-
-    #
-    # Compatibility with other packages
-    #
-
-    def to_galpy_orbit(self, ro=None, vo=None):
-        """Convert this object to a ``galpy.Orbit`` instance.
-
-        Parameters
-        ----------
-        ro : `astropy.units.Quantity` or `astropy.units.UnitBase`
-            "Natural" length unit.
-        vo : `astropy.units.Quantity` or `astropy.units.UnitBase`
-            "Natural" velocity unit.
-
-        Returns
-        -------
-        galpy_orbit : `galpy.orbit.Orbit`
-
-        """
-        from galpy.orbit import Orbit
-        from galpy.util.config import __config__ as galpy_config
-
-        if self.frame is not None:
-            from ..potential import StaticFrame
-            w = self.to_frame(StaticFrame(self.frame.units))
-        else:
-            w = self
-
-        if ro is None:
-            ro = galpy_config.getfloat('normalization', 'ro')
-            ro = ro * u.kpc
-
-        if vo is None:
-            vo = galpy_config.getfloat('normalization', 'vo')
-            vo = vo * u.km/u.s
-
-        # PhaseSpacePosition or Orbit:
-        cyl = w.cylindrical
-
-        R = cyl.rho.to_value(ro).T
-        phi = cyl.phi.to_value(u.rad).T
-        z = cyl.z.to_value(ro).T
-
-        vR = cyl.v_rho.to_value(vo).T
-        vT = (cyl.rho * cyl.pm_phi).to_value(vo, u.dimensionless_angles()).T
-        vz = cyl.v_z.to_value(vo).T
-
-        o = Orbit(np.array([R, vR, vT, z, vz, phi]).T, ro=ro, vo=vo)
-
-        if hasattr(w, 't') and w.t is not None:
-            o.t = w.t.to_value(ro / vo)
-
-        return o
diff --git a/gala/dynamics/orbit.py b/gala/dynamics/orbit.py
index b1a471cd..1c25091d 100644
--- a/gala/dynamics/orbit.py
+++ b/gala/dynamics/orbit.py
@@ -3,6 +3,7 @@
 
 # Third-party
 from astropy import log as logger
+import astropy.coordinates as coord
 import astropy.units as u
 import numpy as np
 from scipy.signal import argrelmax
@@ -932,3 +933,89 @@ def to_frame(self, frame, current_frame=None, **kwargs):
 
         return Orbit(pos=psp.pos, vel=psp.vel, t=self.t,
                      frame=frame, potential=self.potential)
+
+    # ------------------------------------------------------------------------
+    # Compatibility with other packages
+    #
+
+    def to_galpy_orbit(self, ro=None, vo=None):
+        """Convert this object to a ``galpy.Orbit`` instance.
+
+        Parameters
+        ----------
+        ro : `astropy.units.Quantity` or `astropy.units.UnitBase`
+            "Natural" length unit.
+        vo : `astropy.units.Quantity` or `astropy.units.UnitBase`
+            "Natural" velocity unit.
+
+        Returns
+        -------
+        galpy_orbit : `galpy.orbit.Orbit`
+
+        """
+        from galpy.orbit import Orbit
+        from galpy.util.config import __config__ as galpy_config
+
+        if self.frame is not None:
+            from ..potential import StaticFrame
+            w = self.to_frame(StaticFrame(self.frame.units))
+        else:
+            w = self
+
+        if ro is None:
+            ro = galpy_config.getfloat('normalization', 'ro')
+            ro = ro * u.kpc
+
+        if vo is None:
+            vo = galpy_config.getfloat('normalization', 'vo')
+            vo = vo * u.km/u.s
+
+        # PhaseSpacePosition or Orbit:
+        cyl = w.cylindrical
+
+        R = cyl.rho.to_value(ro).T
+        phi = cyl.phi.to_value(u.rad).T
+        z = cyl.z.to_value(ro).T
+
+        vR = cyl.v_rho.to_value(vo).T
+        vT = (cyl.rho * cyl.pm_phi).to_value(vo, u.dimensionless_angles()).T
+        vz = cyl.v_z.to_value(vo).T
+
+        o = Orbit(np.array([R, vR, vT, z, vz, phi]).T, ro=ro, vo=vo)
+        if w.t is not None:
+            o.t = w.t.to_value(ro / vo)
+
+        return o
+
+    @classmethod
+    def from_galpy_orbit(self, galpy_orbit):
+        """Create a Gala ``PhaseSpacePosition`` or ``Orbit`` instance from a
+        ``galpy.Orbit`` instance.
+
+        Parameters
+        ----------
+        galpy_orbit : :class:`galpy.orbit.Orbit`
+
+        Returns
+        -------
+        orbit : :class:`~gala.dynamics.Orbit`
+
+        """
+        ro = galpy_orbit._ro * u.kpc
+        vo = galpy_orbit._vo * u.km/u.s
+        ts = galpy_orbit.t
+
+        rep = coord.CylindricalRepresentation(
+            rho=galpy_orbit.R(ts) * ro,
+            phi=galpy_orbit.phi(ts) * u.rad,
+            z=galpy_orbit.z(ts) * ro
+        )
+        with u.set_enabled_equivalencies(u.dimensionless_angles()):
+            dif = coord.CylindricalDifferential(
+                d_rho=galpy_orbit.vR(ts) * vo,
+                d_phi=galpy_orbit.vT(ts) * vo / rep.rho,
+                d_z=galpy_orbit.vz(ts) * vo
+            )
+
+        t = galpy_orbit.t * ro / vo
+        return Orbit(rep, dif, t=t)
diff --git a/gala/potential/common.py b/gala/potential/common.py
index 91c3fb3e..63af1fe7 100644
--- a/gala/potential/common.py
+++ b/gala/potential/common.py
@@ -39,6 +39,9 @@ def __init__(self, name, physical_type="dimensionless", default=None,
         self.default = default
         self.equivalencies = equivalencies
 
+    def __repr__(self):
+        return f"<PotentialParameter: {self.name} [{self.physical_type}]>"
+
 
 class CommonBase:
 
diff --git a/gala/potential/potential/__init__.py b/gala/potential/potential/__init__.py
index f727f6cf..7e8dfe83 100644
--- a/gala/potential/potential/__init__.py
+++ b/gala/potential/potential/__init__.py
@@ -4,3 +4,4 @@
 from .builtin import *
 from .io import *
 from .util import *
+from .interop import *
diff --git a/gala/potential/potential/builtin/builtin_potentials.c b/gala/potential/potential/builtin/builtin_potentials.c
index fb0b9265..3f820ef9 100644
--- a/gala/potential/potential/builtin/builtin_potentials.c
+++ b/gala/potential/potential/builtin/builtin_potentials.c
@@ -1,5 +1,6 @@
 #include <math.h>
 #include <string.h>
+#include <stdio.h>
 #include "extra_compile_macros.h"
 
 #if USE_GSL == 1
@@ -406,7 +407,7 @@ double jaffe_density(double t, double *pars, double *q, int n_dim) {
     */
     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]);
+    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));
 }
 
@@ -540,7 +541,7 @@ double powerlawcutoff_density(double t, double *pars, double *q, int n_dim) {
     */
     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]));
+    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]));
 }
 
@@ -1222,6 +1223,59 @@ void satoh_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
     hess[8] = hess[8] + -tmp_13*(-tmp_12*z - z) - tmp_7*(a*tmp_2/pow(tmp_3, 3.0/2.0) - tmp_12 - 1);
 }
 
+/* ---------------------------------------------------------------------------
+    Kuzmin potential
+*/
+double kuzmin_value(double t, double *pars, double *q, int n_dim) {
+    /*  pars:
+            - G (Gravitational constant)
+            - m (mass scale)
+            - a (length scale 1) TODO
+    */
+    double S2 = q[0]*q[0] + q[1]*q[1] + pow(pars[2] + fabs(q[2]), 2);
+    return -pars[0] * pars[1] / sqrt(S2);
+}
+
+void kuzmin_gradient(double t, double *pars, double *q, int n_dim,
+                     double *grad) {
+    /*  pars:
+            - G (Gravitational constant)
+            - m (mass scale)
+            - a (length scale 1) TODO
+    */
+
+    double S2 = q[0]*q[0] + q[1]*q[1] + pow(pars[2] + fabs(q[2]), 2);
+    double fac = pars[0] * pars[1] * pow(S2, -1.5);
+    double zsign;
+
+    if (q[2] > 0) {
+        zsign = 1.;
+    } else if (q[2] < 0) {
+        zsign = -1.;
+    } else {
+        zsign = 0.;
+    }
+
+    grad[0] = grad[0] + fac * q[0];
+    grad[1] = grad[1] + fac * q[1];
+    grad[2] = grad[2] + fac * zsign * (pars[2] + fabs(q[2]));
+}
+
+double kuzmin_density(double t, double *pars, double *q, int n_dim) {
+    /*  pars:
+            - G (Gravitational constant)
+            - m (mass scale)
+            - a (length scale 1) TODO
+    */
+    if (q[2] != 0.) {
+        return 0.;
+    } else {
+        return pars[1] * pars[2] / (2 * M_PI) *
+            pow(q[0]*q[0] + q[1]*q[1] + pars[2]*pars[2], -1.5);
+   }
+
+}
+
 /* ---------------------------------------------------------------------------
     Miyamoto-Nagai flattened potential
 */
@@ -1457,6 +1511,30 @@ double logarithmic_value(double t, double *pars, double *q, int n_dim) {
                                      z*z/(pars[5]*pars[5]));
 }
 
+double logarithmic_density(double t, double *pars, double *q, int n_dim) {
+    /*  pars:
+            - G (Gravitational constant)
+            - v_c (velocity scale)
+            - r_h (length scale)
+            - q1
+            - q2
+            - q3
+    */
+    double tmp_0 = pow(pars[3], 2);
+    double tmp_1 = pow(pars[4], 2);
+    double tmp_2 = tmp_0*tmp_1;
+    double tmp_3 = tmp_2*pow(q[2], 2);
+    double tmp_4 = pow(pars[5], 2);
+    double tmp_5 = tmp_0*tmp_4;
+    double tmp_6 = tmp_5*pow(q[1], 2);
+    double tmp_7 = tmp_1*tmp_4;
+    double tmp_8 = tmp_7*pow(q[0], 2);
+    double tmp_9 = pow(pars[2], 2)*tmp_2*tmp_4;
+    double tmp_10 = tmp_6 + tmp_8 + tmp_9;
+    double tmp_11 = tmp_3 + tmp_9;
+    return pow(pars[1], 2)*(tmp_2*(tmp_10 - tmp_3) + tmp_5*(tmp_11 - tmp_6 + tmp_8) + tmp_7*(tmp_11 + tmp_6 - tmp_8))/pow(tmp_10 + tmp_3, 2) / (4*M_PI*pars[0]);
+}
+
 void logarithmic_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
     /*  pars:
             - G (Gravitational constant)
diff --git a/gala/potential/potential/builtin/builtin_potentials.h b/gala/potential/potential/builtin/builtin_potentials.h
index e5624db9..4650bd95 100644
--- a/gala/potential/potential/builtin/builtin_potentials.h
+++ b/gala/potential/potential/builtin/builtin_potentials.h
@@ -65,6 +65,10 @@ extern void satoh_gradient(double t, double *pars, double *q, int n_dim, double
 extern double satoh_density(double t, double *pars, double *q, int n_dim);
 extern void satoh_hessian(double t, double *pars, double *q, int n_dim, double *hess);
 
+extern double kuzmin_value(double t, double *pars, double *q, int n_dim);
+extern void kuzmin_gradient(double t, double *pars, double *q, int n_dim, double *grad);
+extern double kuzmin_density(double t, double *pars, double *q, int n_dim);
+
 extern double miyamotonagai_value(double t, double *pars, double *q, int n_dim);
 extern void miyamotonagai_gradient(double t, double *pars, double *q, int n_dim, double *grad);
 extern void miyamotonagai_hessian(double t, double *pars, double *q, int n_dim, double *hess);
@@ -77,6 +81,7 @@ extern double leesuto_density(double t, double *pars, double *q, int n_dim);
 extern double logarithmic_value(double t, double *pars, double *q, int n_dim);
 extern void logarithmic_gradient(double t, double *pars, double *q, int n_dim, double *grad);
 extern void logarithmic_hessian(double t, double *pars, double *q, int n_dim, double *hess);
+extern double logarithmic_density(double t, double *pars, double *q, int n_dim);
 
 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);
diff --git a/gala/potential/potential/builtin/cybuiltin.pyx b/gala/potential/potential/builtin/cybuiltin.pyx
index 035ed906..e97f971c 100644
--- a/gala/potential/potential/builtin/cybuiltin.pyx
+++ b/gala/potential/potential/builtin/cybuiltin.pyx
@@ -100,6 +100,10 @@ cdef extern from "potential/potential/builtin/builtin_potentials.h":
     double satoh_density(double t, double *pars, double *q, int n_dim) nogil
     void satoh_hessian(double t, double *pars, double *q, int n_dim, double *hess) nogil
 
+    double kuzmin_value(double t, double *pars, double *q, int n_dim) nogil
+    void kuzmin_gradient(double t, double *pars, double *q, int n_dim, double *grad) nogil
+    double kuzmin_density(double t, double *pars, double *q, int n_dim) nogil
+
     double miyamotonagai_value(double t, double *pars, double *q, int n_dim) nogil
     void miyamotonagai_gradient(double t, double *pars, double *q, int n_dim, double *grad) nogil
     void miyamotonagai_hessian(double t, double *pars, double *q, int n_dim, double *hess) nogil
@@ -112,6 +116,7 @@ cdef extern from "potential/potential/builtin/builtin_potentials.h":
     double logarithmic_value(double t, double *pars, double *q, int n_dim) nogil
     void logarithmic_gradient(double t, double *pars, double *q, int n_dim, double *grad) nogil
     void logarithmic_hessian(double t, double *pars, double *q, int n_dim, double *hess) nogil
+    double logarithmic_density(double t, double *pars, double *q, int n_dim) nogil
 
     double longmuralibar_value(double t, double *pars, double *q, int n_dim) nogil
     void longmuralibar_gradient(double t, double *pars, double *q, int n_dim, double *grad) nogil
@@ -121,7 +126,7 @@ cdef extern from "potential/potential/builtin/builtin_potentials.h":
 __all__ = ['NullPotential', 'HenonHeilesPotential', # Misc. potentials
            'KeplerPotential', 'HernquistPotential', 'IsochronePotential', 'PlummerPotential',
            'JaffePotential', 'StonePotential', 'PowerLawCutoffPotential', # Spherical models
-           'SatohPotential', 'MiyamotoNagaiPotential', # Disk models
+           'SatohPotential', 'KuzminPotential', 'MiyamotoNagaiPotential', # Disk models
            'NFWPotential', 'LeeSutoTriaxialNFWPotential', 'LogarithmicPotential',
            'LongMuraliBarPotential', # Triaxial models
            ]
@@ -544,6 +549,46 @@ class SatohPotential(CPotentialBase):
         return expr, v, p
 
 
+cdef class KuzminWrapper(CPotentialWrapper):
+
+    def __init__(self, G, parameters, q0, R):
+        self.init([G] + list(parameters),
+                  np.ascontiguousarray(q0),
+                  np.ascontiguousarray(R))
+        self.cpotential.value[0] = <energyfunc>(kuzmin_value)
+        self.cpotential.density[0] = <densityfunc>(kuzmin_density)
+        self.cpotential.gradient[0] = <gradientfunc>(kuzmin_gradient)
+        self.cpotential.hessian[0] = <hessianfunc>(null_hessian)
+
+@format_doc(common_doc=_potential_docstring)
+class KuzminPotential(CPotentialBase):
+    r"""
+    KuzminPotential(m, a, units=None, origin=None, R=None)
+
+    Kuzmin potential for a flattened mass distribution.
+
+    Parameters
+    ----------
+    m : :class:`~astropy.units.Quantity`, numeric [mass]
+        Mass.
+    a : :class:`~astropy.units.Quantity`, numeric [length]
+        Flattening parameter.
+    {common_doc}
+    """
+    m = PotentialParameter('m', physical_type='mass')
+    a = PotentialParameter('a', physical_type='length')
+
+    Wrapper = KuzminWrapper
+
+    @myclassmethod
+    @sympy_wrap
+    def to_sympy(cls, v, p):
+        import sympy as sy
+        denom = sy.sqrt(v['x']**2 + v['y']**2 + (p['a'] + sy.Abs(v['z']))**2)
+        expr = - p['G'] * p['m'] / denom
+        return expr, v, p
+
+
 cdef class MiyamotoNagaiWrapper(CPotentialWrapper):
 
     def __init__(self, G, parameters, q0, R):
@@ -787,6 +832,7 @@ cdef class LogarithmicWrapper(CPotentialWrapper):
         self.cpotential.value[0] = <energyfunc>(logarithmic_value)
         self.cpotential.gradient[0] = <gradientfunc>(logarithmic_gradient)
         self.cpotential.hessian[0] = <hessianfunc>(logarithmic_hessian)
+        self.cpotential.density[0] = <energyfunc>(logarithmic_density)
 
 @format_doc(common_doc=_potential_docstring)
 class LogarithmicPotential(CPotentialBase):
@@ -927,9 +973,9 @@ class LongMuraliBarPotential(CPotentialBase):
         y = -v['x'] * sy.sin(p['alpha']) + v['y'] * sy.cos(p['alpha'])
         z = v['z']
 
-        Tm = sy.sqrt((p['a']-x)**2 + y**2 +
+        Tm = sy.sqrt((p['a'] - x)**2 + y**2 +
                      (p['b'] + sy.sqrt(p['c']**2 + z**2))**2)
-        Tp = sy.sqrt((p['a']+x)**2 + y**2 +
+        Tp = sy.sqrt((p['a'] + x)**2 + y**2 +
                      (p['b'] + sy.sqrt(p['c']**2 + z**2))**2)
 
         expr = (p['G'] * p['m'] / (2*p['a']) *
diff --git a/gala/potential/potential/builtin/pybuiltin.py b/gala/potential/potential/builtin/pybuiltin.py
index fbda4884..dfede461 100644
--- a/gala/potential/potential/builtin/pybuiltin.py
+++ b/gala/potential/potential/builtin/pybuiltin.py
@@ -5,7 +5,7 @@
 from gala.potential.potential.util import sympy_wrap
 from gala.potential.common import PotentialParameter
 
-__all__ = ["HarmonicOscillatorPotential", "KuzminPotential"]
+__all__ = ["HarmonicOscillatorPotential"]
 
 
 class HarmonicOscillatorPotential(PotentialBase):
@@ -88,48 +88,3 @@ def action_angle(self, w):
     #     """
     #     from ...dynamics.analyticactionangle import harmonic_oscillator_aa_to_xv
     #     return harmonic_oscillator_aa_to_xv(actions, angles, self)
-
-
-class KuzminPotential(PotentialBase):
-    r"""
-    The Kuzmin flattened disk potential.
-
-    .. math::
-
-        \Phi = -\frac{Gm}{\sqrt{x^2 + y^2 + (a + |z|)^2}}
-
-    Parameters
-    ----------
-    m : numeric
-        Mass.
-    a : numeric
-        Flattening parameter.
-    units : iterable
-        Unique list of non-reducible units that specify (at minimum) the
-        length, mass, time, and angle units.
-
-    """
-    m = PotentialParameter('m', physical_type='mass')
-    a = PotentialParameter('a', physical_type='length')
-
-    def _energy(self, q, t):
-        x, y, z = q
-        m = self.parameters['m']
-        a = self.parameters['a']
-        val = -self.G * m / np.sqrt(x**2 + y**2 + (a + np.abs(z))**2)
-        return val
-
-    def _gradient(self, q, t):
-        x, y, z = q
-        m = self.parameters['m']
-        a = self.parameters['a']
-        fac = self.G * m / (x**2 + y**2 + (a + np.abs(z))**2)**1.5
-        return fac[None, ...] * q
-
-    @classmethod
-    @sympy_wrap()
-    def to_sympy(cls, v, p):
-        import sympy as sy
-        denom = sy.sqrt(v['x']**2 + v['y']**2 + (p['a'] + sy.Abs(v['z']))**2)
-        expr = p['G'] * p['m'] / denom
-        return expr, v, p
diff --git a/gala/potential/potential/core.py b/gala/potential/potential/core.py
index fe62edbc..a577ee12 100644
--- a/gala/potential/potential/core.py
+++ b/gala/potential/potential/core.py
@@ -12,9 +12,11 @@
 from astropy.utils import isiterable
 try:
     from scipy.spatial.transform import Rotation
-except ImportError:
-    raise ImportError("Gala requires scipy>=1.2: make sure you have updated "
-                      "your version of scipy and try importing gala again.")
+except ImportError as exc:
+    raise ImportError(
+        "Gala requires scipy>=1.2: make sure you have updated your version of "
+        "scipy and try importing gala again."
+    ) from exc
 
 # Project
 from ..common import CommonBase
@@ -29,12 +31,11 @@ class PotentialBase(CommonBase, metaclass=abc.ABCMeta):
     """
     A baseclass for defining pure-Python gravitational potentials.
 
-    Subclasses must define (at minimum) a method that evaluates
-    the potential energy at a given position ``q``
-    and time ``t``: ``_energy(q, t)``. For integration, the subclasses
-    must also define a method to evaluate the gradient,
-    ``_gradient(q, t)``. Optionally, they may also define methods
-    to compute the density and hessian: ``_density()``, ``_hessian()``.
+    Subclasses must define (at minimum) a method that evaluates the potential
+    energy at a given position ``q`` and time ``t``: ``_energy(q, t)``. For
+    integration, the subclasses must also define a method to evaluate the
+    gradient, ``_gradient(q, t)``. Optionally, they may also define methods to
+    compute the density and hessian: ``_density()``, ``_hessian()``.
     """
     ndim = 3
 
@@ -120,7 +121,7 @@ def to_latex(cls):
         import sympy as sy
         return sy.latex(expr)
 
-    # ========================================================================
+    ###########################################################################
     # Abstract methods that must be implemented by subclasses
     #
     @abc.abstractmethod
@@ -138,7 +139,7 @@ def _density(self, q, t=0.):
     def _hessian(self, q, t=0.):
         raise NotImplementedError("This Potential has no implemented Hessian.")
 
-    # ========================================================================
+    ###########################################################################
     # Utility methods
     #
     def _remove_units(self, x):
@@ -175,7 +176,7 @@ def _remove_units_prepare_shape(self, x):
 
         return x
 
-    # ========================================================================
+    ###########################################################################
     # Core methods that use the above implemented functions
     #
     def energy(self, q, t=0.):
@@ -281,7 +282,7 @@ def hessian(self, q, t=0.):
         return hess.reshape((orig_shape[0],
                              orig_shape[0]) + orig_shape[1:]) * ret_unit
 
-    # ========================================================================
+    ###########################################################################
     # Convenience methods that make use the base methods
     #
     def acceleration(self, q, t=0.):
@@ -375,7 +376,7 @@ def circular_velocity(self, q, t=0.):
 
         return self.units.decompose(np.sqrt(r * np.abs(dPhi_dr)))
 
-    # ========================================================================
+    ###########################################################################
     # Python special methods
     #
     def __call__(self, q):
@@ -399,8 +400,8 @@ def __add__(self, other):
         if isinstance(other, CompositePotential):
             for k, v in self.items():
                 if k in new_pot:
-                    raise KeyError('Potential component "{}" already exists --'
-                                   'duplicate key provided in potential '
+                    raise KeyError(f'Potential component "{k}" already exists '
+                                   '-- duplicate key provided in potential '
                                    'addition')
                 new_pot[k] = v
 
@@ -410,7 +411,7 @@ def __add__(self, other):
 
         return new_pot
 
-    # ========================================================================
+    ###########################################################################
     # Convenience methods that do fancy things
     #
     def plot_contours(self, grid, filled=True, ax=None, labels=None,
@@ -462,8 +463,10 @@ def plot_contours(self, grid, filled=True, ax=None, labels=None,
 
         # if ndim > 2, don't know how to handle this!
         if ndim > 2:
-            raise ValueError("ndim > 2: you can only make contours on a 2D grid. For other "
-                             "dimensions, you have to specify values to slice.")
+            raise ValueError(
+                "ndim > 2: you can only make contours on a 2D grid. For other "
+                "dimensions, you have to specify values to slice."
+            )
 
         if ax is None:
             # default figsize
@@ -522,7 +525,9 @@ def plot_density_contours(self, grid, filled=True, ax=None, labels=None,
         Plot density contours. Computes the density on a grid
         (specified by the array `grid`).
 
-        .. warning:: Right now the grid input must be arrays and must already be in
+        .. warning::
+
+            For now, the grid input must be arrays and must already be in
             the unit system of the potential. Quantity support is coming...
 
         Parameters
@@ -565,8 +570,10 @@ def plot_density_contours(self, grid, filled=True, ax=None, labels=None,
 
         # if ndim > 2, don't know how to handle this!
         if ndim > 2:
-            raise ValueError("ndim > 2: you can only make contours on a 2D grid. For other "
-                             "dimensions, you have to specify values to slice.")
+            raise ValueError(
+                "ndim > 2: you can only make contours on a 2D grid. For other "
+                "dimensions, you have to specify values to slice."
+            )
 
         if ax is None:
             # default figsize
@@ -735,7 +742,7 @@ def replace_units(self, units, copy=True):
 
         return pot
 
-    # ========================================================================
+    ###########################################################################
     # Deprecated methods
     #
     def _value(self, q, t=0.):
@@ -747,6 +754,20 @@ def value(self, *args, **kwargs):
         warnings.warn("Use `energy()` instead.", DeprecationWarning)
         return self.energy(*args, **kwargs)
 
+    ###########################################################################
+    # Interoperability with other packages
+    #
+    def to_galpy_potential(self, ro=None, vo=None):
+        """Convert a Gala potential to a Galpy potential instance
+
+        Parameters
+        ----------
+        ro : quantity-like (optional)
+        vo : quantity-like (optional)
+        """
+        from .interop import gala_to_galpy_potential
+        return gala_to_galpy_potential(self, ro=ro, vo=vo)
+
 
 class CompositePotential(PotentialBase, OrderedDict):
     """
diff --git a/gala/potential/potential/interop.py b/gala/potential/potential/interop.py
new file mode 100644
index 00000000..62bc28bf
--- /dev/null
+++ b/gala/potential/potential/interop.py
@@ -0,0 +1,293 @@
+"""Interoperability with other dynamics packages"""
+
+import inspect
+
+from astropy.constants import G
+import astropy.units as u
+import numpy as np
+
+import gala.potential.potential.builtin as gp
+from gala.potential.potential.ccompositepotential import CCompositePotential
+from gala.potential.potential.core import CompositePotential
+from gala.units import galactic
+from gala.tests.optional_deps import HAS_GALPY
+
+__all__ = ['gala_to_galpy_potential', 'galpy_to_gala_potential']
+
+###############################################################################
+# Galpy interoperability
+#
+
+if HAS_GALPY:
+    from scipy.special import gamma
+    import galpy.potential as galpy_gp
+
+    def _powerlaw_amp_to_galpy(pars, ro, vo):
+        # I don't really remember why this is like this, but it might be related
+        # to the difference between GSL gamma and scipy gamma??
+        fac = ((1/(2*np.pi) * pars['r_c'].to_value(ro)**(pars['alpha'] - 3) /
+               (gamma(3/2 - pars['alpha']/2))))
+        amp = fac * (G * pars['m']).to_value(vo**2 * ro)
+        return amp
+
+    def _powerlaw_m_from_galpy(pars, ro, vo):
+        # See note above!
+        fac = ((1/(2*np.pi) * pars['rc']**(pars['alpha'] - 3) /
+               (gamma(3/2 - pars['alpha']/2))))
+        amp = pars['amp'] * vo**2 * ro
+        m = amp / G / fac
+        return m
+
+    # TODO: some potential conversions drop parameters. Might want to add an
+    # option for a custom validator function or something to raise warnings?
+    _gala_to_galpy = {
+        gp.HernquistPotential: (
+            galpy_gp.HernquistPotential, {
+                'a': 'c',
+                'amp': lambda pars, ro, vo: (G*2*pars['m']).to_value(ro*vo**2)
+            }
+        ),
+        gp.IsochronePotential: (
+            galpy_gp.IsochronePotential, {
+                'b': 'b'
+            }
+        ),
+        gp.JaffePotential: (
+            galpy_gp.JaffePotential, {
+                'a': 'c'
+            }
+        ),
+        gp.KeplerPotential: (galpy_gp.KeplerPotential, {}),
+        gp.KuzminPotential: (
+            galpy_gp.KuzminDiskPotential, {
+                'a': 'a',
+            }
+        ),
+        gp.LogarithmicPotential: (
+            galpy_gp.LogarithmicHaloPotential, {
+                'amp': lambda pars, ro, vo: pars['v_c'].to_value(vo)**2,
+                'core': 'r_h',
+                'q': 'q3'
+            }
+        ),
+        gp.LongMuraliBarPotential: (
+            galpy_gp.SoftenedNeedleBarPotential, {
+                'a': 'a',
+                'b': 'b',
+                'c': 'c',
+                'pa': 'alpha'
+            }
+        ),
+        gp.MiyamotoNagaiPotential: (
+            galpy_gp.MiyamotoNagaiPotential, {
+                'a': 'a',
+                'b': 'b'
+            }
+        ),
+        gp.NFWPotential: (
+            galpy_gp.TriaxialNFWPotential, {
+                'a': 'r_s',
+                'b': lambda pars, *_: pars['b'] / pars['a'],
+                'c': lambda pars, *_: pars['c'] / pars['a'],
+            }
+        ),
+        gp.PlummerPotential: (
+            galpy_gp.PlummerPotential, {
+                'b': 'b'
+            }
+        ),
+        gp.PowerLawCutoffPotential: (
+            galpy_gp.PowerSphericalPotentialwCutoff, {
+                'amp': _powerlaw_amp_to_galpy,
+                'rc': 'r_c',
+                'alpha': 'alpha'
+            }
+        ),
+    }
+
+    _galpy_to_gala = {}
+    for gala_cls, (galpy_cls, pars) in _gala_to_galpy.items():
+        galpy_pars = {v: k for k, v in pars.items()
+                      if isinstance(v, (str, int, float, np.ndarray))}
+        _galpy_to_gala[galpy_cls] = (gala_cls, galpy_pars)
+
+    # Special cases:
+    _galpy_to_gala[galpy_gp.HernquistPotential][1]['m'] = \
+        lambda pars, ro, vo: (pars['amp'] * ro * vo**2 / G / 2)
+
+    _galpy_to_gala[galpy_gp.LogarithmicHaloPotential][1]['v_c'] = \
+        lambda pars, ro, vo: np.sqrt(pars['amp'] * vo**2)
+
+    _galpy_to_gala[galpy_gp.TriaxialNFWPotential][1]['m'] = \
+        lambda pars, ro, vo: (
+            pars['amp'] * ro * vo**2 / G * 4*np.pi*pars['a']**3)
+    _galpy_to_gala[galpy_gp.TriaxialNFWPotential][1]['a'] = 1.
+    _galpy_to_gala[galpy_gp.TriaxialNFWPotential][1]['b'] = 'b'
+    _galpy_to_gala[galpy_gp.TriaxialNFWPotential][1]['c'] = 'c'
+
+    _galpy_to_gala[galpy_gp.PowerSphericalPotentialwCutoff][1]['m'] = \
+        _powerlaw_m_from_galpy
+
+    _galpy_to_gala[galpy_gp.NFWPotential] = (
+        gp.NFWPotential, {
+            'r_s': 'a',
+        }
+    )
+
+
+def _get_ro_vo(ro, vo):
+    # If not specified, get the default ro, vo from Galpy
+    if ro is None or vo is None:
+        from galpy.potential import Force
+        f = Force()
+
+        if ro is None:
+            ro = f._ro * u.kpc
+        if vo is None:
+            vo = f._vo * u.km/u.s
+
+    return u.Quantity(ro), u.Quantity(vo)
+
+
+def gala_to_galpy_potential(potential, ro=None, vo=None):
+
+    if not HAS_GALPY:
+        raise ImportError(
+            "Failed to import galpy.potential: Converting a potential to a "
+            "galpy potential requires galpy to be installed.")
+
+    ro, vo = _get_ro_vo(ro, vo)
+
+    if isinstance(potential, CompositePotential):
+        pot = []
+        for k in potential.keys():
+            pot.append(
+                gala_to_galpy_potential(potential[k], ro, vo))
+
+    else:
+        if potential.__class__ not in _gala_to_galpy:
+            raise TypeError(
+                f"Converting potential class {potential.__class__.__name__} "
+                "to galpy is currently not supported")
+
+        galpy_cls, converters = _gala_to_galpy[potential.__class__]
+        gala_pars = potential.parameters.copy()
+
+        galpy_pars = {}
+        if 'amp' not in converters and 'm' not in gala_pars:
+            raise ValueError("Gala potential has no mass parameter, so "
+                             "converting to a Galpy potential is currently "
+                             "not supported.")
+
+        converters.setdefault(
+            'amp', lambda pars, ro, vo: (G * pars['m']).to_value(ro * vo**2))
+
+        for galpy_par_name, conv in converters.items():
+            if isinstance(conv, str):
+                galpy_pars[galpy_par_name] = gala_pars[conv]
+            elif hasattr(conv, '__call__'):
+                galpy_pars[galpy_par_name] = conv(gala_pars, ro, vo)
+            elif isinstance(conv, (int, float, u.Quantity, np.ndarray)):
+                galpy_pars[galpy_par_name] = conv
+            else:
+                # TODO: invalid parameter??
+                print(f"FAIL: {galpy_par_name}, {conv}")
+
+            par = galpy_pars[galpy_par_name]
+            if hasattr(par, 'unit'):
+                if par.unit.physical_type == 'length':
+                    galpy_pars[galpy_par_name] = par.to_value(ro)
+                elif par.unit.physical_type == 'dimensionless':
+                    galpy_pars[galpy_par_name] = par.value
+                elif par.unit.physical_type == 'angle':
+                    galpy_pars[galpy_par_name] = par.to_value(u.rad)
+                else:
+                    # TODO: raise a warning here??
+                    galpy_pars[galpy_par_name] = par.value
+
+        pot = galpy_cls(**galpy_pars, ro=ro, vo=vo)
+
+    return pot
+
+
+def galpy_to_gala_potential(potential, ro=None, vo=None, units=galactic):
+
+    if not HAS_GALPY:
+        raise ImportError(
+            "Failed to import galpy.potential: Converting a potential to a "
+            "gala potential requires galpy to be installed.")
+
+    ro, vo = _get_ro_vo(ro, vo)
+
+    if potential._roSet:
+        ro = potential._ro * u.kpc
+    if potential._voSet:
+        vo = potential._vo * u.km/u.s
+
+    if isinstance(potential, list):
+        pot = CCompositePotential()
+        for i, sub_pot in enumerate(potential):
+            pot[str(i)] = galpy_to_gala_potential(sub_pot, ro, vo)
+
+    else:
+        if potential.__class__ not in _galpy_to_gala:
+            raise TypeError(
+                f"Converting galpy potential {potential.__class__.__name__} "
+                "to gala is currently not supported")
+
+        gala_cls, converters = _galpy_to_gala[potential.__class__]
+
+        exclude = ['self', 'normalize', 'ro', 'vo']
+        spec = inspect.getfullargspec(potential.__class__)
+        par_names = [arg for arg in spec.args if arg not in exclude]
+
+        # UGH!
+        galpy_pars = {}
+        for name in par_names:
+            galpy_pars[name] = getattr(potential,
+                                       '_' + name,
+                                       getattr(potential, name, None))
+
+        if isinstance(potential, galpy_gp.LogarithmicHaloPotential):
+            galpy_pars['core'] = np.sqrt(potential._core2)
+
+        elif isinstance(potential, galpy_gp.SoftenedNeedleBarPotential):
+            galpy_pars['c'] = np.sqrt(potential._c2)
+
+        if 'm' in inspect.getfullargspec(gala_cls).args:
+            converters.setdefault(
+                'm', lambda pars, ro, vo: pars['amp'] * ro * vo**2 / G
+            )
+
+        gala_pars = {}
+        for gala_par_name, conv in converters.items():
+            if isinstance(conv, str):
+                gala_pars[gala_par_name] = galpy_pars[conv]
+            elif hasattr(conv, '__call__'):
+                gala_pars[gala_par_name] = conv(galpy_pars, ro, vo)
+            elif isinstance(conv, (int, float, u.Quantity, np.ndarray)):
+                gala_pars[gala_par_name] = conv
+            else:
+                # TODO: invalid parameter??
+                print(f"FAIL: {gala_par_name}, {conv}")
+
+            if hasattr(gala_pars[gala_par_name], 'unit'):
+                continue
+
+            gala_par = gala_cls._parameters[gala_par_name]
+            if gala_par.physical_type == 'mass':
+                gala_pars[gala_par_name] = gala_pars[gala_par_name] * u.Msun
+            elif gala_par.physical_type == 'length':
+                gala_pars[gala_par_name] = gala_pars[gala_par_name] * ro
+            elif gala_par.physical_type == 'speed':
+                gala_pars[gala_par_name] = gala_pars[gala_par_name] * vo
+            elif gala_par.physical_type == 'angle':
+                gala_pars[gala_par_name] = gala_pars[gala_par_name] * u.radian
+            elif gala_par.physical_type == 'dimensionless':
+                pass
+            else:
+                print("TODO")
+
+        pot = gala_cls(**gala_pars, units=units)
+
+    return pot
diff --git a/gala/potential/potential/util.py b/gala/potential/potential/util.py
index 34eefdbd..c560519f 100644
--- a/gala/potential/potential/util.py
+++ b/gala/potential/potential/util.py
@@ -83,13 +83,15 @@ def from_equation(expr, vars, pars, name=None, hessian=False):
     ndim = len(vars)
 
     # Energy / value
-    energyfunc = lambdify(vars + pars, expr, dummify=False, modules='numpy')
+    energyfunc = lambdify(vars + pars, expr, dummify=False,
+                          modules=['numpy', 'sympy'])
 
     # Gradient
     gradfuncs = []
     for var in vars:
         gradfuncs.append(lambdify(vars + pars, sympy.diff(expr, var),
-                                  dummify=False, modules='numpy'))
+                                  dummify=False,
+                                  modules=['numpy', 'sympy']))
 
     parameters = {}
     for _name in par_names:
@@ -133,7 +135,8 @@ def _gradient(self, w, t=0.):
             for var2 in vars:
                 hessfuncs.append(lambdify(vars + pars,
                                           sympy.diff(expr, var1, var2),
-                                          dummify=False, modules='numpy'))
+                                          dummify=False,
+                                          modules=['numpy', 'sympy']))
 
         def _hessian(self, w, t):
             kw = self.parameters.copy()
@@ -212,11 +215,12 @@ def wrapper(cls, *func_args, **func_kwargs):
                 raise ImportError("Converting to a latex expression requires "
                                   "the sympy package to be installed")
 
-            _var = sy.symbols(self.var, seq=True)
+            _var = sy.symbols(self.var, seq=True, real=True)
             _var = {v.name: v for v in _var}
 
             if cls._parameters:
-                par = sy.symbols(' '.join(cls._parameters.keys()), seq=True)
+                par = sy.symbols(' '.join(cls._parameters.keys()),
+                                 seq=True, real=True)
                 par = {v.name: v for v in par}
             else:
                 par = {}
diff --git a/tox.ini b/tox.ini
index 17bc92a7..ec226449 100644
--- a/tox.ini
+++ b/tox.ini
@@ -69,7 +69,7 @@ commands =
     cov: pytest --pyargs gala {toxinidir}/docs --cov gala --cov-config={toxinidir}/setup.cfg {posargs} --cov-report=xml:{toxinidir}/coverage.xml --durations=16
 
 # Runs pip install -e . instead of building an sdist and installing
-usedevelop = True
+usedevelop = False
 
 [testenv:local_test]
 changedir = .tmp/{envname}


Test patch:
diff --git a/gala/dynamics/nbody/tests/test_nbody.py b/gala/dynamics/nbody/tests/test_nbody.py
index 0729ff5f..29e21800 100644
--- a/gala/dynamics/nbody/tests/test_nbody.py
+++ b/gala/dynamics/nbody/tests/test_nbody.py
@@ -5,7 +5,7 @@
 
 # Custom
 from ....potential import (NullPotential, NFWPotential,
-                           HernquistPotential, KuzminPotential,
+                           HernquistPotential,
                            ConstantRotatingFrame, StaticFrame)
 from ....dynamics import PhaseSpacePosition, combine
 from ....units import UnitSystem, galactic
@@ -73,11 +73,6 @@ def test_directnbody_init(self):
         with pytest.raises(ValueError):
             DirectNBody(self.w0, particle_potentials=[None, None])
 
-        py_ext_pot = KuzminPotential(m=1e10*u.Msun, a=0.5*u.kpc, units=galactic)
-        with pytest.raises(ValueError):
-            DirectNBody(self.w0, particle_potentials=self.particle_potentials,
-                        external_potential=py_ext_pot)
-
     def test_directnbody_integrate(self):
         # TODO: this is really a unit test, but we should have some functional tests
         # that check that the orbit integration is making sense!
diff --git a/gala/dynamics/tests/test_core.py b/gala/dynamics/tests/test_core.py
index 4029ce49..7c0b8189 100644
--- a/gala/dynamics/tests/test_core.py
+++ b/gala/dynamics/tests/test_core.py
@@ -15,7 +15,7 @@
 from ...potential import Hamiltonian, HernquistPotential
 from ...potential.frame import StaticFrame, ConstantRotatingFrame
 from ...units import galactic, solarsystem
-from gala.tests.optional_deps import HAS_H5PY, HAS_GALPY
+from gala.tests.optional_deps import HAS_H5PY
 
 
 def test_initialize():
@@ -403,17 +403,3 @@ def test_io(tmpdir, obj):
     assert u.allclose(obj.xyz, obj2.xyz)
     assert u.allclose(obj.v_xyz, obj2.v_xyz)
     assert obj.frame == obj2.frame
-
-
-@pytest.mark.parametrize('obj', [
-    PhaseSpacePosition([1, 2, 3.]*u.kpc, [1, 2, 3.]*u.km/u.s),
-    PhaseSpacePosition([1, 2, 3.]*u.kpc, [1, 2, 3.]*u.km/u.s,
-                       StaticFrame(units=galactic)),
-])
-@pytest.mark.skipif(not HAS_GALPY,
-                    reason="requires galpy to run this test")
-def test_to_galpy(obj):
-    obj.to_galpy_orbit()
-    obj.to_galpy_orbit(ro=8*u.kpc)
-    obj.to_galpy_orbit(vo=220*u.km/u.s)
-    obj.to_galpy_orbit(ro=8*u.kpc, vo=220*u.km/u.s)
diff --git a/gala/dynamics/tests/test_orbit.py b/gala/dynamics/tests/test_orbit.py
index 689c2174..be12d876 100644
--- a/gala/dynamics/tests/test_orbit.py
+++ b/gala/dynamics/tests/test_orbit.py
@@ -592,3 +592,17 @@ def test_orbit_to_galpy(obj):
     o2 = obj.to_galpy_orbit(ro=8*u.kpc)  # noqa
     o3 = obj.to_galpy_orbit(vo=220*u.km/u.s)  # noqa
     o4 = obj.to_galpy_orbit(ro=8*u.kpc, vo=220*u.km/u.s)  # noqa
+
+
+@pytest.mark.skipif(not HAS_GALPY,
+                    reason="requires galpy to run this test")
+def test_orbit_from_galpy():
+    import galpy.orbit as galpy_o
+    import galpy.potential as galpy_p
+    mp = galpy_p.MiyamotoNagaiPotential(a=0.5, b=0.0375, amp=1., normalize=1.)
+    galpy_orbit = galpy_o.Orbit([1., 0.1, 1.1, 0., 0.1, 1.])
+    ts = np.linspace(0, 100, 10000)
+    galpy_orbit.integrate(ts, mp, method='odeint')
+    gala_orbit = Orbit.from_galpy_orbit(galpy_orbit)
+
+    assert len(gala_orbit.t) == len(ts)
diff --git a/gala/potential/potential/tests/helpers.py b/gala/potential/potential/tests/helpers.py
index 0e0933e5..8c5e0690 100644
--- a/gala/potential/potential/tests/helpers.py
+++ b/gala/potential/potential/tests/helpers.py
@@ -28,13 +28,16 @@ def wraps(a):
     return derivative(wraps, point[dim_ix], **kwargs)
 
 
-class PotentialTestBase(object):
+class PotentialTestBase:
     name = None
     potential = None  # MUST SET THIS
     frame = None
     tol = 1E-5
     show_plots = False
 
+    sympy_hessian = True
+    sympy_density = True
+
     @classmethod
     def setup_class(cls):
         if cls.name is None:
@@ -292,9 +295,12 @@ def test_pickle(self, tmpdir):
     @pytest.mark.skipif(not HAS_SYMPY,
                         reason="requires sympy to run this test")
     def test_against_sympy(self):
+        # TODO: should really split this into separate tests for each check...
+
         import sympy as sy
+        from sympy import Q
 
-        # compare Gala gradient and hessian to sympy values
+        # compare Gala gradient, hessian, and density to sympy values
 
         pot = self.potential
         Phi, v, p = pot.to_sympy()
@@ -304,21 +310,37 @@ def test_against_sympy(self):
         def lowergamma(a, x):  # noqa
             # Differences between scipy and sympy lower gamma
             return gammainc(a, x) * gamma(a)
-        modules = ['numpy',
-                   {'atan': np.arctan,
-                    'lowergamma': lowergamma,
-                    'gamma': gamma}]
-
-        e_func = sy.lambdify(list(p.values()) + list(v.values()), Phi,
-                             modules=modules)
+        modules = [
+            'numpy',
+            {
+                'atan': np.arctan,
+                'lowergamma': lowergamma,
+                'gamma': gamma,
+                're': np.real,
+                'im': np.imag
+            },
+            'sympy'
+        ]
+
+        vars_ = list(p.values()) + list(v.values())
+        assums = np.bitwise_and.reduce([Q.real(x) for x in vars_])
+        Phi = sy.refine(Phi, assums)
+        e_func = sy.lambdify(vars_, Phi, modules=modules)
+
+        if self.sympy_density:
+            dens_tmp = sum([sy.diff(Phi, var, 2)
+                            for var in v.values()]) / (4 * sy.pi * p['G'])
+            dens_tmp = sy.refine(dens_tmp, assums)
+            dens_func = sy.lambdify(vars_, dens_tmp, modules=modules)
 
         grad = sy.derive_by_array(Phi, list(v.values()))
-        grad_func = sy.lambdify(list(p.values()) + list(v.values()), grad,
-                                modules=modules)
+        grad = sy.refine(grad, assums)
+        grad_func = sy.lambdify(vars_, grad, modules=modules)
 
-        Hess = sy.hessian(Phi, list(v.values()))
-        Hess_func = sy.lambdify(list(p.values()) + list(v.values()), Hess,
-                                modules=modules)
+        if self.sympy_hessian:
+            Hess = sy.hessian(Phi, list(v.values()))
+            Hess = sy.refine(Hess, assums)
+            Hess_func = sy.lambdify(vars_, Hess, modules=modules)
 
         # Make a dict of potential parameter values without units:
         par_vals = {}
@@ -332,20 +354,37 @@ def lowergamma(a, x):  # noqa
         f_gala = pot.energy(trial_x).value
         f_sympy = e_func(G=pot.G, **par_vals, **x_dict)
         e_close = np.allclose(f_gala, f_sympy)
+        test_cases = [e_close]
+        vals = [(f_gala, f_sympy)]
+
+        if self.sympy_density:
+            d_gala = pot.density(trial_x).value
+            d_sympy = dens_func(G=pot.G, **par_vals, **x_dict)
+            d_close = np.allclose(d_gala, d_sympy)
+            test_cases.append(d_close)
+            vals.append((d_gala, d_sympy))
 
         G_gala = pot.gradient(trial_x).value
         G_sympy = grad_func(G=pot.G, **par_vals, **x_dict)
         g_close = np.allclose(G_gala, G_sympy)
-
-        H_gala = pot.hessian(trial_x).value
-        H_sympy = Hess_func(G=pot.G, **par_vals, **x_dict)
-        h_close = np.allclose(H_gala, H_sympy)
-
-        if not all([e_close, g_close, h_close]):
-            print(f'{pot}: energy {e_close}, gradient {g_close}, '
-                  f'hessian {h_close}')
-
-        assert all([e_close, g_close, h_close])
+        test_cases.append(g_close)
+        vals.append((G_gala, G_sympy))
+
+        if self.sympy_hessian:
+            H_gala = pot.hessian(trial_x).value
+            H_sympy = Hess_func(G=pot.G, **par_vals, **x_dict)
+            h_close = np.allclose(H_gala, H_sympy)
+            test_cases.append(h_close)
+            vals.append((H_gala, H_sympy))
+
+        if not all(test_cases):
+            names = ['energy', 'density', 'gradient', 'hessian']
+            for name, (val1, val2), test in zip(names, vals, test_cases):
+                if not test:
+                    print(trial_x)
+                    print(f'{pot}: {name}\nGala:{val1}\nSympy:{val2}')
+
+        assert all(test_cases)
 
     def test_regression_165(self):
         if self.potential.ndim == 1:
diff --git a/gala/potential/potential/tests/test_against_galpy.py b/gala/potential/potential/tests/test_against_galpy.py
deleted file mode 100644
index fcd6becd..00000000
--- a/gala/potential/potential/tests/test_against_galpy.py
+++ /dev/null
@@ -1,140 +0,0 @@
-"""Test some builtin potentials against galpy"""
-
-# Third-party
-import numpy as np
-from astropy.constants import G
-import astropy.units as u
-import pytest
-from scipy.special import gamma
-
-# This project
-from ...._cconfig import GSL_ENABLED
-from ....units import galactic
-from ..builtin import (KeplerPotential, MiyamotoNagaiPotential,
-                       NFWPotential, PowerLawCutoffPotential,
-                       BovyMWPotential2014)
-from gala.tests.optional_deps import HAS_GALPY
-
-# Set to arbitrary values for testing
-ro = 8.1 * u.kpc
-vo = 240 * u.km/u.s
-ntest = 128
-
-
-def helper(gala_pot, galpy_pot):
-    Rs = np.random.uniform(1, 15, size=ntest) * u.kpc
-    zs = np.random.uniform(1, 15, size=ntest) * u.kpc
-
-    xyz = np.zeros((3, Rs.size)) * u.kpc
-    xyz[0] = Rs
-
-    assert np.allclose(gala_pot.circular_velocity(xyz).to_value(u.km/u.s),
-                       galpy_pot.vcirc(R=Rs.to_value(ro)))
-
-    xyz[2] = zs
-    assert np.allclose(gala_pot.density(xyz).to_value(u.Msun/u.pc**3),
-                       galpy_pot.dens(R=Rs.to_value(ro), z=zs.to_value(ro)))
-
-    assert np.allclose(gala_pot.energy(xyz).to_value((u.km / u.s)**2),
-                       galpy_pot(R=Rs.to_value(ro), z=zs.to_value(ro)))
-
-    assert np.allclose(gala_pot.gradient(xyz).to_value((u.km/u.s) * u.pc/u.Myr / u.pc)[2],
-                       -galpy_pot.zforce(R=Rs.to_value(ro), z=zs.to_value(ro)))
-
-
-@pytest.mark.skipif(not HAS_GALPY,
-                    reason="requires galpy to run this test")
-def test_kepler():
-    from galpy.potential import KeplerPotential as BovyKeplerPotential
-
-    M = 5e10 * u.Msun
-    gala_pot = KeplerPotential(m=M, units=galactic)
-
-    amp = (G*M).to_value(vo**2 * ro)
-    bovy_pot = BovyKeplerPotential(amp=amp, ro=ro, vo=vo)
-
-    helper(gala_pot, bovy_pot)
-
-
-@pytest.mark.skipif(not HAS_GALPY,
-                    reason="requires galpy to run this test")
-def test_miyamoto():
-    from galpy.potential import MiyamotoNagaiPotential as BovyMiyamotoNagaiPotential
-
-    M = 5e10 * u.Msun
-    gala_pot = MiyamotoNagaiPotential(m=M, a=3.5*u.kpc, b=300*u.pc,
-                                      units=galactic)
-
-    amp = (G*M).to_value(vo**2 * ro)
-    a = gala_pot.parameters['a'].to_value(ro)
-    b = gala_pot.parameters['b'].to_value(ro)
-    bovy_pot = BovyMiyamotoNagaiPotential(amp=amp, a=a, b=b, ro=ro, vo=vo)
-
-    helper(gala_pot, bovy_pot)
-
-
-@pytest.mark.skipif(not HAS_GALPY,
-                    reason="requires galpy to run this test")
-def test_nfw():
-    from galpy.potential import NFWPotential as BovyNFWPotential
-
-    M = 6.5854e10 * u.Msun
-    gala_pot = NFWPotential(m=M, r_s=15*u.kpc, units=galactic)
-    amp = (G*M).to_value(vo**2 * ro)
-    a = gala_pot.parameters['r_s'].to_value(ro)
-    bovy_pot = BovyNFWPotential(amp, a=a, ro=ro, vo=vo)
-
-    helper(gala_pot, bovy_pot)
-
-
-@pytest.mark.skipif(not HAS_GALPY or not GSL_ENABLED,
-                    reason="requires galpy and GSL to run this test")
-def test_powerlawcutoff():
-    from galpy.potential import PowerSphericalPotentialwCutoff
-
-    M = 1.5854e10 * u.Msun
-    alpha = 1.75
-    gala_pot = PowerLawCutoffPotential(m=M, alpha=alpha, r_c=15*u.kpc,
-                                       units=galactic)
-    r_c = gala_pot.parameters['r_c']
-    amp = (G*M).to_value(vo**2 * ro) * ((1/(2*np.pi) * r_c.to_value(ro)**(alpha - 3) /
-                                         (gamma(3/2 - alpha/2))))
-    bovy_pot = PowerSphericalPotentialwCutoff(amp,
-                                              alpha=gala_pot.parameters['alpha'].value,
-                                              rc=r_c.to_value(ro),
-                                              ro=ro, vo=vo)
-
-    helper(gala_pot, bovy_pot)
-
-
-@pytest.mark.skipif(not HAS_GALPY or not GSL_ENABLED,
-                    reason="requires galpy and GSL to run this test")
-def test_mwpotential2014():
-    from galpy.potential import (MWPotential2014,
-                                 evaluateDensities,
-                                 evaluatezforces,
-                                 evaluatePotentials)
-
-    # Here these have to be default:
-    ro = 8 * u.kpc
-    vo = 220 * u.km/u.s
-
-    gala_pot = BovyMWPotential2014()
-    bovy_pot = MWPotential2014
-    for x in bovy_pot:
-        x.turn_physical_on()
-
-    Rs = np.random.uniform(1, 15, size=ntest) * u.kpc
-    zs = np.random.uniform(1, 15, size=ntest) * u.kpc
-
-    xyz = np.zeros((3, Rs.size)) * u.kpc
-    xyz[0] = Rs
-    xyz[2] = zs
-    assert np.allclose(gala_pot.density(xyz).to_value(u.Msun/u.pc**3),
-                       evaluateDensities(bovy_pot, R=Rs.to_value(ro), z=zs.to_value(ro)))
-
-    assert np.allclose(gala_pot.energy(xyz).to_value((u.km / u.s)**2),
-                       evaluatePotentials(bovy_pot, R=Rs.to_value(ro), z=zs.to_value(ro)))
-
-    assert np.allclose(gala_pot.gradient(xyz).to_value((u.km/u.s) * u.pc/u.Myr / u.pc)[2],
-                       -evaluatezforces(bovy_pot, R=Rs.to_value(ro), z=zs.to_value(ro)))
diff --git a/gala/potential/potential/tests/test_all_builtin.py b/gala/potential/potential/tests/test_all_builtin.py
index e663aeca..51dab931 100644
--- a/gala/potential/potential/tests/test_all_builtin.py
+++ b/gala/potential/potential/tests/test_all_builtin.py
@@ -26,6 +26,7 @@
 class TestHarmonicOscillator1D(PotentialTestBase):
     potential = p.HarmonicOscillatorPotential(omega=1.)
     w0 = [1., 0.1]
+    sympy_density = False
 
     def test_plot(self):
         # Skip for now because contour plotting assumes 3D
@@ -35,6 +36,7 @@ def test_plot(self):
 class TestHarmonicOscillator2D(PotentialTestBase):
     potential = p.HarmonicOscillatorPotential(omega=[1., 2])
     w0 = [1., 0.5, 0., 0.1]
+    sympy_density = False
 
     def test_plot(self):
         # Skip for now because contour plotting assumes 3D
@@ -88,6 +90,7 @@ def test_against_sympy(self):
 class TestHenonHeiles(PotentialTestBase):
     potential = p.HenonHeilesPotential()
     w0 = [1., 0., 0., 2*np.pi]
+    sympy_density = False
 
 
 class TestKepler(PotentialTestBase):
@@ -176,6 +179,13 @@ class TestSatoh(PotentialTestBase):
     w0 = [8., 0., 0., 0., 0.22, 0.1]
 
 
+class TestKuzmin(PotentialTestBase):
+    potential = p.KuzminPotential(units=galactic, m=1.E11, a=3.5)
+    w0 = [8., 0., 0., 0., 0.22, 0.1]
+    sympy_hessian = False
+    sympy_density = False
+
+
 class TestStone(PotentialTestBase):
     potential = p.StonePotential(units=galactic, m=1E11, r_c=0.1, r_h=10.)
     w0 = [8., 0., 0., 0., 0.18, 0.1]
@@ -186,6 +196,7 @@ class TestStone(PotentialTestBase):
 class TestPowerLawCutoff(PotentialTestBase):
     w0 = [8., 0., 0., 0., 0.1, 0.1]
     atol = 1e-3
+    sympy_density = False  # weird underflow issues??
 
     def setup(self):
         self.potential = p.PowerLawCutoffPotential(units=galactic,
@@ -202,6 +213,7 @@ class TestTriaxialNFW(PotentialTestBase):
     potential = p.NFWPotential(units=galactic, m=1E11, r_s=12.,
                                a=1., b=0.95, c=0.9)
     w0 = [19.0, 2.7, -6.9, 0.0352238, -0.03579493, 0.075]
+    sympy_density = False  # not defined
 
 
 class TestSphericalNFWFromCircVel(PotentialTestBase):
@@ -261,6 +273,7 @@ class TestNFW(PotentialTestBase):
     potential = p.NFWPotential(m=6E11*u.Msun, r_s=20*u.kpc, a=1., b=0.9, c=0.75,
                                units=galactic)
     w0 = [19.0, 2.7, -0.9, 0.00352238, -0.15134, 0.0075]
+    sympy_density = False  # like triaxial case
 
     def test_compare(self):
 
diff --git a/gala/potential/potential/tests/test_interop_galpy.py b/gala/potential/potential/tests/test_interop_galpy.py
new file mode 100644
index 00000000..264b4957
--- /dev/null
+++ b/gala/potential/potential/tests/test_interop_galpy.py
@@ -0,0 +1,168 @@
+"""
+Test converting the builtin Potential classes to other packages
+"""
+
+# Third-party
+from astropy.coordinates import CylindricalRepresentation
+import astropy.units as u
+import numpy as np
+import pytest
+
+# This project
+import gala.potential as gp
+from gala.units import galactic
+from gala.tests.optional_deps import HAS_GALPY
+from gala.potential.potential.interop import galpy_to_gala_potential
+
+# Set these globally!
+ro = 8.122 * u.kpc
+vo = 245 * u.km/u.s
+
+if HAS_GALPY:
+    import galpy.potential as galpy_gp
+
+    from gala.potential.potential.interop import (
+        _gala_to_galpy,
+        _galpy_to_gala
+    )
+
+
+def pytest_generate_tests(metafunc):
+    # Some magic, semi-random numbers below!
+    gala_pots = []
+    galpy_pots = []
+
+    if not HAS_GALPY:
+        return
+
+    # Test the Gala -> Galpy direction
+    for Potential in _gala_to_galpy.keys():
+        init = {}
+        len_scale = 1.
+        for k, par in Potential._parameters.items():
+            if k == 'm':
+                val = 1.43e10 * u.Msun
+            elif par.physical_type == 'length':
+                val = 5.12 * u.kpc * len_scale
+                len_scale *= 0.5
+            elif par.physical_type == 'dimensionless':
+                val = 1.
+            elif par.physical_type == 'speed':
+                val = 201.41 * u.km/u.s
+            else:
+                continue
+
+            init[k] = val
+
+        pot = Potential(**init, units=galactic)
+        galpy_pot = pot.to_galpy_potential(ro=ro, vo=vo)
+
+        gala_pots.append(pot)
+        galpy_pots.append(galpy_pot)
+
+    # Make a composite potential too:
+    gala_pots.append(gala_pots[0] + gala_pots[1])
+    galpy_pots.append([galpy_pots[0], galpy_pots[1]])
+
+    # Test the Galpy -> Gala direction
+    for Potential in _galpy_to_gala.keys():
+        galpy_pot = Potential(ro=ro, vo=vo)  # use defaults
+        pot = galpy_to_gala_potential(galpy_pot, ro=ro, vo=vo)
+
+        gala_pots.append(pot)
+        galpy_pots.append(galpy_pot)
+
+    test_names = [f'{g1.__class__.__name__}:{g2.__class__.__name__}'
+                  for g1, g2 in zip(gala_pots, galpy_pots)]
+
+    metafunc.parametrize(['gala_pot', 'galpy_pot'],
+                         list(zip(gala_pots, galpy_pots)),
+                         ids=test_names)
+
+
+@pytest.mark.skipif(not HAS_GALPY,
+                    reason="must have galpy installed to run these tests")
+class TestGalpy:
+
+    def setup(self):
+        # Test points:
+        rng = np.random.default_rng(42)
+        ntest = 128
+
+        Rs = rng.uniform(1, 15, size=ntest) * u.kpc
+        phis = rng.uniform(0, 2*np.pi, size=ntest) * u.radian
+        zs = rng.uniform(1, 15, size=ntest) * u.kpc
+
+        cyl = CylindricalRepresentation(Rs, phis, zs)
+        xyz = cyl.to_cartesian().xyz
+
+        self.Rs = Rs.to_value(ro)
+        self.phis = phis.to_value(u.rad)
+        self.zs = zs.to_value(ro)
+        self.Rpz_iter = list(zip(self.Rs, self.phis, self.zs))
+
+        self.xyz = xyz
+
+        Jac = np.zeros((len(cyl), 3, 3))
+        Jac[:, 0, 0] = xyz[0] / cyl.rho
+        Jac[:, 0, 1] = xyz[1] / cyl.rho
+        Jac[:, 1, 0] = (-xyz[1] / cyl.rho**2).to_value(1 / ro)
+        Jac[:, 1, 1] = (xyz[0] / cyl.rho**2).to_value(1 / ro)
+        Jac[:, 2, 2] = 1.
+        self.Jac = Jac
+
+    def test_density(self, gala_pot, galpy_pot):
+        if isinstance(gala_pot, gp.LogarithmicPotential):
+            pytest.skip()
+
+        gala_val = gala_pot.density(self.xyz).to_value(u.Msun / u.pc**3)
+        galpy_val = np.array([galpy_gp.evaluateDensities(galpy_pot,
+                                                         R=RR, z=zz, phi=pp)
+                              for RR, pp, zz in self.Rpz_iter])
+        assert np.allclose(gala_val, galpy_val)
+
+    def test_energy(self, gala_pot, galpy_pot):
+        gala_val = gala_pot.energy(self.xyz).to_value(u.km**2 / u.s**2)
+        galpy_val = np.array([galpy_gp.evaluatePotentials(galpy_pot,
+                                                          R=RR, z=zz, phi=pp)
+                              for RR, pp, zz in self.Rpz_iter])
+
+        if isinstance(gala_pot, gp.LogarithmicPotential):
+            # Logarithms are weird
+            gala_val -= (0.5 * gala_pot.parameters['v_c']**2 *
+                         np.log(ro.value**2)).to_value((u.km / u.s)**2)
+
+        assert np.allclose(gala_val, galpy_val)
+
+    def test_gradient(self, gala_pot, galpy_pot):
+        gala_grad = gala_pot.gradient(self.xyz)
+        gala_grad = gala_grad.to_value(u.km/u.s/u.Myr)
+
+        galpy_dR = np.array([-galpy_gp.evaluateRforces(galpy_pot,
+                                                       R=RR, z=zz, phi=pp)
+                            for RR, pp, zz in self.Rpz_iter])
+        galpy_dp = np.array([-galpy_gp.evaluatephiforces(galpy_pot,
+                                                         R=RR, z=zz, phi=pp)
+                            for RR, pp, zz in self.Rpz_iter])
+        galpy_dz = np.array([-galpy_gp.evaluatezforces(galpy_pot,
+                                                       R=RR, z=zz, phi=pp)
+                            for RR, pp, zz in self.Rpz_iter])
+        galpy_dRpz = np.stack((galpy_dR, galpy_dp, galpy_dz),
+                              axis=1)
+
+        galpy_grad = np.einsum('nij,ni->nj', self.Jac, galpy_dRpz).T
+
+        assert np.allclose(gala_grad, galpy_grad)
+
+    def test_vcirc(self, gala_pot, galpy_pot):
+        tmp = self.xyz.copy()
+        tmp[2] = 0.
+
+        if (not hasattr(galpy_pot, 'vcirc')
+                or isinstance(gala_pot, gp.LongMuraliBarPotential)):
+            pytest.skip()
+
+        gala_vcirc = gala_pot.circular_velocity(tmp).to_value(u.km/u.s)
+        galpy_vcirc = np.array([galpy_pot.vcirc(R=RR)
+                                for RR, *_ in self.Rpz_iter])
+        assert np.allclose(gala_vcirc, galpy_vcirc)