Raise error when integration times (or angles for SOS) are not evenly…

… spaced for integration methods that require this; fixes #700; test extensively

galpy/orbit/Orbits.py

@@ -1342,6 +1342,16 @@ def _check_method_dissipative_compatible(method, pot):
             )
         return method
 
+    @staticmethod
+    def _check_array_evenlyspaced(t, method):  # for ODE integrators
+        if method in ["odeint", "dop853", "dop853_c"]:
+            return True
+        else:
+            diffs = numpy.diff(t, axis=-1)
+            return numpy.all(
+                numpy.isclose(diffs, numpy.expand_dims(diffs[..., 0], axis=-1))
+            )
+
     def integrate(
         self,
         t,
@@ -1358,7 +1368,7 @@ def integrate(
         Parameters
         ----------
         t : list, numpy.ndarray or Quantity
-            List of equispaced times at which to compute the orbit. The initial condition is t[0].
+            List of equispaced times at which to compute the orbit. The initial condition is t[0]. (note that for method='odeint', method='dop853', and method='dop853_c', the time array can be non-equispaced).
         pot : Potential, DissipativeForce or list of such instances
             Gravitational field to integrate the orbit in.
         method : str, optional
@@ -1406,8 +1416,15 @@ def integrate(
             t = conversion.parse_time(t, ro=self._ro, vo=self._vo)
         else:
             self._integrate_t_asQuantity = False
+        # Check that t is evenly spaced if not using odeint
+        if not self._check_array_evenlyspaced(t, method):
+            raise ValueError(
+                f"Input time array must be equally spaced for method {method}, use method='dop853_c', method='dop853', or method='odeint' instead for non-equispaced time arrays"
+            )
+
         if _APY_LOADED and not dt is None and isinstance(dt, units.Quantity):
             dt = conversion.parse_time(dt, ro=self._ro, vo=self._vo)
+
         from ..potential import MWPotential
 
         if pot == MWPotential:
@@ -1575,7 +1592,7 @@ def integrate_SOS(
         Parameters
         ----------
         psi : list, numpy.ndarray or Quantity
-            Equispaced list of increment angles over which to integrate [increments wrt initial angle].
+            Equispaced list of increment angles over which to integrate [increments wrt initial angle].  (note that for method='odeint', method='dop853', and method='dop853_c', the psi array can be non-equispaced).
         pot : Potential, DissipativeForce or list of such instances
             Gravitational field to integrate the orbit in.
         surface : str, optional
@@ -1619,6 +1636,11 @@ def integrate_SOS(
         # Parse psi
         if _APY_LOADED and isinstance(psi, units.Quantity):
             psi = conversion.parse_angle(psi)
+        # Check that psi is evenly spaced
+        if not self._check_array_evenlyspaced(psi, method):
+            raise ValueError(
+                f"Input psi array must be equally spaced for method {method}, use method='dop853_c', method='dop853', or method='odeint' instead for non-equispaced psi arrays"
+            )
         if _APY_LOADED and isinstance(t0, units.Quantity):
             t0 = conversion.parse_time(t0, ro=self._ro, vo=self._vo)
         self._integrate_t_asQuantity = False
@@ -1727,7 +1749,7 @@ def integrate_dxdv(
         dxdv : numpy.ndarray
             Initial conditions for the orbit in cylindrical or rectangular coordinates. The shape of the array should be (\*input_shape, 4).
         t : list, numpy.ndarray or Quantity
-            List of equispaced times at which to compute the orbit. The initial condition is t[0].
+            List of equispaced times at which to compute the orbit. The initial condition is t[0].  (note that for method='odeint', method='dop853', and method='dop853_c', the time array can be non-equispaced).
         pot : Potential, DissipativeForce or list of such instances
             Gravitational field to integrate the orbit in.
         method : str, optional
@@ -1780,6 +1802,11 @@ def integrate_dxdv(
             t = conversion.parse_time(t, ro=self._ro, vo=self._vo)
         else:
             self._integrate_t_asQuantity = False
+        # Check that t is evenly spaced
+        if not self._check_array_evenlyspaced(t, method):
+            raise ValueError(
+                f"Input time array must be equally spaced for method {method}, use method='dop853_c', method='dop853', or method='odeint' instead for non-equispaced time arrays"
+            )
         if not dt is None:
             dt = conversion.parse_time(dt, ro=self._ro, vo=self._vo)
         # Parse dxdv

tests/test_orbit.py

@@ -8075,6 +8075,91 @@ def test_orbinterp_reset_integratedxdv():
     return None
 
 
+# Test that an error is raised when integration time array is not equally spaced (see #700)
+def test_integrate_notevenlyspaced_issue700():
+    from galpy.orbit import Orbit
+    from galpy.potential import MWPotential2014
+
+    times = numpy.concatenate(
+        [numpy.linspace(0, 10, 21), numpy.linspace(12.0, 50.0, 20)]
+    )
+    orb = Orbit()
+    # Test that the correct error is raised when the time array is not equally spaced
+    with pytest.raises(ValueError) as excinfo:
+        orb.integrate(times, MWPotential2014, method="symplec6_c")
+    assert (
+        str(excinfo.value)
+        == "Input time array must be equally spaced for method symplec6_c, use method='dop853_c', method='dop853', or method='odeint' instead for non-equispaced time arrays"
+    ), "Input time array must be equally spaced error not raised"
+    # Also test backwards integration
+    with pytest.raises(ValueError) as excinfo:
+        orb.integrate(-times, MWPotential2014, method="symplec6_c")
+    assert (
+        str(excinfo.value)
+        == "Input time array must be equally spaced for method symplec6_c, use method='dop853_c', method='dop853', or method='odeint' instead for non-equispaced time arrays"
+    ), "Input time array must be equally spaced error not raised"
+    # Also test integrate_dxdv
+    with pytest.raises(ValueError) as excinfo:
+        orb.toPlanar().integrate_dxdv(None, times, MWPotential2014, method="dopr54_c")
+    assert (
+        str(excinfo.value)
+        == "Input time array must be equally spaced for method dopr54_c, use method='dop853_c', method='dop853', or method='odeint' instead for non-equispaced time arrays"
+    ), "Input time array must be equally spaced error not raised"
+    # Also test integrateSOS, just use times for psi...
+    with pytest.raises(ValueError) as excinfo:
+        orb.integrate_SOS(times, MWPotential2014, method="dopr54_c")
+    assert (
+        str(excinfo.value)
+        == "Input psi array must be equally spaced for method dopr54_c, use method='dop853_c', method='dop853', or method='odeint' instead for non-equispaced psi arrays"
+    ), "Input time array must be equally spaced error not raised"
+    return None
+
+
+# Test that integrators that should be fine with unevenly-spaced times are fine with it
+def test_integrate_notevenlyspaced_ok():
+    from galpy.orbit import Orbit
+    from galpy.potential import MWPotential2014
+
+    time_1 = numpy.linspace(0, 50.0, 1001)
+    time_2 = numpy.concatenate(
+        [numpy.linspace(0, 20.0, 1001), numpy.linspace(20.02, 50.0, 1001)]
+    )
+    for integrator in ["odeint", "dop853", "dop853_c"]:
+        o_1 = Orbit()
+        o_1.integrate(time_1, MWPotential2014, method=integrator)
+        o_2 = Orbit()
+        o_2.integrate(time_2, MWPotential2014, method=integrator)
+        assert numpy.all(numpy.fabs(o_1.R(time_1) - o_2.R(time_1)) < 10.0**-5.0), (
+            "Integration with unevenly-spaced times does not work"
+        )
+        assert numpy.all(numpy.fabs(o_1.vR(time_1) - o_2.vR(time_1)) < 10.0**-4.0), (
+            "Integration with unevenly-spaced times does not work"
+        )
+        assert numpy.all(numpy.fabs(o_1.vT(time_1) - o_2.vT(time_1)) < 10.0**-4.0), (
+            "Integration with unevenly-spaced times does not work"
+        )
+    for integrator in [
+        "leapfrog",
+        "leapfrog_c",
+        "symplec4_c",
+        "symplec6_c",
+        "rk4_c",
+        "rk6_c",
+        "dopr54_c",
+        "ias15_c",
+    ]:
+        o_1 = Orbit()
+        o_1.integrate(time_1, MWPotential2014, method=integrator)
+        o_2 = Orbit()
+        with pytest.raises(ValueError) as excinfo:
+            o_2.integrate(time_2, MWPotential2014, method=integrator)
+        assert (
+            str(excinfo.value)
+            == f"Input time array must be equally spaced for method {integrator}, use method='dop853_c', method='dop853', or method='odeint' instead for non-equispaced time arrays"
+        ), f"Input time array must be equally spaced for method{integrator}"
+    return None
+
+
 def test_linear_plotting():
     from galpy.orbit import Orbit
     from galpy.potential.verticalPotential import RZToverticalPotential