Using NRLMSISE-00 with numerical propagator

dl00718 · August 10, 2022, 8:46am

Thank you all for helping out. I think the time step is the issue as @MaximeJ suggests, because reducing the step-size and position tolerance makes it work, but then we lose the advantange of the semi-analytical propagator.

@bcazabonne The error doesn’t occur when i use a numerical propagator, but it does with the DSST propagator using both Harris Priester and JB2008 atmospheric models. @evan.ward I believe this user had a similar issue https://forum.orekit.org/t/dsst-propagator-issues-at-low-altitudes/640

Here’s the last set of orbital elements written to my output file (Angles in degrees, sma in Km).
Date : 2021-06-16T11:55:08.73936Z
SMA : 6552.9290
Ecc: 0.0003071
Incl: 22.97438
RAAN: 162.40034
AoP: -185.05898
MAN: 58.98348

@evan.ward I have attached the test case and a unit test, and heres a copy of the backtrace from the exception.

Exception in thread "main" org.orekit.errors.OrekitException: altitude (-5,432,513.404 m) is below the 100,000 m allowed threshold
	at org.orekit.models.earth.atmosphere.HarrisPriester.getDensity(HarrisPriester.java:283)
	at org.orekit.models.earth.atmosphere.HarrisPriester.getDensity(HarrisPriester.java:393)
	at org.orekit.forces.drag.DragForce.acceleration(DragForce.java:80)
	at org.orekit.propagation.semianalytical.dsst.forces.AbstractGaussianContribution$IntegrableFunction.value(AbstractGaussianContribution.java:1003)
	at org.orekit.propagation.semianalytical.dsst.forces.AbstractGaussianContribution$GaussQuadrature.basicIntegrate(AbstractGaussianContribution.java:1472)
	at org.orekit.propagation.semianalytical.dsst.forces.AbstractGaussianContribution$GaussQuadrature.integrate(AbstractGaussianContribution.java:1383)
	at org.orekit.propagation.semianalytical.dsst.forces.AbstractGaussianContribution.getMeanElementRate(AbstractGaussianContribution.java:388)
	at org.orekit.propagation.semianalytical.dsst.forces.AbstractGaussianContribution.getMeanElementRate(AbstractGaussianContribution.java:300)
	at org.orekit.propagation.semianalytical.dsst.DSSTPropagator$Main.elementRates(DSSTPropagator.java:1173)
	at org.orekit.propagation.semianalytical.dsst.DSSTPropagator$Main.computeDerivatives(DSSTPropagator.java:1152)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator$ConvertedMainStateEquations.computeDerivatives(AbstractIntegratedPropagator.java:758)
	at org.hipparchus.ode.ExpandableODE.computeDerivatives(ExpandableODE.java:134)
	at org.hipparchus.ode.AbstractIntegrator.computeDerivatives(AbstractIntegrator.java:265)
	at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:252)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:477)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:425)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:385)
	at test.reentry(test.java:100)
	at test.main(test.java:106)

Thank you
test.java (5.7 KB)
testTest.java (429 Bytes)

evan.ward · August 10, 2022, 12:05pm

Thanks for the test case! That works well.

I made maxCheck smaller and it did not change the result of the test. So I believe that means we can rule out the event detector as the issue.

From running the test and the back trace I think Maxime’s theory is correct. I’m guessing it is not an issue for NumericalPropagator because the step sizes are much smaller.

So how to fix it? One idea would be to borrow an idea from how maneuvers are modeled. For maneuvers a whole integration step is taken with the thruster on, then the step is later truncated to the time the thruster turned off. This produces better behavior with the integrator. We could do the same for drag forces. Do not throw an exception when computing acceleration. Perhaps just continue using the density at the minimum altitude of the model for the altitudes below it. Then include a force model event detector that will truncate the step. This would allow propagation down to the minimum altitude. Or the event detector could throw an exception to end with an error. Any user event detectors would be processed first because events are handled in chronological order.

Regards,
Evan

dl00718 · August 10, 2022, 1:07pm

Hi @evan.ward, thanks for the feedback.

I believe so too, changing the maxCheck doesn’t have an effect and it depends on the step sizes, which is why it works with the numerical propagator.

How would I go about implementing this, Is there an example you could reference?

Thank you

evan.ward · August 10, 2022, 2:22pm

Good question. The contributing guide is here: Orekit – Contributing to Orekit

Making an issue on GitLab is a good first step. I think you would modify the the getEventDetectors() in DragForce and DsstDragForce to return an event detector for the minimum altitude. Then modify all the implementations of Atmosphere to return a density value instead of throwing an exception for invalid altitudes. Might need to add a method to Atmosphere to query the minimum altitude. Alternatively one could modify the DragForce implementations to catch altitude too low exceptions and return a reasonable acceleration.

I’ll review your merge request when you submit it. Thanks for your contribution!

MaximeJ · August 10, 2022, 2:55pm

Thank you @evan.ward for the investigation and @dl00718 for the test and information,

Evan, your solution would indeed fix the problem !

I’ve played around a bit with dl00718 code and found the following info that can help.
I added a step handler

        dsstProp.getMultiplexer().add(new OrekitStepHandler() {
            
            int nStep = 0;
            
            @Override
            public void handleStep(OrekitStepInterpolator interpolator) {
                SpacecraftState s = interpolator.getCurrentState();
                double stepSize = s.getDate().durationFrom(interpolator.getPreviousState().getDate());
                System.out.format(Locale.US, "Step %d - dt: %10.3f days / step: %10.3f s /  Alt: %10.3f km\n", 
                                  nStep++, s.getDate().durationFrom(initDate) / Constants.JULIAN_DAY, stepSize,
                                  (s.getPVCoordinates().getPosition().getNorm() - Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 1000.);                
            }
        });

Then I’ve tried some (minStep / position tolerance) configurations for the integrator:

minStep = 3600.
dP = 10.

→ (default config) Crashes with HarrisPriester minimum altitude exception

minStep = 3600.
dP = 0.1

→ Crashes with Hipparchus ODE message: MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION. The integrator tries to reduce the step to fulfill the constraint of 0.1m on the position tolerance but cannot because the minStep is too big.

minStep = 60.
dP = 0.001

→ Working combo: After around 1540 days, the altitude detector kicks in with no crash. Last step is around 60s and truncated at ~41s by the detector.

So, ultimately, it is a configuration issue of the integrator.
Reducing minStep and dP does not impair DSST performances because the maximum step size of 2 days is used until altitude gets below ~270km (after 1513 days of propagation) and the integrator starts reducing the steps.

As a user, I would like the integrator to reduce the step duration until the minimum step size is reached and throw an exception like in the second case above.
That way I would know what to do to improve the configuration, which I found is not explicit with the first case (altitude is too low).

To me, this looks like the multiple questions we had on the forum with an exception thrown due to a negative eccentricity (see this thread for example, though it is related to initial time step only and can be fixed differently).

So I was wondering if, on top of Evan’s proposition, we couldn’t implement a mechanism that would invalidate an integrator step given some exceptions (altitude too low, negative eccentricity etc.).
That way the integrator would try to reduce the step until the minimum is reached instead of blindly throwing the exception and leaving the user in the dark.
That would require implementing a specific exception type extending OrekitException and catching it at integrator level.
I haven’t checked if it is doable with current Hipparchus or Orekit APIs though.
Do you think that it makes sense ?

evan.ward · August 10, 2022, 3:21pm

Interesting proposition. How would the integrator know where to end the step?

If it is guess and check that would be slow. Perhaps include that information in the specialized exception? But the case above the end of the step is based on altitude and to compute the time of the altitude crossing an event detector would be needed. In either case the new time to end the step may cause the step duration to be less than minStep, which would cause the integrator to throw an exception. I don’t know the answers…

dl00718 · August 11, 2022, 8:09am

Hi @MaximeJ,

This configuration works for me, thank you very much.

@evan.ward I’ll give this a try and see if i can implement it. Thank you

MaximeJ · August 11, 2022, 1:02pm

Hi @dl00718,

You are welcome !
And thank you for the upcoming contribution !

Good question, I don’t have the answer yet…

I think that is the best way to do it !

That would be, in my opinion, a proper behavior. The exception is explicit and the user can easily fix it by lowering the minimum step.

It’s true that we will probably need a detector to find the proper step. This won’t be the case for the “negative eccentricty” exception for example.
Maybe we could imagine a mix of your solution and something along the lines of what I’ve tried below.

Probably yes, still I’ve tried a very rough implementation of this in EmbeddedRungeKuntaIntegrator#integrate where I divide current step by 2 every time computeDerivatives throws an ALTITUDE_BELOW_ALLOWED_THRESHOLD exception.
Until it gets below minimum step. In that case, I try the minimum step once and if the altitude is still too low, then it throws a MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION exception.

This gives the following results:

minStep = 3600.
dP = 10.

→ (default config) This one still crashes but this time with a MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION exception.
Log below (“Model Exception” in that case is the ALTITUDE_BELOW_ALLOWED_THRESHOLD exception)

Step 771 - dt:   1540.339 days / step: 172800.000 s /  Alt:    176.307 km
	Model Exception caught - Previous step: 172800.0 / New step: 86400.0
Step 772 - dt:   1541.196 days / step:  74006.216 s /  Alt:    152.422 km
	Model Exception caught - Previous step: 106601.41359396781 / New step: 53300.706796983904
	Model Exception caught - Previous step: 53300.706796983904 / New step: 26650.353398491952
	Model Exception caught - Previous step: 26650.353398491952 / New step: 13325.176699245976
Step 773 - dt:   1541.350 days / step:  13325.177 s /  Alt:    144.398 km
	Model Exception caught - Previous step: 21702.721715022282 / New step: 10851.360857511141
Step 774 - dt:   1541.434 days / step:   7279.631 s /  Alt:    131.044 km
	Model Exception caught - Previous step: 12364.10248376542 / New step: 6182.05124188271
	Model Exception caught - Previous step: 6182.05124188271 / New step: 3091.025620941355
	Minimum step size (3600.0) reached !
	Model Exception caught - Previous step: 3600.0 / New step: 1800.0
Exception in thread "main" org.orekit.errors.OrekitException: pas minimal (3,60E03) atteint, l'intégration nécessite 1,80E03

The exception is easier to fix for the user but he won’t know that it originates from the atmosphere model anymore.

minStep = 1200.
dP = 10.

→ Works with several catching of the atmosphere model exception

Step 771 - dt:   1539.038 days / step: 172800.000 s /  Alt:    196.793 km
Step 772 - dt:   1541.038 days / step: 172800.000 s /  Alt:    162.573 km
	Model Exception caught - Previous step: 163465.3780505176 / New step: 81732.6890252588
	Model Exception caught - Previous step: 81732.6890252588 / New step: 40866.3445126294
	Model Exception caught - Previous step: 40866.3445126294 / New step: 20433.1722563147
Step 773 - dt:   1541.275 days / step:  20433.172 s /  Alt:    150.145 km
	Model Exception caught - Previous step: 32637.238297095602 / New step: 16318.619148547801
Step 774 - dt:   1541.440 days / step:  14328.235 s /  Alt:    130.133 km
	Model Exception caught - Previous step: 15690.604490723932 / New step: 7845.302245361966
	Model Exception caught - Previous step: 7845.302245361966 / New step: 3922.651122680983
	Model Exception caught - Previous step: 3922.651122680983 / New step: 1961.3255613404915
Step 775 - dt:   1541.463 days / step:   1961.326 s /  Alt:    124.568 km
	Model Exception caught - Previous step: 2357.391082375647 / New step: 1178.6955411878234
	Minimum step size (1200.0) reached !
Step 776 - dt:   1541.472 days / step:    785.688 s /  Alt:    118.806 km
2016-04-01T23:15:10.93562067578793Z - Altitude is 117.89798498219903 km and the spacecraft will re-enter soon
Re-entry event occurs on 2016-04-01T23:15:10.93562067578793Z
Step 777 - dt:   1541.472 days / step:      0.000 s /  Alt:    118.806 km
Spacecraft Re-entered