Hipparchus-Math / hipparchus

An efficient, general-purpose mathematics components library in the Java programming language
Apache License 2.0
139 stars 41 forks source link

Errors in Orekit tests when using Hipparchus 1.7-SNAPSHOT #90

Closed MaximeJo closed 4 years ago

MaximeJo commented 4 years ago

Hi all,

I ran into errors in the JUnit tests of the Orekit library (which strongly depends on Hipparchus) when linking the new version of Hipparchus (1.7-SNAPSHOT) in the pom. The errors happen in testEquatorialRetrograde methods of classes CircularOrbit and KeplerianOrbit; as well as their "fielded" versions too. The error is always on the derivative of inclination with respect to time (Idot):

java.lang.AssertionError: expected:<-4.615E-5> but was:<4.6153846153846145E-5>

The tested value is opposed to the expected one. @BryanCazabonne managed to isolate the Hipparchus commit that started this. It seems to be @6485c05"Fixed MathArrays.linearCombination when signed zeros are involved."

I don't know exactly what's happening and why the value changed. I hope one of you can find it out. Thanks in advance, Maxime

maisonobe commented 4 years ago

Hi Maxime,

The problem lies in the Orekit tests, you know who to blame for this ;-) As its names suggests, the test correspond to a retrograde equatorial orbit, which is singular for both Keplerian and circular orbits. With this state, tiny variations of velocity along the Z axis would lead to the orbital plane to cross the equatorial plane from South to North or from North to South. This implies that with the test state, the right ascension of ascending node Ω is either 0 or π. The inclination is always equal to π, but as Ω changes, its derivative changes sign too (the i=π is an angular point).

The way Ω is computed in KeplerianOrbit is by calling atan2(y, x) where in this case both y and x are 0 because they come from a cross product between orbital momentum and Z axis. Ω singularity and atan2 singularity are related to each other. For CircularOrbit the call is atan2(hy, hx) with a different meaning for the two arguments, but with similar behaviour. Atan2 properly takes into account the sign of both 0 into account to select if we should return -0, +0, -π or +π. The recent Hipparchus fix changed the output of Vector3D.crossProduct(Vector3D.PLUS_K, momentum) from (0.0, 0.0, 0.0) to (-0.0, 0.0, 0.0), i.e. the x coordinate is now a negative 0 so now Ω is equal to π instead of 0 before the Hipparchus change. As a result, the node switches from one point to the opposite point relative to the orbit, the perigee argument ω changes from π to 0 (the sum ω+Ω is stable despite each individual angle is singular), and the time derivative of inclination i changes sign.

So this is explainable: the test checks the derivative of a parameter (the inclination) just as it rebounds at its π boundary.

One way to solve the problem is to very slightly shift the orbit so that it is not perfectly equatorial and retrograde, but only very close to it. I achieved this by changing the z component of the velocity in the test cases from a perfect 0 to a tiny +1.0e-10m/s (and change the tolerance for the inclination value test from 1.0e-15 to 2.0e-14 as inclination is not π anymore, but only close to it). With this change, the test passes for both the previous Hipparchus version and the current development version. Note that if you change the velocity Z component to -1.0e-10m/s, you will see the time derivative of inclination change its sign, with seems logical to me. Note that if instead of changing the Z component of the velocity you change the Z component of the position (say to +1.0e-7m), then Ω changes to +π/2 and the test case now corresponds to a latitude maximum rather than a node, and at this point the time derivative of i drops drastically to orders of magnitudes in the 1.0e-21 domain. Once again, this shows the test as it is written is ill-conceived.

So I suggest to change the test cases using the Z velocity component as a way to ensure stable behaviour with a "nearly" equatorial retrograde orbit.

I did not check the field versions of the orbits, but I guess the same reasoning applies.

If you agree with this analysis, we can close the report as invalid.

MaximeJo commented 4 years ago

Hi Luc,

Thank you for your swift analysis and answer. I agree with you we should improve the Orekit tests. I just opened an issue on the Orekit forge in that sense. You can close the report as invalid.