Issues regarding numerical propagator

razvicosti · February 15, 2023, 7:06pm

Hello Orekit team,

I am new to Orekit and I am trying to write a code that would propagate backwards the orbit of a satellite which is supposed to be performing an EOR and then comparing the results to the ones obtained from a simple backwards propagated orbit using a keplerian propagator to validate the concept.

The code I wrote for the backwards propagation of the maneuver is the following:

 public static DragForce dragForce(double crossSection, double dragCoef){
        Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010,true);
        IsotropicDrag drag = new IsotropicDrag(crossSection,dragCoef);
        OneAxisEllipsoid earthShape = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,Constants.WGS84_EARTH_FLATTENING,frame);
        SimpleExponentialAtmosphere earthAtmosphere = new SimpleExponentialAtmosphere(earthShape, 0.0949,20000.0,7500);
        return new DragForce(earthAtmosphere,drag);

    }

    public static void main(String[] args) throws IOException {
        File orekitData = new File("/home/razvan/Downloads/orekit-data");
        DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
        manager.addProvider(new DirectoryCrawler(orekitData));
        Frame inertialFrame = FramesFactory.getEME2000();
        TimeScale utc = TimeScalesFactory.getUTC();
        AbsoluteDate initialDate = new AbsoluteDate(2023, 01, 31, 17, 46, 49.924704, utc);
        double mu = 3.986004415e+14;
        double a = 504848; //semi-major axis in metres
        double e = 0.00101990;
        double i = FastMath.toRadians(98.6875);
        double omega = FastMath.toRadians(131.2378);
        double raan = FastMath.toRadians(39.9845);
        double lM = FastMath.toRadians(337.9425);
        double mass = 6000;

        double crossSection = 10;
        double dragCoef = 2.2;

        Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, inertialFrame, initialDate, mu);

        double dP = 1;
        double minStep = 0.001;
        double maxStep = 1000;
        double initStep = 20;
        double[][] tolerance = NumericalPropagator.tolerances(dP,initialOrbit, OrbitType.KEPLERIAN);
        AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]);
        integrator.setInitialStepSize(initStep);
        NumericalPropagator propagator = new NumericalPropagator(integrator);

        propagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.TNW));
        double duration = 2*24*3600;
        ManeuverTriggers triggers = new DateBasedManeuverTriggers(firingDate, duration);
        double thrust = 250*FastMath.pow(10,-3);
        double isp = 2000;
        propagator.addForceModel(new NewtonianAttraction(mu));
        propagator.addForceModel(dragForce(crossSection,dragCoef));
        propagator.addForceModel(new 
ConfigurableLowThrustManeuver(ThrustDirectionAndAttitudeProvider.buildFromFixedDirectionInSatelliteFrame(Vector3D.PLUS_I),triggers,thrust,isp));
        File csvFile = new File("bck_propagation_data.csv");
        FileWriter fileWriter = new FileWriter(csvFile);
        PrintWriter printWriter = new PrintWriter(fileWriter, true);
        String header = "time,pos_X,pos_Y,pos_Z,inclination";
        printWriter.println(header);


       propagator.getMultiplexer().add(20, (state) -> { 
                   LocalOrbitalFrame localFrame = new LocalOrbitalFrame(inertialFrame,LOFType.TNW, state.getOrbit(), "TNWframe");
                   Transform inertialToLocal = inertialFrame.getTransformTo(localFrame,state.getDate());
                    printWriter.printf(Locale.US, "%s, %.3f,  %.3f,  %.3f, %.3f %n",
                            state.getDate(), inertialToLocal.transformVector(state.getPVCoordinates().getPosition()).getX(), inertialToLocal.transformVector(state.getPVCoordinates().getPosition()).getY(),
                            inertialToLocal.transformPVCoordinates(state.getPVCoordinates()).getPosition().getZ(), FastMath.toDegrees(state.getI()));
}
);
        propagator.propagate(initialDate,firingDate );


        System.exit(0);

This is a code that I wrote reading the forum and the tutorials of Orekit and adapting the answers to my needs. However, I have the following questions:

What is the cause for which I receive the following error when running this code and how could I solve the issue:

 Exception in thread "main" org.orekit.errors.OrekitException: NaN appears during integration near time -20
	at org.orekit.errors.OrekitException.unwrap(OrekitException.java:154)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:511)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:425)
	at org.example.Main.main(Main.java:128)
Caused by: org.hipparchus.exception.MathIllegalStateException: NaN appears during integration near time -20
	at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:269)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:477)
	... 2 more

To me it seems a little bit strange that when converting to RTN the value that I obtain when calling .getY() is around the value of the semimajor axis. Is that normal or do I make a mistake when computing RTN? The results I am talking about were obtained when performing the backward propagation without the aerodynamic drag force model added. To me it was strange as well that the Z coordinate was always 0. That I thought to be normal considering that I configured the low thrust maneuver to be in along-track direction, but am I right?

Thank you very much,
Razvan

Serrof · February 15, 2023, 8:30pm

Hi there,

You forgot to add the earth radius to your semi major axis no? This would explain why the atmospheric model does not like to be evaluated
Also, do not add the Keplerian part of the dynamics to the propagator as it is already included by default.

Best,
Romain.

razvicosti · February 22, 2023, 7:56am

Thank you very much for your answer, Romain. Indeed I didn’t add the Earth’s radius to my semimajor axis. Could you tell me as well if it is normal what I asked in question number 2?

All the best,
Razvan

Serrof · February 22, 2023, 8:09am

Hi Razvan,

In your code you’re using the TNW frame, not the RTN one. So the first axis is colinear to the velocity vector. And the second one is near radial for quasi circular orbits. See orekit documentation for more details.

Best,
Romain.

PS: remember that you do not need to add the NewtonianAttraction as it is already there by default

razvicosti · February 27, 2023, 9:12pm

Hello once again Romain,

I would have another question that has been haunting me. So I did as you said, thank you very much for your answer. However, I am note sure that I am defining the along-track maneuver as it should. In the meantime, I have changed the local frame to QSW which is the same as RTN according to CCSDS. I read as well other posts on the forum concerning the topic of maneuvers and I found out that the attitude provider is responsible for orienting the thrust in the specific direction in the local frame. Then, I also configure the thrust direction in the configurable low thrust maneuver. My question would be am I doing it right? Thank you very much once again for all your kind answers.

Best,
Razvan

Serrof · February 28, 2023, 6:11pm

Hello Razvan,

I think your understanding is correct. Replacing TNW by QSW (which is indeed the same as RTN) in your previous code, that should have you thrust along Q (frame vectors are generally called I, J and K) so in the direction of the position vector.

Best,
Romain.