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.