Repository: jobovy/galpy

Issue:
Handle unbound orbits in actionAngle-related functions in new Orbit class
In the new Orbit class, pure Python actionAngle-related functions fail for unbound orbits and cause the entire function to fail when Orbit contains multiple objects. Need to catch the unbound orbits without raising an error (see also #371).

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/jobovy/galpy/pull/712?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Jo+Bovy) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 99.90%. Comparing base [(`6a74587`)](https://app.codecov.io/gh/jobovy/galpy/commit/6a745879491a6fb0fe3bf9ed3f5d3036abf8c6bc?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Jo+Bovy) to head [(`84e9181`)](https://app.codecov.io/gh/jobovy/galpy/commit/84e918132532ead5871b2cfb3f35fa396cfd1103?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Jo+Bovy).

<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##             main     #712   +/-   ##
=======================================
  Coverage   99.90%   99.90%           
=======================================
  Files         202      202           
  Lines       29698    29722   +24     
  Branches      606      607    +1     
=======================================
+ Hits        29671    29695   +24     
  Misses         27       27           
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/jobovy/galpy/pull/712?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Jo+Bovy).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Jo+Bovy).


Patch:
diff --git a/galpy/actionAngle/actionAngleStaeckel.py b/galpy/actionAngle/actionAngleStaeckel.py
index d1235644..a3732d6c 100644
--- a/galpy/actionAngle/actionAngleStaeckel.py
+++ b/galpy/actionAngle/actionAngleStaeckel.py
@@ -214,7 +214,13 @@ class actionAngleStaeckel(actionAngle):
             else:
                 # Set up the actionAngleStaeckelSingle object
                 aASingle = actionAngleStaeckelSingle(
-                    R[0], vR[0], vT[0], z[0], vz[0], pot=self._pot, delta=delta
+                    R[0],
+                    vR[0],
+                    vT[0],
+                    z[0],
+                    vz[0],
+                    pot=self._pot,
+                    delta=delta[0] if hasattr(delta, "__len__") else delta,
                 )
                 return (
                     numpy.atleast_1d(aASingle.JR(**copy.copy(kwargs))),
diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py
index 2c19efcc..976c81c5 100644
--- a/galpy/orbit/Orbits.py
+++ b/galpy/orbit/Orbits.py
@@ -38,6 +38,7 @@ from ..potential import (
     toPlanarPotential,
 )
 from ..potential.DissipativeForce import _isDissipative
+from ..potential.plotEscapecurve import _INF
 from ..potential.Potential import _check_c
 from ..util import conversion, coords, galpyWarning, galpyWarningVerbose, plot
 from ..util._optional_deps import (
@@ -2601,6 +2602,26 @@ class Orbit:
             self._aA = actionAngle.actionAngleSpherical(pot=self._aAPot, **kwargs)
         return None
 
+    def _unbound_indx_and_aAkwargs(self):
+        """Internal function to compute the index of unbound orbits"""
+        if self.dim() == 2:
+            Einf = evaluateplanarPotentials(
+                toPlanarPotential(self._aAPot), _INF, use_physical=False
+            )
+        elif self.dim() == 3:
+            Einf = evaluatePotentials(self._aAPot, _INF, 0.0, use_physical=False)
+        if numpy.isnan(Einf):
+            Einf = numpy.inf  # Just try to proceed as best as possible, don't make assumptions about the potential
+        indx = self.E(pot=self._aAPot, use_physical=False, dontreshape=True) < Einf
+        if hasattr(self._aA, "_delta"):
+            if hasattr(self._aA._delta, "__len__"):
+                aAkwargs = {"delta": self._aA._delta[indx]}
+            else:
+                aAkwargs = {"delta": self._aA._delta}
+        else:
+            aAkwargs = {}
+        return indx, aAkwargs
+
     def _setup_EccZmaxRperiRap(self, pot=None, **kwargs):
         """Internal function to compute e,zmax,rperi,rap and cache it for reuse"""
         self._setupaA(pot=pot, **kwargs)
@@ -2619,19 +2640,34 @@ class Orbit:
             tz = numpy.zeros(self.size)
             tvz = numpy.zeros(self.size)
         # self.dim() == 1 error caught by _setupaA
+        # Exclude unbound orbits (aAkwargs deals with delta processing)
+        indx, aAkwargs = self._unbound_indx_and_aAkwargs()
         (
             self._aA_ecc,
             self._aA_zmax,
             self._aA_rperi,
             self._aA_rap,
-        ) = self._aA.EccZmaxRperiRap(
-            self.R(use_physical=False, dontreshape=True),
-            self.vR(use_physical=False, dontreshape=True),
-            self.vT(use_physical=False, dontreshape=True),
-            tz,
-            tvz,
-            use_physical=False,
+        ) = (
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
         )
+        if numpy.sum(indx) > 0:
+            (
+                self._aA_ecc[indx],
+                self._aA_zmax[indx],
+                self._aA_rperi[indx],
+                self._aA_rap[indx],
+            ) = self._aA.EccZmaxRperiRap(
+                self.R(use_physical=False, dontreshape=True)[indx],
+                self.vR(use_physical=False, dontreshape=True)[indx],
+                self.vT(use_physical=False, dontreshape=True)[indx],
+                tz[indx],
+                tvz[indx],
+                use_physical=False,
+                **aAkwargs,
+            )
         return None
 
     def _setup_actionsFreqsAngles(self, pot=None, **kwargs):
@@ -2652,6 +2688,8 @@ class Orbit:
             tz = numpy.zeros(self.size)
             tvz = numpy.zeros(self.size)
         # self.dim() == 1 error caught by _setupaA
+        # Exclude unbound orbits (aAkwargs deals with delta processing)
+        indx, aAkwargs = self._unbound_indx_and_aAkwargs()
         (
             self._aA_jr,
             self._aA_jp,
@@ -2662,15 +2700,38 @@ class Orbit:
             self._aA_wr,
             self._aA_wp,
             self._aA_wz,
-        ) = self._aA.actionsFreqsAngles(
-            self.R(use_physical=False, dontreshape=True),
-            self.vR(use_physical=False, dontreshape=True),
-            self.vT(use_physical=False, dontreshape=True),
-            tz,
-            tvz,
-            self.phi(use_physical=False, dontreshape=True),
-            use_physical=False,
+        ) = (
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
         )
+        if numpy.sum(indx) > 0:
+            (
+                self._aA_jr[indx],
+                self._aA_jp[indx],
+                self._aA_jz[indx],
+                self._aA_Or[indx],
+                self._aA_Op[indx],
+                self._aA_Oz[indx],
+                self._aA_wr[indx],
+                self._aA_wp[indx],
+                self._aA_wz[indx],
+            ) = self._aA.actionsFreqsAngles(
+                self.R(use_physical=False, dontreshape=True)[indx],
+                self.vR(use_physical=False, dontreshape=True)[indx],
+                self.vT(use_physical=False, dontreshape=True)[indx],
+                tz[indx],
+                tvz[indx],
+                self.phi(use_physical=False, dontreshape=True)[indx],
+                use_physical=False,
+                **aAkwargs,
+            )
         return None
 
     def _setup_actions(self, pot=None, **kwargs):
@@ -2694,15 +2755,28 @@ class Orbit:
         #    tz = numpy.zeros(self.size)
         #    tvz = numpy.zeros(self.size)
         # self.dim() == 1 error caught by _setupaA
-        self._aA_jr, self._aA_jp, self._aA_jz = self._aA(
-            self.R(use_physical=False, dontreshape=True),
-            self.vR(use_physical=False, dontreshape=True),
-            self.vT(use_physical=False, dontreshape=True),
-            tz,
-            tvz,
-            self.phi(use_physical=False, dontreshape=True),
-            use_physical=False,
+        # Exclude unbound orbits
+        indx, aAkwargs = self._unbound_indx_and_aAkwargs()
+        (
+            self._aA_jr,
+            self._aA_jp,
+            self._aA_jz,
+        ) = (
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
+            numpy.zeros_like(self.R(use_physical=False, dontreshape=True)) + numpy.nan,
         )
+        if numpy.sum(indx) > 0:
+            self._aA_jr[indx], self._aA_jp[indx], self._aA_jz[indx] = self._aA(
+                self.R(use_physical=False, dontreshape=True)[indx],
+                self.vR(use_physical=False, dontreshape=True)[indx],
+                self.vT(use_physical=False, dontreshape=True)[indx],
+                tz[indx],
+                tvz[indx],
+                self.phi(use_physical=False, dontreshape=True)[indx],
+                use_physical=False,
+                **aAkwargs,
+            )
         return None
 
     @shapeDecorator
diff --git a/galpy/potential/AnySphericalPotential.py b/galpy/potential/AnySphericalPotential.py
index 35d3ef90..52f3c5a2 100644
--- a/galpy/potential/AnySphericalPotential.py
+++ b/galpy/potential/AnySphericalPotential.py
@@ -83,7 +83,7 @@ class AnySphericalPotential(SphericalPotential):
             lambda r: 4.0
             * numpy.pi
             * integrate.quad(
-                lambda a: a**2 * self._rawdens(a), 0, numpy.atleast_1d(r)[0]
+                lambda a: a**2 * self._rawdens(a), 0, numpy.atleast_1d(r).flatten()[0]
             )[0]
         )
         # The potential at zero, try to figure out whether it's finite
@@ -125,7 +125,9 @@ class AnySphericalPotential(SphericalPotential):
                 - 4.0
                 * numpy.pi
                 * integrate.quad(
-                    lambda a: self._rawdens(a) * a, numpy.atleast_1d(r)[0], numpy.inf
+                    lambda a: self._rawdens(a) * a,
+                    numpy.atleast_1d(r).flatten()[0],
+                    numpy.inf,
                 )[0]
             )
 


Test patch:
diff --git a/tests/test_actionAngle.py b/tests/test_actionAngle.py
index d5e5258f..247feed8 100644
--- a/tests/test_actionAngle.py
+++ b/tests/test_actionAngle.py
@@ -4927,6 +4927,733 @@ def test_orbit_interface_actionAngleIsochroneApprox():
     return None
 
 
+def test_orbit_interface_unbound_simple_adiabatic_noc():
+    # Test that an unbound orbit in a set of orbits is handled correctly
+    from galpy.actionAngle import actionAngleAdiabatic
+    from galpy.orbit import Orbit
+    from galpy.potential import MWPotential2014
+
+    obs = Orbit(
+        [[1.05, 0.02, 1.05, 0.03, 0.0, 2.0], [1.05, 0.02, 10.05, 0.03, 0.0, 2.0]]
+    )
+    aAAnoc = actionAngleAdiabatic(pot=MWPotential2014, c=False)
+    jr, jp, jz, e, zmax, rperi, rap = (
+        obs.jr(pot=MWPotential2014, type="adiabatic", c=False),
+        obs.jp(pot=MWPotential2014, type="adiabatic", c=False),
+        obs.jz(pot=MWPotential2014, type="adiabatic", c=False),
+        obs.e(pot=MWPotential2014, type="adiabatic", analytic=True, c=False),
+        obs.zmax(pot=MWPotential2014, type="adiabatic", analytic=True, c=False),
+        obs.rperi(pot=MWPotential2014, type="adiabatic", analytic=True, c=False),
+        obs.rap(pot=MWPotential2014, type="adiabatic", analytic=True, c=False),
+    )
+    assert numpy.fabs(jr[0] - aAAnoc(obs[0])[0]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jp[0] - aAAnoc(obs[0])[1]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jz[0] - aAAnoc(obs[0])[2]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(e[0] - aAAnoc.EccZmaxRperiRap(obs[0])[0]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(zmax[0] - aAAnoc.EccZmaxRperiRap(obs[0])[1]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(rperi[0] - aAAnoc.EccZmaxRperiRap(obs[0])[2]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(rap[0] - aAAnoc.EccZmaxRperiRap(obs[0])[3]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jr[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jp[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jz[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(e[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(zmax[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(rperi[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(rap[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    return None
+
+
+def test_orbit_interface_unbound_simple_adiabatic_c():
+    # Test that an unbound orbit in a set of orbits is handled correctly
+    from galpy.actionAngle import actionAngleAdiabatic
+    from galpy.orbit import Orbit
+    from galpy.potential import MWPotential2014
+
+    obs = Orbit(
+        [[1.05, 0.02, 1.05, 0.03, 0.0, 2.0], [1.05, 0.02, 10.05, 0.03, 0.0, 2.0]]
+    )
+    aAAc = actionAngleAdiabatic(pot=MWPotential2014, c=True)
+    jr, jp, jz, e, zmax, rperi, rap = (
+        obs.jr(pot=MWPotential2014, type="adiabatic", c=True),
+        obs.jp(pot=MWPotential2014, type="adiabatic", c=True),
+        obs.jz(pot=MWPotential2014, type="adiabatic", c=True),
+        obs.e(pot=MWPotential2014, type="adiabatic", analytic=True, c=True),
+        obs.zmax(pot=MWPotential2014, type="adiabatic", analytic=True, c=True),
+        obs.rperi(pot=MWPotential2014, type="adiabatic", analytic=True, c=True),
+        obs.rap(pot=MWPotential2014, type="adiabatic", analytic=True, c=True),
+    )
+    assert numpy.fabs(jr[0] - aAAc(obs[0])[0]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jp[0] - aAAc(obs[0])[1]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jz[0] - aAAc(obs[0])[2]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(e[0] - aAAc.EccZmaxRperiRap(obs[0])[0]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(zmax[0] - aAAc.EccZmaxRperiRap(obs[0])[1]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(rperi[0] - aAAc.EccZmaxRperiRap(obs[0])[2]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(rap[0] - aAAc.EccZmaxRperiRap(obs[0])[3]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jr[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jp[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jz[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(e[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(zmax[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(rperi[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(rap[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    return None
+
+
+def test_orbit_interface_unbound_simple_staeckel_noc():
+    # Test that an unbound orbit in a set of orbits is handled correctly
+    from galpy.actionAngle import actionAngleStaeckel
+    from galpy.orbit import Orbit
+    from galpy.potential import MWPotential2014
+
+    obs = Orbit(
+        [[1.05, 0.02, 1.05, 0.03, 0.0, 2.0], [1.05, 0.02, 10.05, 0.03, 0.0, 2.0]]
+    )
+    aASnoc = actionAngleStaeckel(pot=MWPotential2014, delta=0.71, c=False)
+    jr, jp, jz, e, zmax, rperi, rap = (
+        obs.jr(pot=MWPotential2014, type="staeckel", delta=0.71, c=False),
+        obs.jp(pot=MWPotential2014, type="staeckel", delta=0.71, c=False),
+        obs.jz(pot=MWPotential2014, type="staeckel", delta=0.71, c=False),
+        obs.e(pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True, c=False),
+        obs.zmax(
+            pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True, c=False
+        ),
+        obs.rperi(
+            pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True, c=False
+        ),
+        obs.rap(
+            pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True, c=False
+        ),
+    )
+    assert numpy.fabs(jr[0] - aASnoc(obs[0])[0]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jp[0] - aASnoc(obs[0])[1]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jz[0] - aASnoc(obs[0])[2]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(e[0] - aASnoc.EccZmaxRperiRap(obs[0])[0]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(zmax[0] - aASnoc.EccZmaxRperiRap(obs[0])[1]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(rperi[0] - aASnoc.EccZmaxRperiRap(obs[0])[2]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(rap[0] - aASnoc.EccZmaxRperiRap(obs[0])[3]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jr[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jp[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jz[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(e[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(zmax[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(rperi[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(rap[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    return None
+
+
+def test_orbit_interface_unbound_simple_staeckel_c():
+    # Test that an unbound orbit in a set of orbits is handled correctly
+    from galpy.actionAngle import actionAngleStaeckel
+    from galpy.orbit import Orbit
+    from galpy.potential import MWPotential2014
+
+    obs = Orbit(
+        [[1.05, 0.02, 1.05, 0.03, 0.0, 2.0], [1.05, 0.02, 10.05, 0.03, 0.0, 2.0]]
+    )
+    aASc = actionAngleStaeckel(pot=MWPotential2014, delta=0.71, c=True)
+    jr, jp, jz, omr, omp, omz, wr, wp, wz, e, zmax, rperi, rap = (
+        obs.jr(pot=MWPotential2014, type="staeckel", delta=0.71, c=True),
+        obs.jp(pot=MWPotential2014, type="staeckel", delta=0.71, c=True),
+        obs.jz(pot=MWPotential2014, type="staeckel", delta=0.71, c=True),
+        obs.Or(pot=MWPotential2014, type="staeckel", delta=0.71, c=True),
+        obs.Op(pot=MWPotential2014, type="staeckel", delta=0.71, c=True),
+        obs.Oz(pot=MWPotential2014, type="staeckel", delta=0.71, c=True),
+        obs.wr(pot=MWPotential2014, type="staeckel", delta=0.71, c=True),
+        obs.wp(pot=MWPotential2014, type="staeckel", delta=0.71, c=True),
+        obs.wz(pot=MWPotential2014, type="staeckel", delta=0.71, c=True),
+        obs.e(pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True, c=True),
+        obs.zmax(
+            pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True, c=True
+        ),
+        obs.rperi(
+            pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True, c=True
+        ),
+        obs.rap(
+            pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True, c=True
+        ),
+    )
+    assert numpy.fabs(jr[0] - aASc(obs[0])[0]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jp[0] - aASc(obs[0])[1]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jz[0] - aASc(obs[0])[2]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(omr[0] - aASc.actionsFreqs(obs[0])[3]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(omp[0] - aASc.actionsFreqs(obs[0])[4]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(omz[0] - aASc.actionsFreqs(obs[0])[5]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(wr[0] - aASc.actionsFreqsAngles(obs[0])[6]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(wp[0] - aASc.actionsFreqsAngles(obs[0])[7]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(wz[0] - aASc.actionsFreqsAngles(obs[0])[8]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(e[0] - aASc.EccZmaxRperiRap(obs[0])[0]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(zmax[0] - aASc.EccZmaxRperiRap(obs[0])[1]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(rperi[0] - aASc.EccZmaxRperiRap(obs[0])[2]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(rap[0] - aASc.EccZmaxRperiRap(obs[0])[3]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jr[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jp[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jz[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(omr[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(omp[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(omz[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(wr[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(wp[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(wz[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    return None
+
+
+def test_orbit_interface_unbound_simple_2d_adiabatic():
+    # Test that an unbound orbit in a set of orbits is handled correctly
+    from galpy.actionAngle import actionAngleSpherical
+    from galpy.orbit import Orbit
+    from galpy.potential import MWPotential2014
+
+    obs = Orbit([[1.05, 0.02, 1.05, 2.0], [1.05, 0.02, 10.05, 2.0]])
+    # in 2D, adiabatic and Staeckel are the same and the same as spherical
+    aAS = actionAngleSpherical(pot=MWPotential2014)
+    jr, jp, jz = (
+        obs.jr(pot=MWPotential2014, type="adiabatic"),
+        obs.jp(pot=MWPotential2014, type="adiabatic"),
+        obs.jz(pot=MWPotential2014, type="adiabatic"),
+    )
+    assert numpy.fabs(jr[0] - aAS(obs[0])[0]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jp[0] - aAS(obs[0])[1]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jz[0] - aAS(obs[0])[2]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jr[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jp[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jz[1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    return None
+
+
+def test_orbit_interface_unbound_simple_2d_staeckel():
+    # Test that an unbound orbit in a set of orbits is handled correctly
+    from galpy.actionAngle import actionAngleSpherical
+    from galpy.orbit import Orbit
+    from galpy.potential import MWPotential2014
+
+    obs = Orbit([[1.05, 0.02, 1.05, 2.0], [1.05, 0.02, 10.05, 2.0]])
+    # in 2D, adiabatic and Staeckel are the same and the same as spherical
+    aAS = actionAngleSpherical(pot=MWPotential2014)
+    jr, jp, jz = (
+        obs.jr(pot=MWPotential2014, type="staeckel", delta=0.71),
+        obs.jp(pot=MWPotential2014, type="staeckel", delta=0.71),
+        obs.jz(pot=MWPotential2014, type="staeckel", delta=0.71),
+    )
+    assert numpy.fabs(jr[0] - aAS(obs[0])[0]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jp[0] - aAS(obs[0])[1]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.fabs(jz[0] - aAS(obs[0])[2]) < 10.0**-10.0, (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jr[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jp[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jz[1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    return None
+
+
+def test_orbit_interface_unbound_complexshape_adiabatic():
+    # Test that an unbound orbit in a set of orbits is handled correctly
+    from galpy.actionAngle import actionAngleAdiabatic
+    from galpy.orbit import Orbit
+    from galpy.potential import MWPotential2014
+
+    obs = Orbit(
+        numpy.array(
+            [
+                [
+                    [1.05, 0.02, 1.05, 0.03, 0.0, 2.0],
+                    [1.05, 0.02, 10.05, 0.03, 0.0, 2.0],
+                ],
+                [
+                    [1.05, 0.02, 1.05, 0.03, 0.0, 2.0],
+                    [1.05, 0.02, 10.05, 0.03, 0.0, 2.0],
+                ],
+                [
+                    [1.05, 0.02, 1.05, 0.03, 0.0, 2.0],
+                    [1.05, 0.02, 10.05, 0.03, 0.0, 2.0],
+                ],
+            ]
+        )
+    )
+    aAA = actionAngleAdiabatic(pot=MWPotential2014)
+    jr, jp, jz, e, zmax, rperi, rap = (
+        obs.jr(pot=MWPotential2014, type="adiabatic"),
+        obs.jp(pot=MWPotential2014, type="adiabatic"),
+        obs.jz(pot=MWPotential2014, type="adiabatic"),
+        obs.e(pot=MWPotential2014, type="adiabatic", analytic=True),
+        obs.zmax(pot=MWPotential2014, type="adiabatic", analytic=True),
+        obs.rperi(pot=MWPotential2014, type="adiabatic", analytic=True),
+        obs.rap(pot=MWPotential2014, type="adiabatic", analytic=True),
+    )
+    assert numpy.all(numpy.fabs(jr[:, 0] - aAA(obs[:, 0])[0]) < 10.0**-10.0), (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.fabs(jp[:, 0] - aAA(obs[:, 0])[1]) < 10.0**-10.0), (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.fabs(jz[:, 0] - aAA(obs[:, 0])[2]) < 10.0**-10.0), (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(e[:, 0] - aAA.EccZmaxRperiRap(obs[:, 0])[0]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(zmax[:, 0] - aAA.EccZmaxRperiRap(obs[:, 0])[1]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(rperi[:, 0] - aAA.EccZmaxRperiRap(obs[:, 0])[2]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(rap[:, 0] - aAA.EccZmaxRperiRap(obs[:, 0])[3]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleAdiabatic does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jr[:, 1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jp[:, 1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jz[:, 1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(e[:, 1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(zmax[:, 1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(rperi[:, 1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(rap[:, 1])), (
+        "Orbit interface for actionAngleAdiabatic does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    return None
+
+
+def test_orbit_interface_unbound_complexshape_staeckel():
+    # Test that an unbound orbit in a set of orbits is handled correctly
+    from galpy.actionAngle import actionAngleStaeckel
+    from galpy.orbit import Orbit
+    from galpy.potential import MWPotential2014
+
+    obs = Orbit(
+        numpy.array(
+            [
+                [
+                    [1.05, 0.02, 1.05, 0.03, 0.0, 2.0],
+                    [1.05, 0.02, 10.05, 0.03, 0.0, 2.0],
+                ],
+                [
+                    [1.05, 0.02, 1.05, 0.03, 0.0, 2.0],
+                    [1.05, 0.02, 10.05, 0.03, 0.0, 2.0],
+                ],
+                [
+                    [1.05, 0.02, 1.05, 0.03, 0.0, 2.0],
+                    [1.05, 0.02, 10.05, 0.03, 0.0, 2.0],
+                ],
+            ]
+        )
+    )
+    aAS = actionAngleStaeckel(pot=MWPotential2014, delta=0.71)
+    jr, jp, jz, omr, omp, omz, wr, wp, wz, e, zmax, rperi, rap = (
+        obs.jr(pot=MWPotential2014, type="staeckel", delta=0.71),
+        obs.jp(pot=MWPotential2014, type="staeckel", delta=0.71),
+        obs.jz(pot=MWPotential2014, type="staeckel", delta=0.71),
+        obs.Or(pot=MWPotential2014, type="staeckel", delta=0.71),
+        obs.Op(pot=MWPotential2014, type="staeckel", delta=0.71),
+        obs.Oz(pot=MWPotential2014, type="staeckel", delta=0.71),
+        obs.wr(pot=MWPotential2014, type="staeckel", delta=0.71),
+        obs.wp(pot=MWPotential2014, type="staeckel", delta=0.71),
+        obs.wz(pot=MWPotential2014, type="staeckel", delta=0.71),
+        obs.e(pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True),
+        obs.zmax(pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True),
+        obs.rperi(pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True),
+        obs.rap(pot=MWPotential2014, type="staeckel", delta=0.71, analytic=True),
+    )
+    assert numpy.all(numpy.fabs(jr[:, 0] - aAS(obs[:, 0])[0]) < 10.0**-10.0), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.fabs(jp[:, 0] - aAS(obs[:, 0])[1]) < 10.0**-10.0), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.fabs(jz[:, 0] - aAS(obs[:, 0])[2]) < 10.0**-10.0), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(omr[:, 0] - aAS.actionsFreqs(obs[:, 0])[3]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(omp[:, 0] - aAS.actionsFreqs(obs[:, 0])[4]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(omz[:, 0] - aAS.actionsFreqs(obs[:, 0])[5]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(wr[:, 0] - aAS.actionsFreqsAngles(obs[:, 0])[6]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(wp[:, 0] - aAS.actionsFreqsAngles(obs[:, 0])[7]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(wz[:, 0] - aAS.actionsFreqsAngles(obs[:, 0])[8]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(e[:, 0] - aAS.EccZmaxRperiRap(obs[:, 0])[0]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(zmax[:, 0] - aAS.EccZmaxRperiRap(obs[:, 0])[1]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(rperi[:, 0] - aAS.EccZmaxRperiRap(obs[:, 0])[2]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(rap[:, 0] - aAS.EccZmaxRperiRap(obs[:, 0])[3]) < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jr[:, 1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jp[:, 1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jz[:, 1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    return None
+
+
+def test_orbit_interface_unbound_staeckeldelta_handling():
+    # Test that the automagically determined delta is handled correctly when there are unbound orbits
+    # Use a complex shape
+    from galpy.actionAngle import actionAngleStaeckel
+    from galpy.orbit import Orbit
+    from galpy.potential import MWPotential2014
+
+    obs = Orbit(
+        numpy.array(
+            [
+                [
+                    [1.15, 0.02, 1.15, 0.03, 0.0, 2.0],
+                    [1.05, 0.02, 10.05, 0.03, 0.0, 2.0],
+                ],
+                [
+                    [1.02, 0.02, 0.95, 0.03, 0.0, 2.0],
+                    [1.05, 0.02, 10.05, 0.03, 0.0, 2.0],
+                ],
+                [
+                    [0.97, 0.02, 1.25, 0.03, 0.0, 2.0],
+                    [1.05, 0.02, 10.05, 0.03, 0.0, 2.0],
+                ],
+            ]
+        )
+    )
+    # Compute the actions with the automagically determined delta using the orbit interface
+    jr, jp, jz, omr, omp, omz, wr, wp, wz, e, zmax, rperi, rap = (
+        obs.jr(pot=MWPotential2014, type="staeckel"),
+        obs.jp(pot=MWPotential2014, type="staeckel"),
+        obs.jz(pot=MWPotential2014, type="staeckel"),
+        obs.Or(pot=MWPotential2014, type="staeckel"),
+        obs.Op(pot=MWPotential2014, type="staeckel"),
+        obs.Oz(pot=MWPotential2014, type="staeckel"),
+        obs.wr(pot=MWPotential2014, type="staeckel"),
+        obs.wp(pot=MWPotential2014, type="staeckel"),
+        obs.wz(pot=MWPotential2014, type="staeckel"),
+        obs.e(pot=MWPotential2014, type="staeckel", analytic=True),
+        obs.zmax(pot=MWPotential2014, type="staeckel", analytic=True),
+        obs.rperi(pot=MWPotential2014, type="staeckel", analytic=True),
+        obs.rap(pot=MWPotential2014, type="staeckel", analytic=True),
+    )
+    # Now do the same with the actionAngle interface
+    aAS = actionAngleStaeckel(pot=MWPotential2014, delta=0.71)  # just a dummy delta
+    bound_indx = numpy.array([True, False, True, False, True, False])
+    assert numpy.all(
+        numpy.fabs(jr[:, 0] - aAS(obs[:, 0], delta=obs._aA._delta[bound_indx])[0])
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(jp[:, 0] - aAS(obs[:, 0], delta=obs._aA._delta[bound_indx])[1])
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(jz[:, 0] - aAS(obs[:, 0], delta=obs._aA._delta[bound_indx])[2])
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(
+            omr[:, 0] - aAS.actionsFreqs(obs[:, 0], delta=obs._aA._delta[bound_indx])[3]
+        )
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(
+            omp[:, 0] - aAS.actionsFreqs(obs[:, 0], delta=obs._aA._delta[bound_indx])[4]
+        )
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(
+            omz[:, 0] - aAS.actionsFreqs(obs[:, 0], delta=obs._aA._delta[bound_indx])[5]
+        )
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(
+            wr[:, 0]
+            - aAS.actionsFreqsAngles(obs[:, 0], delta=obs._aA._delta[bound_indx])[6]
+        )
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(
+            wp[:, 0]
+            - aAS.actionsFreqsAngles(obs[:, 0], delta=obs._aA._delta[bound_indx])[7]
+        )
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(
+            wz[:, 0]
+            - aAS.actionsFreqsAngles(obs[:, 0], delta=obs._aA._delta[bound_indx])[8]
+        )
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(
+            e[:, 0]
+            - aAS.EccZmaxRperiRap(obs[:, 0], delta=obs._aA._delta[bound_indx])[0]
+        )
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(
+            zmax[:, 0]
+            - aAS.EccZmaxRperiRap(obs[:, 0], delta=obs._aA._delta[bound_indx])[1]
+        )
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(
+            rperi[:, 0]
+            - aAS.EccZmaxRperiRap(obs[:, 0], delta=obs._aA._delta[bound_indx])[2]
+        )
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(
+        numpy.fabs(
+            rap[:, 0]
+            - aAS.EccZmaxRperiRap(obs[:, 0], delta=obs._aA._delta[bound_indx])[3]
+        )
+        < 10.0**-10.0
+    ), (
+        "Orbit interface for actionAngleStaeckel does not return the same as actionAngle interface for bound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jr[:, 1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jp[:, 1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    assert numpy.all(numpy.isnan(jz[:, 1])), (
+        "Orbit interface for actionAngleStaeckel does not return NaN for unbound orbit in a collection with an unbound orbit"
+    )
+    return None
+
+
 # Test physical output for actionAngleStaeckel
 def test_physical_staeckel():
     from galpy.actionAngle import actionAngleStaeckel