Below is the code I think should reproduce the error (apologies for how long it is). I believe it comes down to a timing error. When investigating, the scheduled burn is set for a time before the next time step.
public class simpleAtmosphereDrag {
// ----- Simulation Constants
private static final String START_TIME = "2025-12-29T00:00:00";
private static final double SIMULATION_TIME_S = 4.9 * 3600.;
private static final double ALTITUDE_START_KM = 175, ALTITUDE_END_KM = 175, ALTITUDE_DELTA = 25;
// ----- Spacecraft detail
private final static double MASS_KG = 100.0; // kg
private final static double CROSS_SECTION_M2 = 0.5;
private final static double CD = 2.5; // for constant Cd model,
// ----- Thruster
private final static double EFFICIENCY = 0.5;
private final static double THRUST = 0.6; // N
private final static double ISP = 3000.0; // s
private final static double DALTITUDE_THRUST = 5; // How far sc falls below starting alt before thrusting (km)
// ----- Propagator details
private static final double MIN_STEP = 1.0;
private static final double MAX_STEP = 100.0;
private static final double POSITION_TOLERANCE = 10.0;
private simpleAtmosphereDrag() {}
public static void main(String[] args) {
// ====== Initialize Orekit data (needs a folder with Earth orientation, ephemeris, etc.)
initOrekit.initializeOrekitData();
for (double ALTITUDE_KM = ALTITUDE_START_KM; ALTITUDE_KM <= ALTITUDE_END_KM; ALTITUDE_KM += ALTITUDE_DELTA){
System.out.printf("Alt: %.0f\n", ALTITUDE_KM);
// ====== Setup Reference Frames and Time Parameters
Frame emeFrame = FramesFactory.getEME2000();
Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true); // Earth Frame
AbsoluteDate startDate = new AbsoluteDate(START_TIME, TimeScalesFactory.getUTC());
AbsoluteDate endDate = startDate.shiftedBy(SIMULATION_TIME_S);
// ====== Generate the Earth and Sun
BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
earthFrame);
PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
// ====== Arbitrary Orbit (Keplerian)
Orbit orbit = buildOrbit(ALTITUDE_KM, startDate);
// ====== Atmosphere and Drag Models
NRLMSISE00InputParameters weatherData = new CssiSpaceWeatherData("SpaceWeather-All.*\\.txt");
Atmosphere atmosphere = new NRLMSISE00(weatherData, sun, earth);
IsotropicDrag spacecraft = new IsotropicDrag(CROSS_SECTION_M2, CD);
DragForce dragModel = new DragForce(atmosphere, spacecraft);
// ====== Gravity Force Model
ForceModel gravityModel = createGravityModel(earthFrame);
// ====== Attitude Provider (Nadir Locking for now)
AttitudeProvider attitudeProvider = new LofOffset(emeFrame, LOFType.TNW);
// ====== Thruster Model
ThrusterModel thrusterModel = new ThrusterModel(EFFICIENCY, THRUST, ISP);
// ====== Propagator
NumericalPropagator propagator = createNumericalPropagator(orbit, POSITION_TOLERANCE, MIN_STEP, MAX_STEP);
// ====== Propulsion Detectors
double thresholdSMAm = (ALTITUDE_KM - DALTITUDE_THRUST) * 1e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS; // m
double targetSMAm = ALTITUDE_KM * 1e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS; // m
// Create shared guidance state
GuidanceState guidanceState = new GuidanceState();
// ---- SMA Detector - when SC falls below threshold
SMAHandler smaHandler = new SMAHandler(guidanceState);
SMADetector smaDetector = new SMADetector(smaHandler, thresholdSMAm);
// ---- Perigee Detector - when at perigee to perform orbit raising burn
PerigeeHandlerFull perigeeHandlerFull = new PerigeeHandlerFull(guidanceState, propagator, thrusterModel, targetSMAm);
PositionAngleDetector perigeeDetector = new PositionAngleDetector(OrbitType.KEPLERIAN, PositionAngle.TRUE, 0.).withHandler(perigeeHandlerFull);
// ====== Add everything to the Propagator
// ---- 1. Set Initial Spacecraft State
SpacecraftState initialState = new SpacecraftState(orbit, MASS_KG);
propagator.setInitialState(initialState);
propagator.setAttitudeProvider(attitudeProvider);
// ---- 2. Add all force models
propagator.addForceModel(dragModel);
propagator.addForceModel(gravityModel);
// ---- 3. Add propulsion detectors
propagator.addEventDetector(smaDetector);
propagator.addEventDetector(perigeeDetector);
// ---- 4. Add step handlers/output
StepHandlerMultiplexer multiplexer = new StepHandlerMultiplexer();
// Debugging when forces are added
OrekitStepHandler forcesPrinter = PrintPropagatorForces(propagator);
multiplexer.add(forcesPrinter);
propagator.setStepHandler(multiplexer);
// ====== Propagate
SpacecraftState currentState = propagator.propagate(endDate);
System.out.printf("Simulated to: %s", currentState.getDate());
}
}
// ================================ Helper Methods ================================
private static Orbit buildOrbit(double altitudeKm, AbsoluteDate startDate) {
double ra = altitudeKm * 1000.0; // altitude at apogee
double rp = ra; // altitude at perigee
double ae = Constants.WGS84_EARTH_EQUATORIAL_RADIUS; // Earth's radius
double a = (rp + ra + 2* ae) / 2.0; // Semi-major axis in meters
double e = 0; // circular orbit
double i = FastMath.toRadians(98.);
double mu = Constants.WGS84_EARTH_MU;
Frame emeFrame = FramesFactory.getEME2000();
Orbit kepOrbit = new KeplerianOrbit(a, e, i, 0.0, 0.0, 0.0, PositionAngle.TRUE, emeFrame, startDate, mu);
return kepOrbit;
}
private static ForceModel createGravityModel(Frame earthFrame) {
NormalizedSphericalHarmonicsProvider gravityField =
GravityFieldFactory.getNormalizedProvider(0,0);
ForceModel gravityModel = new HolmesFeatherstoneAttractionModel(earthFrame, gravityField);
return gravityModel;
}
private static NumericalPropagator createNumericalPropagator(Orbit initialOrbit, double positionTolerance, double minStep, double maxStep) {
double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, OrbitType.CARTESIAN);
AdaptiveStepsizeIntegrator integrator = new DormandPrince54Integrator(
minStep,
maxStep,
tolerances[0], // absolute tolerances (array)
tolerances[1] // relative tolerances (array)
);
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(OrbitType.CARTESIAN);
return propagator;
}
private static OrekitStepHandler PrintPropagatorForces(
NumericalPropagator propagator) {
return new OrekitStepHandler() {
@Override
public void handleStep(OrekitStepInterpolator interpolator) {
SpacecraftState s = interpolator.getCurrentState();
List<ForceModel> forces = propagator.getAllForceModels();
// Only print if more than 3 forces active
if (forces.size() > 3) {
// Check all ForceModels
double a = s.getA();
double altKm = (a - Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 1e3; // km
System.out.printf("t=%s, SMA Altitude = %.2f km%n", s.getDate(), altKm);
for (ForceModel fm : propagator.getAllForceModels()) {
double[] params = fm.getParameters(); // Get the parameters to receive accelerations
Vector3D acc = fm.acceleration(s, params);
System.out.printf(" ForceModel active: %s, acc=%.2e m/s2%n",
fm.getClass().getSimpleName(),
acc.getNorm()
);
}
System.out.println();
}
}
@Override
public void finish(SpacecraftState finalState) {}
};
}
}
public class GuidanceState {
public boolean belowThreshold = false;
public BurnPhase phase = BurnPhase.IDLE;
}
public class PerigeeHandlerFull implements EventHandler<PositionAngleDetector>{
private final GuidanceState guidanceState;
private final NumericalPropagator propagator;
private final ThrusterModel thruster;
private final double desiredRad;
public PerigeeHandlerFull(GuidanceState guidanceState,
NumericalPropagator propagator,
ThrusterModel thruster,
double desiredRad) {
this.guidanceState = guidanceState;
this.propagator = propagator;
this.thruster = thruster;
this.desiredRad = desiredRad;
}
@Override
public Action eventOccurred(SpacecraftState s,
PositionAngleDetector detector,
boolean increasing) {
if (!increasing) {
return Action.CONTINUE;
}
if (guidanceState.belowThreshold &&
guidanceState.phase == BurnPhase.WAIT_PERIGEE) {
// Constants
double mu = Constants.WGS84_EARTH_MU;
// State Variables
PVCoordinates perigeePV = s.getPVCoordinates();
double perigeePosNorm = perigeePV.getPosition().getNorm();
Vector3D perigeeVel = perigeePV.getVelocity();
double perigeeVelNorm = perigeeVel.getNorm();
// Thruster Quantities
double perigeeMass = s.getMass();
double thrust = thruster.getThrust();
double isp = thruster.getISP();
// Burn values - assuming starting from ~circular orbit
double aTrans = (perigeePosNorm + desiredRad) / 2;
double vTrans = FastMath.sqrt(2 * mu / perigeePosNorm - mu / aTrans);
double dvPerigee = vTrans - perigeeVelNorm;
// Calcualte the burn time - bad approximation, will need to do something else eventually
AbsoluteDate startPerigeeBurn = s.getDate();
double perigeeBurnTime = dvPerigee * perigeeMass / thrust;
Vector3D burnDir = perigeeVel.normalize();
ConstantThrustManeuver perigeeBurn = new ConstantThrustManeuver(startPerigeeBurn, perigeeBurnTime, thrust, isp, burnDir);
propagator.addForceModel(perigeeBurn);
guidanceState.phase = BurnPhase.WAIT_APOGEE;
System.out.println(" --- Made it to Perigee, Scheduled Burn");
System.out.printf(" --- ConstantThrustManuever(%s, %.2f, %.2f, %.1f, [%.2f, %.2f, %.2f])%n", startPerigeeBurn, perigeeBurnTime, thrust, isp, burnDir.getX(), burnDir.getY(), burnDir.getZ());
System.out.println();
guidanceState.phase = BurnPhase.WAIT_APOGEE;
}
return Action.CONTINUE;
}
}
public class SMAHandler implements EventHandler<SMADetector>{
private final GuidanceState guidanceState;
public SMAHandler(GuidanceState guidanceState){
this.guidanceState = guidanceState;
}
@Override
public Action eventOccurred(SpacecraftState s,
SMADetector detector,
boolean increasing) {
if (!increasing) { // crossed downward -> below threshold
guidanceState.belowThreshold = true;
guidanceState.phase = BurnPhase.WAIT_PERIGEE;
// Debugging print infor
double a = s.getA();
double altKm = (a - Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 1e3; // km
System.out.printf("Spacecraft at %.1f km -> Waiting for Perigee%n", altKm);
}
return Action.CONTINUE;
}
}
public class SMADetector extends AbstractDetector<SMADetector> {
private final double thresholdSMAm; // the threshold semi-major axis (m), altitude + equitorial radius
public SMADetector(EventHandler<SMADetector> handler, double thresholdSMAm) {
super(DEFAULT_MAXCHECK, DEFAULT_THRESHOLD, DEFAULT_MAX_ITER, handler);
this.thresholdSMAm = thresholdSMAm;
}
public SMADetector(double maxcheck, double threshold, int max_iter, EventHandler<SMADetector> handler, double thresholdAltitudeKm) {
super(maxcheck, threshold, max_iter, handler);
this.thresholdSMAm = thresholdAltitudeKm;
}
public double getThresholdAltitudeM() {
return thresholdSMAm;
}
@Override
public double g(SpacecraftState s) {
double aM = s.getA();
// g positive -> negative when SMA is less than lower threshold, thrust ON
// g negative -> positive when SMA is less than higher threshold, thrust OFF
double g = aM - thresholdSMAm;
return g;
}
@Override
protected SMADetector create(double newMaxCheck, double newThreshold, int newMaxIter, EventHandler newHandler) {
return new SMADetector(newMaxCheck, newThreshold, newMaxIter, newHandler, thresholdSMAm);
}
}
public class ThrusterModel {
/** Input values */
private final double efficiency;
private final double thrust; // Thrust (N)
private final double isp; // Specific impulse (s) - currently for tracking fuel consumption
/** Calculated values */
private final double power; // Power (kW) - f(thrust, isp, efficiency)
private final double mdot; // Exhaust mass flow rate - f(thrust, isp)
private final double exhaustVelocity; // Thrust exhaust velocity (km/s) - f(thrust, mdot)
public ThrusterModel(double efficiency, double thrust, double isp) {
this.efficiency = efficiency;
this.thrust = thrust;
this.isp = isp;
double g0 = Constants.G0_STANDARD_GRAVITY;
this.power = this.thrust * g0 * this.isp / (2 * this.efficiency) / 1000.; // Feldman2022
this.mdot = this.thrust / (this.isp * g0);
this.exhaustVelocity = this.thrust / this.mdot / 1000.;
}
public double getThrust() {
return this.thrust;
}
public double getMassFlowRate() {
return this.mdot;
}
public double getPower() {
return this.power;
}
public double getExhaustVelocity() {
return this.exhaustVelocity;
}
public double getISP() {
return this.isp;
}
}
The output when I run this is:
Alt: 175
Spacecraft at 170.0 km -> Waiting for Perigee
--- Made it to Perigee, Scheduled Burn
--- ConstantThrustManuever(2025-12-29T04:51:50.00006656309779Z, 303.06, 0.60, 3000.0, [-0.91, 0.06, -0.42])
t=2025-12-29T04:51:51.000Z, SMA Altitude = 168.92 km
ForceModel active: DragForce, acc=1.53e-04 m/s2
ForceModel active: HolmesFeatherstoneAttractionModel, acc=0.00e+00 m/s2
ForceModel active: ConstantThrustManeuver, acc=0.00e+00 m/s2
ForceModel active: NewtonianAttraction, acc=9.30e+00 m/s2
t=2025-12-29T04:53:31.000Z, SMA Altitude = 168.89 km
ForceModel active: DragForce, acc=1.69e-04 m/s2
ForceModel active: HolmesFeatherstoneAttractionModel, acc=0.00e+00 m/s2
ForceModel active: ConstantThrustManeuver, acc=0.00e+00 m/s2
ForceModel active: NewtonianAttraction, acc=9.30e+00 m/s2
t=2025-12-29T04:54:00.000Z, SMA Altitude = 168.88 km
ForceModel active: DragForce, acc=1.75e-04 m/s2
ForceModel active: HolmesFeatherstoneAttractionModel, acc=0.00e+00 m/s2
ForceModel active: ConstantThrustManeuver, acc=0.00e+00 m/s2
ForceModel active: NewtonianAttraction, acc=9.30e+00 m/s2
. I’ll try to look into it asap.