From Osculating to Mean Elements with Brouwer-Lyddane Propagator, negative eccentricity

changed the description

Thanks a lot for the information and the test, @alberto-ferrero. I could reproduce the bug and wanted to give you an update.

The issue appears when the eccentricities are too little within the while-loop in ComputeMeanParameters. To give you an idea, in your test, the eccentricity of the osculating orbit that is fed to the method is already 0.000958. Indeed, at some point during the iterations, the eccentricity of the orbit computed by the algorithm drops below zero and the KeplerianOrbit constructor rises an exception. I had a look at the equations and everything looks fine to me.

In ComputeMeanParameters there is already a guard for negative eccentricities that resets the value to zero when a KeplerianOrbit is created, but that does not extend to all the times the KeplerianOrbit constructor is called (hence the exception). Nevertheless, extending this guard does not solve the issue, as the algorithm would not converge anyway.

I have noticed that resetting the eccentricity to the one of the initial osculating orbit (instead of zero) would let the algorithm converge for this specific test. However, that does not solve the issue completely.

A possible solution would be replacing the KeplerianOrbit with CartesianOrbit (or EquinoctialOrbit). Unfortunately this is not straight forward as the algorithm will need to be adapted in several points.

I keep investigating this issue, but I am open to comments or suggestions.

Hi @gc thanks for the update.

I have noticed that resetting the eccentricity to the one of the initial osculating orbit (instead of zero) would let the algorithm converge for this specific test. However, that does not solve the issue completely.

This is a nice hint. I can try to temporarily work from here as intermediate solution, waiting for your investigation to be completed. Much appreciated the effort and the updates. Alberto

@gc Hi Do we have any plan to fix this? Thanks Alberto

mentioned in issue #1122 (closed)

New record of the issue with the addition of a StepHandler: https://forum.orekit.org/t/negative-eccentricity-computed-while-propagating-with-brouwer-lyddane-algorithm/2669?u=bcazabonne

mentioned in merge request !382 (closed)

assigned to @vcucchie

I've taken a look at it and i would like to wish good luck to anyone tackling this issue.

Also, i've noticed an infinite cycle when propagating the test provided in the issue description to Constants.JULIAN_DAY * 0.1.

Cheers, Vincent

Hi @gc Going back to this topic. I am trying to follow your suggestion to return a different orbit type. If I consider EquinoctialOrbit, would be something like this?

// Return a Equinoctial orbit
final double ex = e.getValue() * FastMath.cos(g.getValue() + h.getValue());
final double ey = e.getValue() * FastMath.sin(g.getValue() + h.getValue());
final double hx = FastMath.tan(i.getValue() / 2.0) * FastMath.cos(h.getValue());
final double hy = FastMath.tan(i.getValue() / 2.0) * FastMath.sin(h.getValue());

Basically I derivate ex -> exDot, ey -> eyDot, etc.

d/dt(e cos(g(t) + h(t))) = -e (g'(t) + h'(t)) sin(g(t) + h(t))
d/dt(e sin(g(t) + h(t))) = e (g'(t) + h'(t)) cos(g(t) + h(t))
d/dt(e tan(i(t)/2) cos(h(t))) = 1/2 e cos(h(t)) i'(t) sec^2(i(t)/2) - e h'(t) sin(h(t)) tan(i(t)/2)
d/dt(e tan(i(t)/2) sin(h(t))) = e h'(t) cos(h(t)) tan(i(t)/2) + 1/2 e sin(h(t)) i'(t) sec^2(i(t)/2)

So I can build:

return new EquinoctialOrbit(a.getValue(), ex, ey, hx, hy, l.getValue(),
                a.getFirstDerivative(), exDot, eyDot, hxDot, hyDot, l.getFirstDerivative(),
                PositionAngleType.TRUE, mean.getFrame(), date, mu);

What do you think? This implementation solves my bug, but many other tests in orekit\src\test\java\org\orekit\propagation\analytical\BrouwerLyddanePropagatorTest are actually failing

What about this approach for CartesianOrbit? Should I get X,Y,Z from Keplerian elements and derivate?

Thanks Alberto

Not sure it's worth tagging GC, I don't think he's working on this anymore.

In your code Alberto, did you recompute l as well? Because it's a true anomaly in the original code, not an argument of longitude.

unassigned @vcucchie

oh, @Serrof thanks! yes indeed, I was considering the true anomaly, not true longitude.

//previous code untouched
// True anomaly
final UnivariateDerivative2 l = FastMath.atan2(B, A);

// Argument of perigee
final UnivariateDerivative2 g = g_p_l.subtract(l);

// True longitude
final UnivariateDerivative2 Lv = l.add(g).add(h);

final double sma = a.getValue();
final double ex = e.getValue() * FastMath.cos(g.getValue() + h.getValue());
final double ey = e.getValue() * FastMath.sin(g.getValue() + h.getValue());
final double hx = FastMath.tan(i.getValue() / 2.0) * FastMath.cos(h.getValue());
final double hy = FastMath.tan(i.getValue() / 2.0) * FastMath.sin(h.getValue());
final double Lv = lv.getValue();
final double smaDot = a.getFirstDerivative();
final double exDot = e.getFirstDerivative() * FastMath.cos(g.getValue() + h.getValue())
 - e.getValue() * (g.getFirstDerivative() + h.getFirstDerivative()) * FastMath.sin(g.getValue() + h.getValue());
final double eyDot = e.getFirstDerivative() * FastMath.sin(g.getValue() + h.getValue())
 + e.getValue() * (g.getFirstDerivative() + h.getFirstDerivative()) * FastMath.cos(g.getValue() + h.getValue());
final double hxDot = (FastMath.cos(h.getValue()) * i.getFirstDerivative() - i.getValue() * h.getFirstDerivative() * FastMath.sin(h.getValue())) / (FastMath.cos(i.getValue() * FastMath.cos(h.getValue())) + 1);
final double hyDot = (i.getValue() * h.getFirstDerivative() * FastMath.cos(h.getValue()) + FastMath.sin(h.getValue()) * i.getFirstDerivative()) / (FastMath.cos(i.getValue() * FastMath.sin(h.getValue())) + 1);
final double LvDot = lv.getFirstDerivative();

return new EquinoctialOrbit(sma, ex, ey, hx, hy, Lv,
                smaDot, exDot, eyDot, hxDot, hyDot, LvDot,
                PositionAngleType.TRUE, mean.getFrame(), date, mu);

But still, no regression is failing (an actually also the circular test of this thread does not converge as well.

I see that despite some tests failing because of larger error (so maybe need to adjust the tolerance), some of them does not converge at all. This is wierd. Can it be because of the further conversion Equinoctial -> Keplerian during the call computeMeanParameters?

I will investigate further.

Thanks for your comment.

As I said on the forum, the method propagateParameters is suspicious. It uses UnivariateDerivative2 but using UnivariateDerivative1 would have exactly the same effect. It looks like a misuse of Derivative actually.

@Serrof thanks for coming back to this. This is the comment right #1292 (closed) I will investigate

Hi Alberto, yes it is. Thanks for looking at this

So I fixed the convergence:

proper derivatives on Equinoctial elements to build the mean
usage of Keplerian to calculate residuals during iteration

I realized I was considering the Mean Longitude as True Longitude, so thanks to your issue #1291 (closed) I could resolve it

Opened MR !484 (closed)

mentioned in merge request !484 (closed)

changed milestone to %12.1.1

assigned to @alberto-ferrero

@alberto-ferrero let's wait for 12.1 to be out so we can create a clean MR. I'll help out. Cheers

From Osculating to Mean Elements with Brouwer-Lyddane Propagator, negative eccentricity

Designs

Child items ...

Activity