[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Orekit Developers] status on second order derivatives

Hi Evan,
Evan Ward <evan.ward@nrl.navy.mil> a écrit :

Hi Luc, Paul,

If I understand correctly this issue only affects propagators that do
not use Cartesian representation of the state. When a non-Cartesian
representation is used the propagator's position is combined with the
Keplerian velocity and the resulting combination is not locally tangent
(osculating) to the propagator's orbit. Am I understanding the issue

Yes, this is exactly what happens.

Could this issue be solved by using the propagator's complete state to
determine the instantaneous P/V and use that to build a CartesianOrbit?

Not exactly. The numericalPropagator once only integrated CartesianOrbit and it was considered a limitation. Since a few years, we have added the possibility to directly integrate the orbit type specifed by user. Moving back to CartesianOrbit only would be a regression IMHO.

However what you suggest is also almost what I have in mind since we really use the propagator complete state. In fact, for the numerical propagator the conversion between acceleration and orbit is not done at orbit level after integration, but rather at force level as we use the Jacobian of the (P,V) vs orbital parameters conversion to directly get parameters derivatives. When you look at it, the Gauss equations are really the last columns of the Jacobian. This is done with the TimeDerivatives interface, which for NumericalPropagator is implemented using the Jacobian.

So as far as N?umericalPropagator is involved, I think consistency is automatically preserved, since we do integrated something that is computed from the Jacobian.

The DSST propagator on the other hand propagated directly in equinoctial parameters and uses another dedicated implementation of the TimeDerivatives interface.

Then we wouldn't have to modify the existing set of Orbit classes, and
the user would see the correct osculating P/V. (This might be equivalent
to your second approach.)

As far as where to put the code, it seems like the conversion code would
be specific to the internal state representation used by the propagator,
so I think it makes sense to keep the code as private to the propagator.
Though if you think there would be other uses for the conversion, I
think a public static factory method would work well.

Yes, the conversion code is propagator dependent. The first implementation I played with is Eckstein-Hechler propagator, and it definitely is Eckstein-Hechler specific (it's a simple derivation of the original equations, so each time we did compute something like e = a * c + b * d, now we also have another statement to compute eDot = aDot * c + a * dCot + bDot * d + b * dDot). For numerical propagator we already have the derivatives since we start from the derivatives and integrate them, so its even simpler. For DSST, this will be a mix as the mean elements are integrated and the short periodics terms are computed from Fourier coefficients which are straightforward to differentiate. For ephemeris-based propagator, we will need to compute the derivatives of the underlying polynomials, which is also straightforward.

Thanks Luc for finding this issue and doing the analysis. I can see how
this would be an issue when computing the Doppler as well as time shifting.

your welcome

best regards,

Best Regards,

On 10/29/2014 06:57 AM, MAISONOBE Luc wrote:
Hi Paul,

paulcefo <paulcefo@buffalo.edu> a écrit :


Do I correctly understand that your concern is that Keplerian
transformations do apply outside the osculating space?

The problem I had was that we did use Keplerian-only expression
to set up local Taylor expansions around the current point (a few
seconds away). This was slightly wrong when all the parameters were
time-dependent and not only the anomaly was time-dependent. Of course,
the error increasing with the time offset with respect to the central
date at which the Taylor expansion is built. The fix was simply to
not forget the derivatives of these other parameters.

This Taylor expansion feature is a built-in feature available in all
Orekit orbits, it typically allows to do computation in the vicinity of
an already computed point without needed to trigger a complete
It can even be used for some computation inside the run of a propagator,
as for example when the higher level propagator takes care of the long
term propagation and at each step we need some additional points
surrounding the current step to compute attitude evolution in some
specific modes.

My concern was how to implement this fix in our current architecture,
and more precisely were to put the code: in an existing class or in
a dedicated class which would be used only by propagators.

best regards,


Dr. Paul J. Cefola
Consultant in Aerospace Systems, Spaceflight Mechanics, & Astrodynamics
Adjunct Faculty, Dept. of Mechanical and Aerospace Engineering,
University at Buffalo (SUNY)

4 Moonstone Way
Vineyard Haven, MA 02568

508-696-1884 (phone on Martha's Vineyard)
978-201-1393 (cell)


On 10/29/2014 6:02 am, MAISONOBE Luc wrote:

As some of you may be aware, I have been working for a few months on
second order derivatives in the git branch
position-velocity-acceleration. This work is still ongoing but I hope
to finish it for 7.0 and merge the branch back to master soon. For
now, there are still failing tests so I can't do it.

This change should allow us to reach several goals :

- improved accuracy in shiftedBy methods
- improved accuracy in interpolators (with user-defined
  choices to use or not first and second derivatives
  from the sample)
- improved accuracy in attitude
- removal of ugly hidden finite differences in some classes
  (most notably attitude modes) with hard-coded steps
- hopefully faster Earth transforms, by replacing Hermite
  interpolation with single point extrapolation
- availability of non-Keplerian acceleration everywhere
- availability of angular acceleration in attitude and frames
- proper composition of dynamics in frames
- possibility to propagate orbits in non-inertial frames
- possibility to propagate orbits without a central body
  (interplanetary missions, Lagrange point missions, ...)

There is one point that bothers me right now. As I removed some of the
ugly finite differences, some non-regression tests started to fail. I
finally found the raw cause of these failures and was surprised to
discover an old bug in the way we use the osculating orbits produced
by the Eckstein-Hechler analytical propagator. This propagator takes
zonal terms into account, and produces directly circular parameters a,
ex, ey, ... When we compute anything related to geometry, we compute
Cartesian coordinates using the Orbit getPVCoordinates method. As the
Orbit classes do not know anything about the perturbation, the (P, V)
pair does in fact implicitly relies on Keplerian-only expressions. So
the velocity part is *not* consistent with the derivative of the
position. The real derivative of the position takes the non-Keplerian
effects into account which are ignored by getPVCoordinates. The
difference is small, but as the tests threshold were deliberatly very
tight, the tests started to fail when the various pointing directions
were not computed anymore from finite differences mainly involving
position and when they relied on the computed velocity. So the problem
already happens in the master branch, it is not specific to the
introduction of acceleration (it was just detected here during

The solution is in fact quite simple. If an orbit has been produced by
a non-Keplerian propagator, the propagator already knows about the
derivatives of the orbital elements (which are circular in the
Eckstein-Hechler model case but can be any kind of parameters for
other propagators). The propagator should therefore provide these
derivatives to the orbit so they can be used in the PVCoordinates
conversion. The code is very simple and straightforward. I have
checked this and got very interesting results with
Eckstein-Hechler/Circular, as for example a simple interpolation over
a 900s arc with proper velocity/acceleration has a 88m error with two
base points now whereas it was 5162 m before (and 0.02m vs 650m for 3
points, 1.0e-5m vs 259m for 4 points).

Here is what bothers me:

Should we create specialized classes for perturbed orbits or should we
simply add a constructor to the existing orbits with the parameters
derivatives and set them to 0 when they are not known?

For my tests, I created PerturbedCircularOrbit which extends
CircularOrbit and override the protected initPVCoordinates method and
the public shiftedBy and interpolate methods. I could also have simply
moved everything into CircularOrbit with a new constructor.

I do not like much the PerturbedXxxxOrbit approach, as it forces to
create also additional entries in the OrbitType enum with additional
converters and it becomes awkward if for example a user configures a
NumericalPropagator to generate XxxxOrbit, despite this propagator
will in fact really generate PerturbedXxxOrbit because it is what a
Numerical propagator is for. So there should be either an internal
modification of the user setting from OrbitType.XXXX to
OrbitType.PERTURBED_XXXX or an error triggered which would invalidate
*all* current user code as it would become forbiddent to generate XXXX
orbits now.

On the other hand, the drawback of modifying the existing classes to
hold the non-Keplerian derivatives is that they will consume more
memory. I don't think it is a problem with current computers.

In any case, initial orbits created directly from user code or by
reading files would not include the derivatives and therefore will be
built as usual (by calling the unmodified classes in the first
approach, or by using the already existing constructors in the second
approach, assuming these constructors will automatically set the
derivatives to Keplerian-only values). In any case, full-blown
perturbed orbits will be created internally by Orekit propagators,
which can easily be modified to provide the derivatives they know (by
creating instances of the new derived classes in the first approach,
or by using new constructors with additional parameters in the second

My humble opinion would be to use the second approach to solve this
bug. I will probably do this in the position-velocity-acceleration
branch so it will include accelerations right from the start and will
be merged to master at the same time as the rest of the branch. Of
course, this will be a dedicated commits (Git branches are great!).

What do you think ?

best regards,

This message was sent using IMP, the Internet Messaging Program.

This message was sent using IMP, the Internet Messaging Program.

This message was sent using IMP, the Internet Messaging Program.