diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index 4ad6e801f..2c19efcce 100644 --- a/galpy/orbit/Orbits.py +++ b/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 diff --git a/tests/test_orbit.py b/tests/test_orbit.py index 3d2674234..288ba3435 100644 --- a/tests/test_orbit.py +++ b/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