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).

Hints from issue comments:
@RaynaRampalli: Hi, I'm running into this issue when calculating actions for 3+ million stars using galpy and am looping thru each star one by one. Do you have any suggested simple workarounds for suppressing the error and just returning nans or something for now? I've tried a couple things and no luck. Thank you! @jobovy: > Hi, I'm running into this issue when calculating actions for 3+ million stars using galpy and am looping thru each star one by one. Do you have any suggested simple workarounds for suppressing the error and just returning nans or something for now? I've tried a couple things and no luck. Thank you!

Thanks for reaching out, sounds like we might actually finally have to fix this issue! Quick question: Can you determine the stars for which the action-angle calculation fails just by determining which stars are unbound, that is, have E > E(infinity)? For example, assuming the ``orb`` object is an ``Orbit`` instance with multiple objects, some of which have failing action-angles, do
```python
import numpy
from galpy.potential import MWPotential2014, evaluatePotentials
indx= orb.E(pot=MWPotential2014,use_physical=False) < evaluatePotentials(MWPotential2014,numpy.inf,0.,use_physical=False)
```
Does ``indx`` than have the index for all objects that have successful action-angle coordiantes (and ``True^indx`` the index for all objects for which they fail?). If that method works, then it should be easy for you to do and I can then also use it to create a fix for this issue. Thanks!@jobovy: Hi @RaynaRampalli! Just wanted to check whether you've had a chance to try the fix that I suggested above? I'm getting ready to release a new version of galpy soon and it would be nice to include a proper fix to this issue then. Thanks!@RaynaRampalli: Hi! Sorry for the delay; thank you for responding.  I have been traveling so have not had a chance yet. I will try this out and will have an answer for you in about a week. Thanks! Best,Rayna On Jun 26, 2024, at 10:09 AM, Jo Bovy ***@***.***> wrote:﻿
Hi @RaynaRampalli! Just wanted to check whether you've had a chance to try the fix that I suggested above? I'm getting ready to release a new version of galpy soon and it would be nice to include a proper fix to this issue then. Thanks!

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>@jobovy: This fell off my radar until I came across this issue again just now, but I was wondering whether you ever had a chance to check the fix that I suggested @RaynaRampalli?@RaynaRampalli: So sorry, I got really busy with travel and application season. I just
tried out your solution, and it does work! Will incorporate it into my
action calculations going forward, thank you so much.

On Thu, Nov 14, 2024 at 9:43 AM Jo Bovy ***@***.***> wrote:

> This fell off my radar until I came across this issue again just now, but
> I was wondering whether you ever had a chance to check the fix that I
> suggested @RaynaRampalli <https://github.com/RaynaRampalli>?
>
> —
> Reply to this email directly, view it on GitHub
> <https://github.com/jobovy/galpy/issues/383#issuecomment-2477044052>, or
> unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AEYD372VMJUKM3HGB5WYQHT2ATOM3AVCNFSM6AAAAABJL6MGXOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDINZXGA2DIMBVGI>
> .
> You are receiving this because you were mentioned.Message ID:
> ***@***.***>
>
@jobovy: > So sorry, I got really busy with travel and application season. I just tried out your solution, and it does work! Will incorporate it into my action calculations going forward, thank you so much.

Great, thanks for the feedback!

Solution 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