The reason for this behavior is that solar radiation pressure generates discontinuous acceleration at eclipses umbra entry/exit and penumbra entry/exit. These discontinuities are handled by automatically setting up internal event detectors that trigger a
RESET_DERIVATIVES action when the events occur.
This reset is forwarded from the numerical propagator to the underlying ODE integrator, which enforces a new step to start at this point. This is straightforward for single step integrators like
DormandPrince853Integrator that do not preserve memory of previous steps and for which starting a new step doesn’t involve any fancy computation. This is not as easy for multistep integrators like
AdamsMoultonIntegrator as they need to wipe away information gathered on previous steps and rebuild a new model (this information is stored in a Nordsieck vector). For rebuilding this model, multistep integrators use an underlying starter integrator
which by default is a
DormandPrince853Integrator. This restarting procedure requires the underlying integrator to be able to run for (n + 3) / 2 steps where n is the number for steps of the Adams method used. In your case, n = 8 so the starter integrator must compute 5 steps.
What happens in your case is that due to bad luck, your propagation range goes from 00:00:00.000 to 06:00:00.000 (i.e. 0.25 day) and a penumbra exit event occurs at 05:59:56.973, i.e. 3.027 seconds before propagation end. There are several other umbra/penumbra entry/exit events before, but this last one is the culprit. After this event the Adams-Bashforth integrator experience a restart, but the underlying Dormand-Prince integrator cannot perform 5 steps in the short 3.027 seconds range that remains before integration end. This triggers the exception.
Note also that your setup for step size handling could be improved. You use constants 1e-4 for absolute and relative tolerances on the integrator and do not specify the orbit type to use for integration, which leads to
OrbitType.EQUINOCTIAL to be selected by default. We recommend to set up the orbit type explicitly and to use it for both configuring the propagator and the tolerances in adaptive stepsize. You can do it using this code snippet in your case:
OrbitType type = OrbitType.EQUINOCTIAL;
double tol = NumericalPropagator.tolerances(0.1, orbit, type);
AdamsBashforthIntegrator integrator = new AdamsBashforthIntegrator(8, 1e-15, timeSpan, tol, tol);
NumericalPropagator propagator = new NumericalPropagator(integrator);
Finally, when running your test case, I found there was indeed a bug in the way Hipparchus handles the reset and hence the restarts. The
resetOccurred flag is not reset to false after the reset has been handled and spurious restarts occur between events, greatly increasing computation costs with multistep integrators. Could you open a bug report on https://github.com/Hipparchus-Math/hipparchus/issues for this? I will fix this bug soon but as a workaround, you can just add
resetOccurred = false; before the
boolean doneWithStep = false; statement at line 314 in Hipparchus
So as a conclusion, there are three things you could do:
- extend your propagation time so it does not end just after a discontinuity-inducing event
- configure tolerances in step size management using
and an explicit
AbstractIntegrator to avoid too many restarts