Using NRLMSISE-00 with numerical propagator

I’m new to orekit and I’m currently trying to add the NRLMSISE-00 atmospheric model to my drag force model so i can propagate an orbit.Not sure how to implement it so that the integrator automatically gets the atmospheric density at the integration step epoch and altitude. My current approach is shown below.

// Drag
NRLMSISE00 atm = new NRLMSISE00(null, sun, earth);
IsotropicDrag spacecraft = new IsotropicDrag(area,Cd);
DragForce drag = new DragForce(atm,spacecraft);
propagatorBuilder.addForceModel(drag);

Hi @dl00718

Welcome to the Orekit forum! :slight_smile:

Your implementation is close to be perfect. However, you missed the initialization of the provider for solar activity data. This provider is mandatory because it is used to retrieve the needed parameters used to compute the atmosphere density.

Orekit has two providers: MarshallSolarActivityFutureEstimation and CssiSpaceWeatherData. Please find below an example showing how initializing the first one.

final MarshallSolarActivityFutureEstimation msafe =
                            new MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
                                                                      MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
            final DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
            manager.feed(msafe.getSupportedNames(), msafe);

The initialization of the other one is very close. Like that you can replace your null parameter by msafe

An important information. These two provider use external data. They are available in the orekit-data folder available here https://gitlab.orekit.org/orekit/orekit-data. The orekit-data can be updated with the last available values using the update.sh script available in the folder.

I hope this will help you.

Best regards,
Bryan

Hi @bcazabonne

Thanks alot for the help. So this is how i went about it. Is it right?
File orekitData = new File(“path to orekit data folder”);
DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
manager.addProvider(new DirectoryCrawler(orekitData));
// Initialize Solar Activity Provider Data Provider
final CssiSpaceWeatherData cssi = new CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES,manager,timescale);

    // Drag
    final NRLMSISE00 atm  = new NRLMSISE00(cssi,sun,earth);
    IsotropicDrag spacecraft = new IsotropicDrag(area,Cd);

    DragForce drag =   new DragForce(atm,spacecraft);

I have another question, although it might be a silly one. When i use the tlePropagator.getInitialState() method, does it give a state vector computed from the oscullating orbital elements or do i get that after converting the tlePropagator to a numericalPropagator and calling numericalPropagator.getInitialState()?

Thank you

Hi,

Yes, it looks good!

That’s not a silly question. Calling tlePropagator.getInitialState() will give you the osculating orbital state of the satellite according to SGP4 model, which is the dynamical model used to convert TLE to orbital state. Calling numericalPropagator.getInitialState() gives you the initial state of the numerical propagator. Just be careful that after an orbit propagation, the numerical propagator updates its initial state.
Now, if you initialized the initial state of the numerical propagator using the the results of tlePropagator.getInitialState(), it should give you the same result at time t = t0.

Bryan

Hi Bryan,
Thanks alot for the explanation. So if i want to propagate an orbit from a TLE epoch to re-entry, there’s no need to convert the TLe propagator to a numerical propagator. Instead I should set up a numerical propagator and use the initial state from

as my initial state, set up an event detector, then propagate forward. Is that the right approach to take?

Thank you

Hi @dl00718,

I think you will get a large error by doing so.
See this former answer from Luc for more insights on why “instantaneous position given by a TLE propagator cannot be used as an initial orbit for another propagator, regardless of the model you use”.

You need to first convert the TLEpropagator to a NumericalPropagator, and then use this converted propagator to propagate until re-entry.
See this post or this tutorial for converting a propagator…

Hi @MaximeJ

Thank you for the clarification Maxime. I’ll implement it that way.

Hi @MaximeJ @bcazabonne,

I was reading through the FramesFactory class and it recommends not propagating in the TEME frame unless using the SGP4 propagator. Should i retrieve the initial state of the converted numerical propagator and rotate it to another frame like EME2000, then create a new propagator with the rotated initial state? Or is it alright to use the converted numerical propagator in the TEME frame.

You don’t need to, just transform the initial state and use propagator.resetInitialState(myNewState)

As you read it it’s not recommended so I would avoid it.
I wasn’t aware of that, thank you for pointing it (the explanation is page 57 of the CCSDS ODM blue book mentioned in the Javadoc)
However, ince the frame is ambiguous, your initial state converted to EME2000 will still bear this ambiguity (since we’re not 100% sure how the original TEME frame you used should be defined). So this small uncertainty will be propagated through.
That being said, I would still follow CCSDS and Orekit recommendation of propagating in a non-ambiguous inertial frame like EME2000 or GCRF

@MaximeJ

Thank you for this tip Maxime, I’ve implemented it that way and it worked. I tried to replicate the procedure with the DSST propagator using the JB2008 atmospheric model but i get the following error.
“Exception in thread “main” org.orekit.errors.OrekitException: altitude (-1,031,992.212 m) is below the 90,000 m allowed threshold”
I have an altitude detector implemented that’s meant to stop the propagation when the altitude gets below 120km but it doesn’t seem to work with the DSST propagator. My code is shown below.

    // Set path to orekit data folder
    File orekitData = new File("Path to Orekit Data folder");
    DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
    manager.addProvider(new DirectoryCrawler(orekitData));

    // TimeScale
    TimeScale timescale = TimeScalesFactory.getUTC();

    // Enter TLE
    String line1 = "1 23757U 95074A   12012.49666365  .00007180  00000-0  20301-3 0  9995";
    String line2 = "2 23757 022.9885 165.8605 0007544 261.0176 098.9451 15.32981482889705";
    TLE inputTLE = new TLE(line1, line2, timescale);
    System.out.println(inputTLE.getDate());

    // Retrieve Propagator used for TLE
    Propagator tlePropagator = TLEPropagator.selectExtrapolator(inputTLE);
    Orbit initialOrbit = new CartesianOrbit(tlePropagator.getInitialState().getOrbit());
    SpacecraftState initialState = new SpacecraftState(initialOrbit, mass);

    // Initial reference Orbit for the builder
    Orbit referenceOrbit = tlePropagator.getInitialState().getOrbit();
    System.out.println("Initial orbit");
    System.out.println(initialState.getDate());
    System.out.println(initialState.getOrbit());

    // Set up Integrator Parameters
    double minStep = 360;  // Minimum TimeStep
    double maxStep = 86400;  // Maximum TimeStep
    double dP = 10;  // Position Tolerance

    // Builder
    DormandPrince853IntegratorBuilder builder = new DormandPrince853IntegratorBuilder(minStep, maxStep, dP); // Propagator Factory
    DSSTPropagatorBuilder dsstPropbuilder = new DSSTPropagatorBuilder(referenceOrbit,builder,dP, PropagationType.MEAN,PropagationType.MEAN);
    dsstPropbuilder.setMass(mass);

    // Add Force models to Fitted Propagator Builder
    // Initialize Solar Activity Provider Data Provider
    final JB2008SpaceEnvironmentData jbdata = new JB2008SpaceEnvironmentData(JB2008SpaceEnvironmentData.DEFAULT_SUPPORTED_NAMES_SOLFSMY,JB2008SpaceEnvironmentData.DEFAULT_SUPPORTED_NAMES_DTC,manager,timescale);

    // Initialize JB2008 Atmospheric Model
    final JB2008 jb2008  = new JB2008(jbdata,sun, earth);

    dsstPropagatorBuilder.addForceModel(new DSSTAtmosphericDrag(jb2008, Cd, area, provider.getMu()));

    //  SGP4 Propagator Converter
    FiniteDifferencePropagatorConverter dsstFitter = new FiniteDifferencePropagatorConverter(dsstPropBuilder, 1.0e-3, 1000);

    // Fitting Parameters
    double duration = 86400.0; // Fitting Timespan
    double stepSize = 300.0;   // Fitting Stepsize
    Double points = duration / stepSize;  // Number of fitting points

    // Convert SGP4 Propagator to DSST Propagator
    System.out.println("START DSST PROPAGATOR FITTING");
    dsstFitter.convert(tlePropagator, duration, points.intValue()); // Conversion from SGP4 to DSST Propagator
    System.out.println("END DSST PROPAGATOR FITTING");
    DSSTPropagator dsstPropagator = (DSSTPropagator) dsstFitter.getAdaptedPropagator();  // Retrieving the converted numerical propagator
    System.out.println(dsstPropagator.getInitialState().getOrbit());

    // Converting Initial State to EME2000 Frame
    TimeStampedPVCoordinates tspvc1 = dsstPropagator.getInitialState().getPVCoordinates(FramesFactory.getEME2000());
    CartesianOrbit corbit = new CartesianOrbit(tspvc1,FramesFactory.getEME2000(),dsstPropagator.getMu());
    SpacecraftState newState = new SpacecraftState(corbit,mass);
    dsstPropagator.resetInitialState(newState);


    // Create Earth Model
    OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
            FramesFactory.getITRF(IERSConventions.IERS_2010, false));

    // Re-entry Detector
    EventDetector altdet = new AltitudeDetector(maxCheck, threshold, alt, earth).withHandler(((s, detector, decreasing) -> {
        System.out.println("Altitude is " + FastMath.min(((s.getA() * (1 + s.getE())) - Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 1000, ((s.getA() * (1 - s.getE())) - Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 1000) + " km and the spacecraft will re-enter soon");
        System.out.println("Re-entry event occurs on " + s.getDate().toString());
        return Action.STOP;
    }));


    // Add Event detector
    dsstPropagator.addEventDetector(altdet);

    // Setup Step Handlers
    final StatesHandler FittedStepHandler   = new StatesHandler();

    // Add fixed step handlers to the propagator
    dsstPropagator.getMultiplexer().add(3600., FittedStepHandler);

    // Get Final States
    final SpacecraftState finalState = dsstPropagator.propagate(initialState.getDate().shiftedBy(10. * Constants.JULIAN_YEAR));

    // Retrieve the states
    final List<SpacecraftState>  FittedStates = FittedStepHandler.getStates();

    final File output = new File("Output path" + name + ".txt");
    try (PrintStream stream = new PrintStream(output, StandardCharsets.UTF_8)) {
        stream.println("          Date                   SMA        Ecc        Incl        RAAN         AoP          MAN   ");
        stream.println("   [YYYY-MM-DD hh:mm:ss]         (km)                  (deg)       (deg)        (deg)        (deg)   ");

        for (SpacecraftState numstate: FittedStates){
            KeplerianOrbit Orb = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(numstate.getOrbit());
            String text = " %s   %.4f   %.7f   %.5f   %.5f   %.5f   %.5f";
            text = String.format(text,Orb.getDate(),Orb.getA()/1000,Orb.getE(),FastMath.toDegrees(Orb.getI()),FastMath.toDegrees(Orb.getRightAscensionOfAscendingNode()),FastMath.toDegrees(Orb.getPerigeeArgument()),FastMath.toDegrees(Orb.getMeanAnomaly()));
            stream.println(text);
        }
    }
    final String resultSavedAs = "Output file saved to file " + name + ".txt in specified path";
    System.out.println(resultSavedAs);

    KeplerianOrbit Neworbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit());
    System.out.println("Final Fitted orbit keplerian elements");
    System.out.println(Neworbit.toString());
    double years = (finalState.getDate().durationFrom((initialState.getDate()))) / (Constants.JULIAN_YEAR);
    System.out.println("The remaining Lifetime is " + years + "years");
    System.out.println("The remaining Lifetime with a 5% margin is " + years * 1.05 + "years");



}

private static class StatesHandler implements OrekitFixedStepHandler {

    // Points
    private final List<SpacecraftState> states;

    // Simple Constructor//
    StatesHandler() {
        // prepare an empty list of SpacecraftStates
        states = new ArrayList<SpacecraftState>();
    }

    public void handleStep(final SpacecraftState currentState) {
        // add the current state
        states.add(currentState);
    }

    /**
     * Get the list of handled orbital elements.
     *
     * @return orbital elements list
     */
    public List<SpacecraftState> getStates() {
        return states;
    }
}

Hi @dl00718,

I reckon getting such a negative altitude is weird, I’ll try to look into it.

I don’t have data for JB2008 so I’ll use NRLMSISE-00 instead, hoping to raise the same issue.

Could you please provide the data used for mass, Cd, area, the degree/order used for the gravity potential and the maxCheck and threshold used for the altitude detector ?

Hi Maxime,
Thanks for the timely response, the values are provided below.
Mass = 3007.93kg
Cd = 2.2
Area = 11.3162
Gravity potential degree and order = 6 , 6
I also had Sun and Moon third body perturbations
maxCheck = 60
threshold = 0.001

Thank you

Thank you for the clarifications !

Does the problem occur during the conversion or after during the 10 years propagation ?

If it occurs during the propagation could you simply share the converted initial state of the propagator (the conversion does not seem to work very well with MSIS00 here)

Hi Maxime,
The problem occured during propagation. I changed the maxCheck time to 7200 and that seems to fix it.
The initial state in the EME2000 frame is
SpacecraftState{orbit=Cartesian parameters: {P(-6638439.662248914, 1692150.1913326886, 8046.025780697497), V(-1733.267305451603, -6804.2916712986425, 2979.0833184996322)}, attitude=org.orekit.attitudes.Attitude@5fd4f8f5, mass=3007.93, additional={}, additionalDot={}}

By the way is my method for converting the initial state from TEME to EME2000 correct?
Thank you.

Hi,

Thanks for the info on the initial state.

Great, I suspected something related to a bad configuration of the altitude detector. Although I don’t see why increasing the maximum check interval would fix the problem. I’ll try to take a look when I have time.

It seems fine yes.

@MaximeJ
Actually I take that back. The error still occurs depending on the value of Position Tolerance i set when implementing the DormandPrinceIntegratorBuilder and the DSSTPropagator Builder. Using a value of 1e-3m works but switching to 1m or 10m as you used in an example on another thread, doesn’t work and leads to the negative altitude error.

Ok thanks for the update.

Quick question, we noted that you didn’t add any zonal or tesseral forces to your DSST propagator.
Is that true or does it just doesn’t appear on the snippet of code you sent ?

Could you send the last orbital elements written in your output file before the crash. So we don’t have to run the whole 10 years propagation to inspect the potential bug ?
And do you have a stack trace of the exception ?

My guess is that the integrator is trying steps that are too big and crashes when it evaluates the next state because the altitude is too low.
When reentering, if you’re flying downward at a few km/s it takes only tens of seconds to go from an altitude above 120km to an altitude below 90km.
So before the crash, the previous integration step has its previous and current state (start and end of the integration step) above 120km.
But when evaluating the next step, it tries to compute a state with an altitude below 90km and crashes because of the atmospheric model.
And that happens before the altitude detector is activated on this step.
That would explain why lowering the position tolerance solves the problem.

That being said, I’m not very familiar with the interconnection between integration steps and event detection, @luc and @evan.ward, do you have any ideas on this ?

I’ll take a look, do we have a failing test case?

I’ve another question. Does the error occur using a numerical propagator?

I agree that for 10 years propagation, the DSST is interesting. But I just would like to know if the error also occurs with the numerical method.

Bryan

I tried using the mass of code from above, but it was missing some variable definitions and other data.

I’ll take a look at the issue if you can make a JUnit test case that runs and fails. Preferably a minimal test case that only includes what is necessary to reproduce the issue and runs quickly. (That greatly helps with debugging!) I saw that you’re using JB2008. If that is necessary to reproduce the issue then you’ll also need to include those input files. Attaching a patch is a great way to include everything.

Back trace from the exception would also be useful.

I looked through the open issues, but nothing jumped out as being immediately related. Maxime’s theory sounds plausible, but I’m not sure why changing the maxCheck for the event detector would fix it. The step is computed before the event detector is called. Perhaps there is something extra going on with the DSST?