[Orekit Developers] status on second order derivatives

[MAISONOBE Luc <luc.maisonobe@c-s.fr>]

Hello,

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 testing).

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  
approach).

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,
Luc

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