Hey, folks!
My 2-month orbit propagation (60s steps) runs in ~30min with gravity alone but jumps to 103min when adding atmospheric drag. Including other forces (SRP, third-body), I suppose it gets worse.
I ran alone the forces and got the following times:
- Just propagation: 12.46s
- Gravitational Forces: 1876s
- Atmospheric Drag: 66.65s
- SRP : 48.45s
- Third Bodies (Sun & Moon): 16.01s
How can I reduce computation time without sacrificing critical accuracy?
public class SpacecraftDecay {
public static void main (String[] args) throws Exception {
long startTime = System.currentTimeMillis();
OrekitInitializer.loadFiles();
DataContext context = DataContext.getDefault();
TimeScale tai = context.getTimeScales().getTAI();
// Epoch start and end
AbsoluteDate startEpoch = new AbsoluteDate(2025, 1, 21, 0, 0, 37.0, tai);
AbsoluteDate endEpoch = new AbsoluteDate(2025, 3, 21, 0, 0, 37.0, tai);
double totalDuration = endEpoch.durationFrom(startEpoch); // Total duration in seconds
// Initial orbital parameters
final double a = 6918.601 * 1000; // Semi Major Axis in metros
final double e = 0.001; // eccentricity
final double i = FastMath.toRadians(97.556); // inclination
final double raan = FastMath.toRadians(-102.078);
final double w = FastMath.toRadians(90.0); // arg of Periapsis
final double v = FastMath.toRadians(0.0); // True anomaly
// Frames
Frame inertialFrame = FramesFactory.getGCRF();
Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
// Earth Body Shape
OneAxisEllipsoid earthShape = new OneAxisEllipsoid(
Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
earthFrame
);
// gravitation coefficient
double mu = Constants.WGS84_EARTH_MU;
// Orbit construction as Keplerian
Orbit initialOrbit = new KeplerianOrbit(a, e, i, w, raan, v, PositionAngleType.TRUE, inertialFrame, startEpoch, mu);
// Attitude provider (Nadir pointing)
AttitudeProvider attitudeProvider = new NadirPointing(inertialFrame, earthShape);
// Create the initial state of the spacecraft
double mass = 280;
SpacecraftState initialState = new SpacecraftState(initialOrbit, attitudeProvider.getAttitude(initialOrbit, startEpoch, inertialFrame), mass);
// Numerical integrator configuration with tolerances
double minStep = 1e-11; // minimum step in s //GODOT_DEFAULT:
double maxStep = 300.0; // maximum step in s
double positionTolerance = 1e-11; // Position tolerance
// Tolerance estimate based on initial orbit
double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, OrbitType.KEPLERIAN);
AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(OrbitType.KEPLERIAN);
propagator.setAttitudeProvider(attitudeProvider);
propagator.setInitialState(initialState);
/*
// Perturbations
// 1. Gravitational Perturbations (Spherical Harmonics JGM-3 4x4)
NormalizedSphericalHarmonicsProvider gravityProvider = GravityFieldFactory.getNormalizedProvider(4,4);
ForceModel earthGravity = new HolmesFeatherstoneAttractionModel(earthFrame, gravityProvider);
propagator.addForceModel(earthGravity);
// 2. Atmospheric Drag
MarshallSolarActivityFutureEstimation solarActivity =
new MarshallSolarActivityFutureEstimation(
MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE
);
double dragArea = 3.0;
double cd = 2.2;
Atmosphere atmosphere = new NRLMSISE00(solarActivity, CelestialBodyFactory.getSun(), earthShape);
IsotropicDrag isotropicDrag = new IsotropicDrag(dragArea, cd);
ForceModel drag = new DragForce(atmosphere, isotropicDrag);
propagator.addForceModel(drag);
// 3. Solar Radiation Pressure
double srpArea = 4.0;
double cr = 1.2;
RadiationSensitive spacecraftSRP = new IsotropicRadiationSingleCoefficient(srpArea, cr);
SolarRadiationPressure srp = new SolarRadiationPressure(CelestialBodyFactory.getSun(), earthShape, spacecraftSRP);
propagator.addForceModel(srp);
// 4. Third bodies: Sun and Moon
propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
*/
// Stop criterion when semi-axis less than or equal to 520 km above the surface (~6891 km SMA)
double earthRadius = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
double targetSMA = earthRadius + 520000.0; // in metres
// Track last printed percentage to avoid flooding the console
final int[] lastPercent = { -1 };
// Step handler to track orbit evolution and stop at specific semi-major axis
propagator.setStepHandler(60.0, spacecraftState -> {
final KeplerianOrbit o = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(spacecraftState.getOrbit());
//System.out.println(o.getA());
/*System.out.format(Locale.US, "%s %12.11f %10.8f %10.6f %10.6f %10.6f %10.6f%n",
spacecraftState.getDate(),
o.getA(), o.getE(),
FastMath.toDegrees(o.getI()),
FastMath.toDegrees(o.getPerigeeArgument()),
FastMath.toDegrees(o.getRightAscensionOfAscendingNode()),
FastMath.toDegrees(o.getTrueAnomaly()));
*/
// Calculate progress
double elapsed = spacecraftState.getDate().durationFrom(startEpoch);
double progressPercent = (elapsed / totalDuration) * 100.0;
int currentPercent = (int) progressPercent;
// Update progress bar only when percentage increases
if (currentPercent > lastPercent[0]) {
lastPercent[0] = currentPercent;
int barWidth = 50; // Width of the progress bar
int filled = (int) (barWidth * (currentPercent / 100.0));
StringBuilder bar = new StringBuilder("[");
// Inside the lambda's progress bar loop:
for (int j = 0; j < barWidth; j++) {
if (j < filled) bar.append("=");
else if (j == filled) bar.append(">");
else bar.append(" ");
}
bar.append("] ").append(currentPercent).append("%");
System.out.print("\r" + bar.toString());
System.out.flush();
}
// Check condition
if (o.getA() <= targetSMA) {
throw new RuntimeException("Target semi-major axis reached");
}
});
// Propagate until event or end date
SpacecraftState finalState = propagator.propagate(endEpoch);
/*
At the end, finalState will be the state when the SMA reaches ~520 km,
or the end date if it doesn't reach it before.
*/
System.out.println("Final date: " + finalState.getDate());
System.out.println("Final SMA (km): " + (finalState.getA() / 1000.0));
// Add runtime calculation at the end
long endTime = System.currentTimeMillis();
double durationSeconds = (endTime - startTime) / 1000.0;
System.out.printf(Locale.US, "\nTotal runtime: %.2f seconds%n", durationSeconds);
}
}