Re: [Orekit Developers] status on second order derivatives

[Evan Ward <evan.ward@nrl.navy.mil>]

On 11/10/2014 11:51 AM, MAISONOBE Luc wrote:
> Evan Ward <evan.ward@nrl.navy.mil> a écrit :
>
>>
>> On 11/01/2014 04:23 AM, MAISONOBE Luc wrote:
>>> Hi Evan,
>>> Evan Ward <evan.ward@nrl.navy.mil> a écrit :
>>>
>>>>
>>>> On 10/30/2014 07:10 AM, Luc Maisonobe wrote:
>>>>> Hi Evan,
>>>>>
>>>>> Le 29/10/2014 21:34, Evan Ward a écrit :
>>>>>> On 10/29/2014 03:14 PM, MAISONOBE Luc wrote:
>>>>>>> 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
>>>>>>>> correctly?
>>>>>>> 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.
>>>>>> I guess I don't understand why we need more than 6 parameters to
>>>>>> specify
>>>>>> the osculating Keplerian orbit.
>>>>> For Keplerian orbit, 6 parameters are enough and the relationship
>>>>> between (a, e, i, pa, raan, M) and (x, y, z, vx, vy, vz) is a clear
>>>>> bijection. If we restrict ourselves to Keplerian motion, this
>>>>> relationship holds not only at t0 but also at t0+dt for every
>>>>> value we
>>>>> want. If we use Keplerian parameters evolution (i.e. only M evolving)
>>>>> and recompute at several steps the corresponding Cartesian, we
>>>>> will see
>>>>> vx, vy, vz are really the derivatives or x, y, z, because the
>>>>> relationship does take Kepler laws into account.
>>>>>
>>>>> The problem arise when we introduce non-Keplerian parts. Then if we
>>>>> blindly apply the same relationship for converting between the two
>>>>> sets
>>>>> of parameters, we see that they are not consistent anymore.
>>>>>
>>>>> What happens is that the mapping:
>>>>>
>>>>>   (a, e, i, pa, raan, M)     <---> (x, y, z, vx, vy, vz)
>>>>>
>>>>> is such that when da/dt = 0, de/dt = 0, ... dM/dt = n, then the
>>>>> computed
>>>>> vx is dx/dt, vy is dy/dt and vz is dz/dt.
>>>>>
>>>>> If now we decide to set arbitrarily da/dt to some made up non null
>>>>> value
>>>>> and so on for all parameters, we still compute the same 6 values
>>>>> (x, y,
>>>>> z, vx, vy, vz) as before (the initial value of a, e, i ... did not
>>>>> change), but now differentiating the mapping shows that vx will
>>>>> not be
>>>>> equal to dx/dt anymore, vy will not be dy/dt anymore and vz will
>>>>> not be
>>>>> dz/dt anymore. The made up da/dt, de/dt will contribute to the
>>>>> results.
>>>>> This is what happens for Eckstein-Hechler as it computes the left
>>>>> part
>>>>> of the mapping and we (Orekit team) compute the right part using the
>>>>> Keplerian-based mapping. If we want proper computation of vx, vy,
>>>>> vz, we
>>>>> need the additional parameters.
>>>>
>>>> I agree that we need to account for the element rate terms to
>>>> correctly
>>>> compute the velocity. I think we could do the (a, e, i, pa, raan, M,
>>>> rates, other propagator state) --> (p, v) mapping within the
>>>> Propagator
>>>> and then return a Orbit built from the correctly computed PV. This way
>>>> the returned orbit would be the osculating Keplerian orbit that
>>>> matches
>>>> the propagator's position and velocity. If the user then translates
>>>> the
>>>> osculating PV to osculating Keplerian elements they will be different
>>>> from the propagator's elements, but that is to be expected since the
>>>> propagator includes perturbations and osculating elements don't.
>>>> (Maybe
>>>> this is the source of my confusion; when hear "osculating orbit" I
>>>> think
>>>> of a conic orbit that matches the instantaneous PV at a single instant
>>>> in time.)
>>>
>>> OK, I think this time I understood your concerns.
>>> I'll give this a try, were the only modification to CircularOrbit will
>>> be to add a constructor that accepts both the circular and Cartesian
>>> parameters, hence the class will not do the conversion by itself
>>> anymore and only the propagator need to be aware of the circular
>>> derivatives.
>>>
>>> My current setting was to use either circular or Cartesian but never
>>> build an orbit from both, so the conversion could be done only by the
>>> orbit class and therefore it should have known about the circular
>>> derivatives. With your proposal, it is simpler. Thanks for the
>>> suggestion.
>>
>> Why not just return a CartesianOrbit from the EHPropagator? Something
>> like:
>>
>> //in EHPropagator.propagateOrbit
>> PVCoordinates pv = ...; //correctly compute instantaneous PV
>> return new CartesianOrbit(pv, frame, date, mu);
>>
>> This way no modifications are necessary for the Orbit classes, and the
>> user gets the correct PV at the requested time. One caveat is that if
>> the user needs circular parameters they will have to convert using the
>> CircularOrbit(orbit) constructor, and the circular parameters under the
>> new scheme would not match those returned by the current implementation.
>
> This looked very interesting, so I tried it ... see below.
>
>>
>> Maybe this is a good time to add some documentation to the
>> Propagator/SpacecraftState classes to specify exactly how the returned
>> orbit is defined. The way I see it we have three options:
>>
>> 1. The current state in EH, where the EH parameters are used to create a
>> Keplerian orbit. The user sees the EH parameters, but the velocity and
>> local extrapolation is incorrect. This is what we are trying to replace.
>
> Yes. It was wrong and it should be fixed.
>
>>
>> 2. Create and return a conic orbit computed from the propagator's
>> position and velocity. This matches the definition of osculating orbit
>> provided several authors. (e.g. Battin, Montenbruck & Gill) I think this
>> approach is the simplest for the user to understand and the simplest to
>> implement. (This is the option I outlined above.)
>
> I tried it. It is very simple to implement and removes lots of code
> with respect to the third option below.
>
> The problems I got after making this change is that the orbit returned
> by the propagator seems different from the initial orbit used to
> create the propagator (or used in resetInitialState). Typically, the
> unit tests sameDateKeplerian, sameDateCartesian and
> propagatedEquinoctial in EcksteinHechlerPropagatorTest fail with semi
> major axis being a few hundreds of meters away from expected value,
> despite position was the same at 2.5e-8m. Of course this can be
> expected as we know velocity is inconsistent.
>
> However, as the input osculating orbit should be considered perfect I
> don't know how to handle this.
>
> I first tried to force the mean parameters computation to converge to
> the same initial parameters, which worked great for the specific tests
> outlined above but created very wrong propagated orbits. The
> propagation comparison with numerical propagators
> HolmesFeatherstoneAttractionModelTest.testEcksteinHechlerReference
> failed with big errors (a few thousand meters after 12h and diverging,
> when the expected error is periodic and stable at roughly 100m).

This result surprised me since I would expect that two propagators that
purport to implement the same model would produce the same ephemeris
given the same initial position and velocity. After experimenting with
your code for a while I think that the 100 meter differences we see with
options 1 and 3 are the cause of the 9km differences we see with option
2. When the EH propagator determines its mean orbit it either fits EH
p/v to the initial p/v (option 2) or it fits the EH elements to the
initial Keplerian elements (option 3). Option 3 seems to fit the orbit
as a whole better, but the EH velocity and the initial velocity are
different by a small amount, where as option 2 forces the velocity to
match exactly even when the EH velocity is in error by a small amount.
For example in  testEcksteinHechlerReference the velocity difference
between the two options at the initial epoch is ~0.15 m/s which seems
large enough to produce the observed km of difference over 50 000 s.

To check that the initial orbit matching is causing the errors, I found
that if the initial orbit is shifted by a quarter period (or -0.75
period) that the EH ephemeris stays within 70 m of the numerical propagator.

So IMHO we can either fit the velocity exactly at the epoch or fit the
orbit parameters and tolerate a small velocity discrepancy at the epoch
in return for an orbit closer to the numerical propagator. Both can work
with option 2.

Do you have access the reference publication for the EH propagator? I
think it may help to see how the propagator is initialized there.

Best Regards,
Evan

>
> So up to now, I was not able to have both initial orbit, propagated
> orbit and velocity correct. In all my tests, I had only two correct
> among the three.
>
>>
>> 3. Create a combined orbit that has both the EH elements and the
>> position and velocity. This may be a bit confusing for the users since
>> the EH orbital elements are only define w.r.t. the EH propagation model.
>
> In fact, as seen in the previous case, the orbital elements are also
> visible to the user as they are used to initialize the propagator from
> an osculating orbit (notwithstanding velocity). The hidden parameters
> are the mean parameters from wich EH computes the osculating ones.
>
>> In other words, if a user writes a method to use a Propagator (the base
>> interface, not one of its sub classes) then the user wouldn't be able to
>> use the elements (a, e, i, ...) in the returned orbit since they
>> wouldn't know how the elements are defined. On the other hand, if the
>> user needs the EH parameters then I think this is the definition we
>> should choose. (IIUC this is what Luc is describing.)
>
> I tried it for the last few days (before switching to your suggestion
> number 2 above). It almost works but is quite complicated and I did
> not yet manage to solve a last problem.
>
> The way I've done it was to create an internal EHOrbit class extending
> CircularOrbit and overriding a few methods. At construction time it
> preserve the value of the circular osculating parameters (so initial
> orbit is preserved when propagating to same date), and it does not
> modify the mean parameters computation loop (so propagated orbit is
> consistent with other numerical propagators), and it additionally
> computes Cartesian parameters knowing about the derivatives (so it
> solves the initial problem). It's a bunch of code, but is quite
> straightforward thanks to the DerivativeStructure class which greatly
> simplifies derivatives computation.
>
> This however is not enough, and a few methods should be overriden too:
> shitfetBy and interpolate. Failing to overriding thme creates another
> type of inconsistency: the velocity of shifted or interpolated
> instances computed directly by the base CircularOrbit class is wrong.
> Overriding shiftedBy is trivial and works well. Overriding interpolate
> is difficult as it could be fed from several types of orbit in the
> sample, so I tried to rebuild the derivatives if the sample did not
> include them. This is a nightmare and I utterly failed. The problem is
> interpolate is automatically used for example when one creates an
> ephemeris from an EcksteinHechler propagator. It is a feature shared
> among all propagators and widely used.
>
> So I am confused now.
>
> I could attempt to reset the specialized EHOrbit class that knows
> about non-Keplerian effect when converting back and forth between
> circular and Cartesian, but I have yet to find how interpolate should
> work. Perhaps we could simply say that the specific implementation of
> interpolate in EHOrbit works *only* when the sample contains only
> EHOrbit instances and throws an exception otherwise. It seems fair to
> me as we cannot rebuild missing derivatives, and it would still work
> correctly in the regular case where people use consistently EHOrbit
> for a sample. In fact, it could probably be considered a general
> contract of the interpolate method to use only one type of orbits at a
> time and not to convert anything under the hood.
>
> Any thought ?
> Luc
>
>>>
>>>
>>>>
>>>>> The problem may not hold for NumericalPropagator (I still have to
>>>>> check), because basically we do the computation the other way
>>>>> round. We
>>>>> start from x, y, z, vx, vy, vz and deduce a, e, ... from the
>>>>> mapping. So
>>>>> as long as our vx=dx/dt, vy=dy/dt, vz=dz/dt, the initial mapping will
>>>>> still hold.
>>>>>
>>>>>> It seems that if we add rates to all the elements then the resulting
>>>>>> class could be classified as a Propagator.
>>>>> It is only a way to get a consistent (P, V), and as a propagator it
>>>>> would be really limited as it is a Taylor expansion. I would simply
>>>>> qualify it as a perturbed orbit allowing local expansion.
>>>>
>>>> Extending the Taylor series expansion to higher order derivatives
>>>> makes
>>>> sense if we need to return the propagator's internal state, and the
>>>> satellite's PV in one consistent object. Is there a use case for
>>>> providing the propagator's internal state? I don't know, but I can't
>>>> think of any.
>>>>
>>>> If we need to provide the user with acceleration as well as PV then I
>>>> agree that we need a bigger container class to hold the extra
>>>> information. Maybe the new PVACoordinates would be the right choice
>>>
>>> I agree this is a good place. However, I did not create an extra class
>>> for this either, I simply added a third vector to PVCoordinates. We
>>> can discuss this later, when the merge of the branch will be attempted.
>>>
>>>> instead of a PerturbedOrbit class since I don't think we would get the
>>>> same extrapolation from a PerturbedOrbit stored in Keplerian vs
>>>> Cartesian elements. To put it another way, If I converted a
>>>> PerturbedKeplerianOrbit to a PerturbedCartesianOrbit would they both
>>>> follow the same path?
>>>
>>> At second order level, in the neighborhood of the matching date, yes
>>> they should. Differences should appear at the first ignored derivative
>>> level and build up from there.
>>>
>>>>
>>>> I'm not against the PerturbedOrbit classes; I just want to make
>>>> sure we
>>>> don't choose a complex solution to a simple problem.
>>>
>>> You are right. I have a clear tendency to over-engineer things, this
>>> is why I ask for advice from time to time. This is the great force of
>>> collaborative development : people don't get stuck in their own errors
>>> too long.
>>
>> I think these constructive discussions combine the best ideas from
>> everyone.
>>
>>>
>>> So to be clear, here is the consensus :
>>>
>>>  - I will simplify what I have done to limit as much as possible
>>>    the changes to CircularOrbit.
>>>  - In this setting, orbits will not know about non-Keplerian effects.
>>>  - As this knowledge is needed for accurate and more importantly
>>>    consistent velocity computation, the mapping between orbital
>>> parameters
>>>    and Cartesian parameters will be precomputed by the propagators
>>>  - If the mapping is not precomputed and orbits are built from one
>>>    set of parameters only, then the mapping will be Keplerian only
>>>    and will be performed byt orbit classes just as it is now
>>>
>>
>> If we decide to define the orbit returned from a Propagator according to
>> option 3 above then I think this is a good plan to implement it.
>>
>> Best Regards,
>> Evan
>>
>>> best regards,
>>> Luc
>>>
>>>>
>>>>>> There might be a case for
>>>>>> this intermediate level of fidelity/speed, but is it worth the added
>>>>>> complexity? Especially since we already have one that is fast (the
>>>>>> Orbit
>>>>>> class) and one that is precise (the Propagator classes). Maybe I'm
>>>>>> still
>>>>>> not understanding your proposal.
>>>>> I detected the problem when checking some pointing attitude modes
>>>>> that
>>>>> are related to spacecraft velocity (alignment of the spacecraft axis
>>>>> with ground drift for Earth observation). This mode transforms
>>>>> something
>>>>> that is a derivative (velocity) into a regular non-differentiated
>>>>> value
>>>>> (an angle). So we go up one order of derivation. For such modes, a
>>>>> wrong
>>>>> velocity induces a wrong angle and a wrong acceleration induces a
>>>>> wrong
>>>>> angular rate. As the angular rate is used for time shift, I needed
>>>>> it to
>>>>> be accurate.
>>>>>
>>>>> You also mentioned issues with Doppler which would also occur.
>>>>>
>>>>> We do need accurate velocity, and we do need a velocity that is
>>>>> consistent with the derivative of the position.
>>>>
>>>> Definitely agree that getting the velocity correct is important.
>>>> Thanks
>>>> again for all the work you've put into solving this bug.
>>>>
>>>> Best Regards,
>>>> Evan
>>>>
>>>>>
>>>>> best regards,
>>>>> Luc
>>>>>
>>>>>> Best Regards,
>>>>>> Evan
>>>>>>
>>>>>>>> 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,
>>>>>>> Luc
>>>>>>>
>>>>>>>> Best Regards,
>>>>>>>> Evan
>>>>>>>>
>>>>>>>> On 10/29/2014 06:57 AM, MAISONOBE Luc wrote:
>>>>>>>>> Hi Paul,
>>>>>>>>>
>>>>>>>>> paulcefo <paulcefo@buffalo.edu> a écrit :
>>>>>>>>>
>>>>>>>>>> Luc,
>>>>>>>>>>
>>>>>>>>>> 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
>>>>>>>>> propagator.
>>>>>>>>> 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,
>>>>>>>>> Luc
>>>>>>>>>
>>>>>>>>>> Paul
>>>>>>>>>>
>>>>>>>>>> -- 
>>>>>>>>>> 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
>>>>>>>>>> USA
>>>>>>>>>>
>>>>>>>>>> 508-696-1884 (phone on Martha's Vineyard)
>>>>>>>>>> 978-201-1393 (cell)
>>>>>>>>>>
>>>>>>>>>> paulcefo@buffalo.edu
>>>>>>>>>> paul.cefola@gmail.com
>>>>>>>>>>
>>>>>>>>>> On 10/29/2014 6:02 am, MAISONOBE Luc wrote:
>>>>>>>>>>> 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.
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> ----------------------------------------------------------------
>>>>>>>>> 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.
>>>
>>>
>>
>>
>
>
>
> ----------------------------------------------------------------
> This message was sent using IMP, the Internet Messaging Program.
>
>