I tried to switch 10.3 today, my numeric propagation code goes wrong when I in 10.3 or 10.2,which is not in 10.1 and before.
the Exception is as following:
org.orekit.errors.OrekitException: invalid parameter eccentricity: -0.005 not in range [0, ∞]
at org.orekit.orbits.KeplerianOrbit.checkParameterRangeInclusive(KeplerianOrbit.java:1624)
at org.orekit.orbits.KeplerianOrbit.(KeplerianOrbit.java:206)
at org.orekit.orbits.KeplerianOrbit.(KeplerianOrbit.java:164)
at org.orekit.orbits.OrbitType$4.mapArrayToOrbit(OrbitType.java:492)
at org.orekit.orbits.OrbitType$4.mapArrayToOrbit(OrbitType.java:448)
I checked release notes of 10.2，there is a bugfix about negative eccentricity. I tried to make some changes to the code to fix this. I found that when I switch the integrator to ClassicalRungeKuttaIntegrator, the code works well. On the contrary,it goes wrong with DormandPrince853Integrator( minStep_num = 0.001 maxstep_num = 20; positionTolerance = 0.1;),which is in my code before.
Did I code something wrong?
p.s.
the init orbit elements:
epoch:2018-12-27T12:00:00.000
A:6872866.0
E:0.0016002208
I:97.58207017
Raan:72.15268469
mean anormaly:59.44509482
PA:97.42471457
Hi @youngcle
Negative eccentricity should never happen, that’s why we added this exception.
Could you provide a runnable script in order to see if there is something wrong? (i.e. definition of the initial state, propagator initialization, etc.)
Specifically, could you provide the script that throws the negative eccentricity exception?
Thank you,
Bryan
the code of setting up the orbit is as following:
Blockquote
public void createKeplerianOrbit(double ORBIT_KEPLERIAN_A, double ORBIT_KEPLERIAN_E, double ORBIT_KEPLERIAN_I,
double ORBIT_KEPLERIAN_PA, double ORBIT_KEPLERIAN_RAAN, double ORBIT_KEPLERIAN_ANOMALY,
AbsoluteDate ORBIT_DATE) {
// All dates in UTC
Frame frame_orbit = FramesFactory.getEME2000();
PositionAngle angleType = PositionAngle.MEAN;
Orbit orbit = new KeplerianOrbit(ORBIT_KEPLERIAN_A, ORBIT_KEPLERIAN_E, FastMath.toRadians(ORBIT_KEPLERIAN_I),
FastMath.toRadians(ORBIT_KEPLERIAN_PA), FastMath.toRadians(ORBIT_KEPLERIAN_RAAN),
FastMath.toRadians(ORBIT_KEPLERIAN_ANOMALY), angleType, frame_orbit, ORBIT_DATE, mu);
orbit_loaded = orbit;
}
AbsoluteDate ORBIT_DATE = orbitPropagatorOrekitBugshow.convertFullStrToAbsDateTime(“2018-12-27 12:00:00”);
double ORBIT_KEPLERIAN_A = 6872866.0;
double ORBIT_KEPLERIAN_E = 0.0016002208;
double ORBIT_KEPLERIAN_I = 97.58207017;
double ORBIT_KEPLERIAN_RAAN = 72.15268469;
double ORBIT_KEPLERIAN_PA = 97.42471457;
double ORBIT_KEPLERIAN_ANOMALY = 59.44509482;
orbitPropagatorOrekitBugshow.createKeplerianOrbit(ORBIT_KEPLERIAN_A, ORBIT_KEPLERIAN_E, ORBIT_KEPLERIAN_I,
ORBIT_KEPLERIAN_PA, ORBIT_KEPLERIAN_RAAN, ORBIT_KEPLERIAN_ANOMALY, ORBIT_DATE);
there is a main in the code and when you run it,you will get exception:
Blockquote
org.orekit.errors.OrekitException: invalid parameter eccentricity: -0.007 not in range [0, ∞]
at org.orekit.orbits.KeplerianOrbit.checkParameterRangeInclusive(KeplerianOrbit.java:1624)
at org.orekit.orbits.KeplerianOrbit.(KeplerianOrbit.java:206)
at org.orekit.orbits.KeplerianOrbit.(KeplerianOrbit.java:164)
at org.orekit.orbits.OrbitType$4.mapArrayToOrbit(OrbitType.java:492)
at org.orekit.orbits.OrbitType$4.mapArrayToOrbit(OrbitType.java:448)
at
OrbitPropagatorOrekitBugshow.java (22.6 KB)
I don’t see an important error in your code. I’m just a bit surprised to see the values of the degree and order uses for ocean tides (i.e. 50x50). I’m more used to seeing values like 6x6 or 8x8.
Maybe you could also try to change the values of the max step and the position tolerance for the integrator. You could try to use 300 seconds for the max step and 10 meters for the position tolerance.
Unfortunately, I can’t do some tests in your code this week because I’m on Holidays and I don’t have my computer. I can do some runs next Monday.
Bryan
Hi @youngcle,
It’s an issue in the initialization of the step size of the Dormand-Prince 853 integrator.
It is very similar to this closed issue (628) where a quick fix is provided at the very end of the thread.
To fix your code, in method createPropagator, after the line:
AdaptiveStepsizeIntegrator integrator_adpative = new DormandPrince853Integrator(minStep_num, maxstep_num,
tolerances[0], tolerances[1]);
Add: integrator_adpative.setInitialStepSize(maxstep_num);
This will bypass the step initialization method in AdaptiveStepsizeIntegrator.initializeStep#246 and fix your issue.
The core of the problem is a bit long to explain.
To sum it up, with the initial parameters you provided (orbit, step size, tolerance etc.), the integrator is initializing a step size of h~16261s (which is way too big) and computes the new state y1 as y1 = y0 + h * dy0/dt. (y0 is the initial state you provided)
Since eccentricity e is very small and de/dt is negative here, then e(y1) is negative and an exception is raised.
Your initial parameters look fine to me, so the problem seems to come from the calculation of the initial step in Hipparchus.
Moreover, I don’t see why the initial step is allowed to get bigger than the user-defined max step size.
Maybe we should ask Hipparchus developers to limit the initial step size to the max step size. It’s a very simple check to add and it would avoid this loophole in the step size initialization. What do you think @bcazabonne and @youngcle ?
Maxime
1 Like
First of all, I want to say “thank you very much” to @bcazabonne and @MaximeJ
The problem was sovled with this:
integrator_adpative.setInitialStepSize(maxstep_num);
I checked the comment of the method SetInitialStepSize:
/** Set the initial step size.
* <p>This method allows the user to specify an initial positive
* step size instead of letting the integrator guess it by
* itself. If this method is not called before integration is
* started, the initial step size will be estimated by the
* integrator.</p>
* @param initialStepSize initial step size to use (must be positive even
* for backward integration ; providing a negative value or a value
* outside of the min/max step interval will lead the integrator to
* ignore the value and compute the initial step size by itself)
*/
in my opinion,as the base of orekit, the “org.hipparchus.ode” need not change. the numericpropagator of orekit should check this situation.
Before the code be changed, maybe the tutorial should go first. this will encourage the user community to follow your step.
Anyway, by now,I can switch to 10.3, thanks a lot. You are very kind and great guys.
Good catch!
I didn’t not remember this issue with the initial step.
I also think it could be interesting to add a check in Hipparchus to limit the initial step size to the max size. It is something we can discuss in the forum with other Hipparchus developers, in the category “Hipparchus development”.
Do you want to open a thread @MaximeJ ?
Bryan
I just opened the thread.
Hi guys,
unluckily，I found an exception show up again today,
after I add the code
integrator_adpative.setInitialStepSize(maxStep_num);
the init state as following
Propagating…
Propagator:org.orekit.propagation.numerical.NumericalPropagator@343f4d3d frame:EME2000
orbit:Init Stat orbit:
epoch:2021-04-13T12:00:00.000
A:7461477.1
E:8.27298499643803E-4
I:63.5100154876709
Raan:73.557959318161
mean anormaly:203.919706344604
PA:85.564591884613
Propagate time start:2021-04-13T16:00:00.000 end:2021-04-14T16:00:00.000
geo wgs84(init):{lat: -57.6006156639 deg, lon: 0.2284269178 deg, alt: 1,104,196.3164409997}
the exception:
org.orekit.errors.OrekitException: invalid parameter eccentricity: -0 not in range [0, ∞]
If i propagate eph from 2021-04-13T12:00:00.000 to 2021-04-14T16:00:00.000, the code goes well,no exception pop up.
Hi @youngcle,
The problem is very similar to the previous but seems to happen in the middle of the propagation instead of at the beginning.
This can be quick-fixed by replacing (l.185 of OrbitPropagatorOrekitSecondBugshow):
final OrbitType propagationType = OrbitType.KEPLERIAN;
With:
final OrbitType propagationType = OrbitType.CIRCULAR;
This will avoid the problem of negative eccentricity, and it suits well with your almost-circular orbit.
You will need to convert your final satPos_list orbits from circular to Keplerian if you want to have the same output.
I’m not fully satisfied with this answer though. I don’t see why the propagation handles differently when the starting date you ask is not the initial state date.
The only difference I’ve seen is that the propagator handler is not activated from 2021-04-13T12:00:00.000 to 2021-04-13T16:00:00.000; but that doesn’t explain the exception.
I’ll have to dig into it a little bit more, maybe open an “investigation” issue if I get stuck.
Have a good day,
Maxime
ok，i tried OrbitType.CIRCULAR，it worked.
how about OrbitType.CARTESIAN, maybe the X\Y\Z has no limit,not like “E”
yes Cartesian (or equinoctial) should work too