Discontinuity in NRLMSISE00?

Orekit usage
1

Hi!

I propagate with NRLMSISE00 for reentry purposes, and I have sometimes a problem with the integrator. I have the exception

Caused by: org.hipparchus.exception.MathIllegalArgumentException: minimal step size (1.00E-03) reached, integration needs 9.61E-04
	at org.hipparchus.ode.nonstiff.StepsizeHelper.filterStep(StepsizeHelper.java:179)
	at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:259)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:550)
	... 23 more

I have investigated, and it looks like it comes from a discontinuity in the model MSIS00. Indeed if I execute the code below, it shows that the density juste before and after midnight are not exactly the same and differ by 0.1%. With such density, this relative error is enough to force the integrator (I use DormandPrince853Integrator) to reduce its time step lower and lower until the error is raised.

Does someone know if such a discontinuity in MSIS is something that is expected?

Many thanks,

Dorian

import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.junit.jupiter.api.Test;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.LazyLoadedDataContext;
import org.orekit.models.earth.atmosphere.NRLMSISE00;
import org.orekit.models.earth.atmosphere.NRLMSISE00InputParameters;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;

public class MSISTest {

    @Test
    public void testDiscontinuity() {

        LazyLoadedDataContext dataContext = new LazyLoadedDataContext();
        BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
                                               dataContext.getFrames().getITRF(IERSConventions.IERS_2010, false));
        NRLMSISE00 msis = new NRLMSISE00(newSolarActivity(), dataContext.getCelestialBodies().getSun(), earth, dataContext.getTimeScales().getUT1(IERSConventions.IERS_2010, true));
        AbsoluteDate date1 = new AbsoluteDate(2041, 4, 17, 0, 0, 0.0, dataContext.getTimeScales().getUTC());
        AbsoluteDate date2 = date1.shiftedBy(-1e-9);
        Vector3D position = new Vector3D(5788625.05094290, 1550341.676937541, 2427500.326308472);

        System.out.println(msis.getDensity(date1, position, dataContext.getFrames().getEME2000()));
        System.out.println(msis.getDensity(date2, position, dataContext.getFrames().getEME2000()));
    }

    private NRLMSISE00InputParameters newSolarActivity() {
        return new NRLMSISE00InputParameters() {
            @Override
            public double getDailyFlux(AbsoluteDate date) {
                return 71.96;
            }

            @Override
            public double getAverageFlux(AbsoluteDate date) {
                return 71.95;
            }

            @Override
            public double[] getAp(AbsoluteDate date) {
                return new double[]{9.05, 9.05, 9.05, 9.05, 9.05, 9.04833, 9.045};
            }

            @Override
            public AbsoluteDate getMinDate() {
                return AbsoluteDate.PAST_INFINITY;
            }

            @Override
            public AbsoluteDate getMaxDate() {
                return AbsoluteDate.FUTURE_INFINITY;
            }
        };
    }
}
2

Hi Dorian,

I’m working on re-entry scenarios as well using orekit and I had encountered the same type of error where my integrator was exploding. If you don’t mind can you once check how many propagators are being built and executed simultaneously. That was the issue for me. Also you can look into using GraggBulirschStoerIntegrator, and see if your error is being solved. Let me know!

3

Hello @dorian,

Could you also precise how you initialize your integrator ?

Cheers,
Vincent

4

Hi,

I used the Dormant Prince Integrator Builder to initialize the integrator, as follow:

new DormandPrince853IntegratorBuilder(0.001, 300.0, 0.00001);

The settings on the integrator come from studies that found the best compromize between precision and computation time.

Dorian

5

Hi,

Thanks for the anwer, I’ll have a look at GraggBulirschStoerIntegrator, see if it solves my problem.

However, I still feel that it is a little bit strange that NRLMSISE00 has such discontinuities. I was expecting the solar activity to create discontinuities in the model, but the test in my first post shows that there are also temporal discontinuities in the model. And it seems that these discontinuities cause the integrator to raise error.

But maybe it is something that is expected. I am not fully familiar with this model.

Many thanks,

Dorian

6

Many atmospheric models do have discontinuities, either at day boundaries or at altitude layers boundaries. This is also true for tropospheric effects, solar activity, ionosphere, EOP…

7

Many thanks for your answer. As suggested I will try to use another integrator to avoid the exception

Have a good day!

Dorian