Hi @luc,
I found that Orekit’s InertialForces formula is correct, but the transform it receives is incomplete. In Orekit 13.1, the scalar MODProvider.getTransform(AbsoluteDate) returns the time-dependent precession rotation without its angular rate and angular acceleration. TOD is built through MOD, so the initial EME2000→TOD PVA transformation and every subsequent InertialForces evaluation miss the same derivatives. The missing initial GEO velocity contribution is about 3.258e-4 m/s.
Correcting only this initial PVA is not sufficient because the derivatives are also required throughout propagation.
The essential Orekit-side change is small (imports omitted):
- final double[] angles = precessionFunction.value(date);
- final Rotation precession = r4.compose(
- new Rotation(RotationOrder.ZXZ, RotationConvention.FRAME_TRANSFORM,
- -angles[0], -angles[1], angles[2]),
- RotationConvention.FRAME_TRANSFORM);
- return new Transform(date, precession);
+ final FieldAbsoluteDate<UnivariateDerivative2> derivativeDate =
+ new FieldAbsoluteDate<>(UnivariateDerivative2Field.getInstance(), date)
+ .shiftedBy(new UnivariateDerivative2(0.0, 1.0, 0.0));
+ final UnivariateDerivative2[] angles = precessionFunction.value(derivativeDate);
+ final FieldRotation<UnivariateDerivative2> derivativeR4 =
+ new FieldRotation<>(UnivariateDerivative2Field.getInstance(), r4);
+ final FieldRotation<UnivariateDerivative2> precession = derivativeR4.compose(
+ new FieldRotation<>(RotationOrder.ZXZ, RotationConvention.FRAME_TRANSFORM,
+ angles[0].negate(), angles[1].negate(), angles[2]),
+ RotationConvention.FRAME_TRANSFORM);
+ return new Transform(date, new AngularCoordinates(precession));
My application-side CompleteInertialForces applies this derivative recovery
only to MOD precession, composes Orekit’s normal MOD→TOD transform, and otherwise
uses the same Coriolis, centrifugal, Euler and translational formulas as
InertialForces.
Results with Cartesian propagation and 36Ă—36 gravity:
-
GEO, seven-day fixed-step RK4 (10 s): 196.759 m with Orekit versus 0.0038 m corrected.
-
LEO (a=7000 km, e=0.001, i=51.6 deg), seven days: 29.971 m with Orekit versus 0.0035 m corrected.
-
Correcting only the initial PVA gives 112.753 m in the LEO case.
The optimized workaround used essentially the same integration steps/evaluations and was about 8–9% slower in the LEO timing. Propagating in EME2000 and transforming only the output remains the no-workaround-overhead option. No third-body forces were included.
I have also implemented a more general method to retrieve the full transform, you can find it at the end of this post. The drawbacks is that the propagation takes 1.6x time as the original one.
NOTE: this problem is also highlighted in the classic .getPVCoordinates(Frame) method, as expected I might add since it use the same Frame transform methods that don’t include the all the terms needed.
Here the graph with the overall results for a GEO orbit, I have also tested a LEO satellite and the results are good as well.
Here my more general implementation, the completeTransform is where the “magic” happens:
@Override
public Vector3D acceleration(final SpacecraftState state, final double[] parameters) {
// During propagation these are the coordinate (fictitious) terms needed
// when the equations are integrated in state.getFrame().
final Transform transform = forceTransform(state.getFrame(), state.getDate());
final Vector3D omega = transform.getRotationRate();
final Vector3D position = state.getPosition();
final Vector3D velocity = state.getPVCoordinates().getVelocity();
return transform.getRotation().applyTo(transform.getAcceleration())
.add(Vector3D.crossProduct(omega, velocity).scalarMultiply(-2.0))
.subtract(Vector3D.crossProduct(omega, Vector3D.crossProduct(omega, position)))
.subtract(Vector3D.crossProduct(transform.getRotationAcceleration(), position));
}
@Override
public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(
final FieldSpacecraftState<T> state, final T[] parameters) {
final FieldTransform<FieldUnivariateDerivative2<T>> transform =
referenceInertialFrame.getTransformTo(state.getFrame(), state.getDate().toFUD2Field());
final FieldPVCoordinates<T> cartesian = new FieldPVCoordinates<>(transform.getCartesian().getPosition());
final FieldAngularCoordinates<T> angular = new FieldAngularCoordinates<>(transform.getRotation());
final FieldVector3D<T> omega = angular.getRotationRate();
final FieldVector3D<T> position = state.getPosition();
final FieldVector3D<T> velocity = state.getPVCoordinates().getVelocity();
return angular.getRotation().applyTo(cartesian.getAcceleration())
.add(FieldVector3D.crossProduct(omega, velocity).scalarMultiply(-2.0))
.subtract(FieldVector3D.crossProduct(omega, FieldVector3D.crossProduct(omega, position)))
.subtract(FieldVector3D.crossProduct(angular.getRotationAcceleration(), position));
}
@Override
public List<ParameterDriver> getParametersDrivers() {
return Collections.emptyList();
}
private Transform completeTransform(final Frame targetFrame, final AbsoluteDate date) {
final FieldAbsoluteDate<UnivariateDerivative2> derivativeDate =
new FieldAbsoluteDate<>(UnivariateDerivative2Field.getInstance(), date)
.shiftedBy(new UnivariateDerivative2(0.0, 1.0, 0.0));
final FieldTransform<UnivariateDerivative2> transform =
referenceInertialFrame.getTransformTo(targetFrame, derivativeDate);
return new Transform(date,
new PVCoordinates(transform.getCartesian().getPosition()),
new AngularCoordinates(transform.getRotation()));
}
In conclusion, I wonder if the best path forward is to open an issue on GitLab. @luc , what do you think ?
Thanks,
@DDega
