DerivativeStructures with AdaptiveStepsizeIntegrator for reentry

Hi everyone!

I am currently trying to implement Orekit’s Taylor Algebra functionality into a reentry propagating tool. This tool, in it’s current (non-Taylor Algebra) version, uses a DormandPrince853Integrator, which is generally a good choice for a variable stepsize method. I have successfully used this integrator (or rather its “field” version) for Taylor Algebra-based propagations in the past, but never for reentry.

The problem comes when propagating until an AltitudeDetector (Alt = 0.) event, where the derivatives in the DerivativeStructure variables become null (real order-zero value is still kept). I have made several observations about this:

  • If the AltitudeDetector is set to a considerably higher altitude (e.g. 10 km), derivatives are still there, suggesting that numerics are not stable when the orbit becomes too vertical.
  • Every other variable stepsize method integrator available in Hipparchus suffers the same problem
  • If, however, a fixed stepsize method is chosen instead (RungeKutta in any of its variations), propagation can take place until Alt = 0 and derivatives are kept. This obviously lacks all the benefits of a variable stepsize integrator.

For context, I use a FieldNumericalPropagator with a Holmes-Featherstone attraction model, Sun and Moon as 3rd bodies, solar radiation pressure and an implementation of IsotropicDrag that let’s me input a DerivativeStructure as drag coefficient. Minimum and maximum stepsize for the DormandPrince integrator are respectively [1e-20 s, 3e02 s] and position tolerance is set to 1e-8 m.

Is there any chance my config is wrong here? If not, I would really appreciate suggestions to avoid or improve the fixed-step method usage.

Thank you and cheers,


I don’t know for sure, but my wild guess would be that the behaviour shown by variable stepsize algorithms is normal and behaviour shown by fixed step size is just good luck.

From my experience, computing reentry trajectories often needs special handling for two force models: gravity field and atmospheric drag. Spherical harmonics expansion are singular when crossing the Brillouin sphere. Indeed, I would even have guessed you would trigger an exception as soon as the distance to Earth center comes below equatorial radius, even if you are still above the ellipsoid if reentry is not at equator. Atmospheric drag relies on models that often have a lower altitude limitation in the range of 90km or 120km. The NRLMSISE2000 model can go down to ground as can the simple exponential model, though.

As you put the position tolerance to a really small (I would say unrealistic) value of 1.0e-8m, there are chances that when approaching the singularity, the step sizes decreases a lot and you end up evaluating points too close to the singularity and see the derivatives go weird. On the other hand, a fixed step size integrator would most probably not evaluate anything exactly at singularity, but only points slightly away from the singularity, hence not noticing it.

As Luc said. I found exactly the same thing happening. A custom ForceModel below the reentry interface altitude could help.