Hello,
I think I found two problems with the use of the “optimized” accelerationWrtState method in the Field version of HolmesFeatherstoneAttractionModel#acceleration.
- There is a bug in the check that is performed to decide whether to use the faster method or not.
isGradientStateDerivativedoes not check the partials of the date of the state, which however do impact the partials of the acceleration through the rotation matrix between state frame and ITRF- the gradient of
muis not checked. However,accelerationWrtStateuses only its real part.
- The check performed on the components of the position of the state is too restrictive. Indeed, it would be possible to call
accelerationWrtStatewith a position vector of the form[x+dx, y+dy, z+dz]as is it the case now, and then compose the output with the actual derivatives of the position.
To demonstrate (1), tests #2 and #3 fail:
@Test
void testHolmesFeatherstoneAttractionModelPartials() {
testHolmesFeatherstoneAttractionModelPartials(false, false);
}
@Test
void testHolmesFeatherstoneAttractionModelPartialsWrtEpoch() {
testHolmesFeatherstoneAttractionModelPartials(true, false);
}
@Test
void testHolmesFeatherstoneAttractionModelPartialsWrtMu() {
testHolmesFeatherstoneAttractionModelPartials(false, true);
}
void testHolmesFeatherstoneAttractionModelPartials(final boolean wrtDate, final boolean wrtMu) {
final File orekitData = new File(System.getProperty("user.home") + "/Resources/orekit-data");
final DataProvider dirCrawler = new DirectoryCrawler(orekitData);
DataContext.getDefault().getDataProvidersManager().addProvider(dirCrawler);
final HolmesFeatherstoneAttractionModel model = new HolmesFeatherstoneAttractionModel(
FramesFactory.getITRF(IERSConventions.IERS_2010, false),
GravityFieldFactory.getNormalizedProvider(8, 8)
);
final int prs = 8;
final DSFactory factory = new DSFactory(prs, 1);
final double mu = Constants.EGM96_EARTH_MU;
final Gradient gMu = wrtMu ? Gradient.variable(prs, 7, 0.0).multiply(1e6).add(mu) : Gradient.constant(prs, mu);
final DerivativeStructure dsMu = wrtMu ? factory.variable(7, 0.0).multiply(1e6).add(mu) : factory.constant(mu);
final double dt = 0.0;
final Gradient gDt = wrtDate ? Gradient.variable(prs, 6, dt) : Gradient.constant(prs, dt);
final DerivativeStructure dsDt = wrtDate ? factory.variable(6, dt) : factory.constant(dt);
final Frame frame = FramesFactory.getGCRF();
final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
final FieldAbsoluteDate<Gradient> gDate = new FieldAbsoluteDate<>(GradientField.getField(prs), date).shiftedBy(gDt);
final FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(factory.getDerivativeField(), date).shiftedBy(dsDt);
final Orbit orbit = new KeplerianOrbit(7e6, 0.2, 0.7, -1.0, 0.5, 0.7, PositionAngleType.TRUE, frame, date, mu);
final PVCoordinates pv = orbit.getPVCoordinates();
final FieldPVCoordinates<Gradient> gPv = new FieldPVCoordinates<>(
new FieldVector3D<>(
Gradient.variable(prs, 0, pv.getPosition().getX()),
Gradient.variable(prs, 1, pv.getPosition().getY()),
Gradient.variable(prs, 2, pv.getPosition().getZ())
),
new FieldVector3D<>(
Gradient.variable(prs, 3, pv.getVelocity().getX()),
Gradient.variable(prs, 4, pv.getVelocity().getY()),
Gradient.variable(prs, 5, pv.getVelocity().getZ())
)
);
final FieldSpacecraftState<Gradient> gState = new FieldSpacecraftState<>(new FieldCartesianOrbit<>(gPv, frame, gDate, gMu));
final FieldVector3D<Gradient> gAcc = model.acceleration(gState, new Gradient[]{gMu});
final FieldPVCoordinates<DerivativeStructure> dsPv = new FieldPVCoordinates<>(
new FieldVector3D<>(
factory.variable(0, pv.getPosition().getX()),
factory.variable(1, pv.getPosition().getY()),
factory.variable(2, pv.getPosition().getZ())
),
new FieldVector3D<>(
factory.variable(3, pv.getVelocity().getX()),
factory.variable(4, pv.getVelocity().getY()),
factory.variable(5, pv.getVelocity().getZ())
)
);
final FieldSpacecraftState<DerivativeStructure> dsState = new FieldSpacecraftState<>(new FieldCartesianOrbit<>(dsPv, frame, dsDate, dsMu));
final FieldVector3D<DerivativeStructure> dsAcc = model.acceleration(dsState, new DerivativeStructure[]{dsMu});
Assertions.assertArrayEquals(dsAcc.toVector3D().toArray(), gAcc.toVector3D().toArray(), Precision.EPSILON);
Assertions.assertArrayEquals(Arrays.copyOfRange(dsAcc.getX().getAllDerivatives(), 1, prs + 1), gAcc.getX().getGradient(), Precision.EPSILON);
Assertions.assertArrayEquals(Arrays.copyOfRange(dsAcc.getY().getAllDerivatives(), 1, prs + 1), gAcc.getY().getGradient(), Precision.EPSILON);
Assertions.assertArrayEquals(Arrays.copyOfRange(dsAcc.getZ().getAllDerivatives(), 1, prs + 1), gAcc.getZ().getGradient(), Precision.EPSILON);
}
while for (2), something like this works:
public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration2(final FieldSpacecraftState<T> s, final T[] parameters) {
final T mu = parameters[0];
// check for faster computation dedicated to derivatives with respect to state
if (isGradientStateDerivative2(s) && ((Gradient) mu).getAddendum().isZero()) {
@SuppressWarnings("unchecked")
final FieldVector3D<Gradient> p = (FieldVector3D<Gradient>) s.getPosition();
final FieldVector3D<Gradient> redPos = new FieldVector3D<>(
Gradient.variable(3, 0, p.getX().getValue()),
Gradient.variable(3, 0, p.getY().getValue()),
Gradient.variable(3, 0, p.getZ().getValue())
);
final FieldVector3D<Gradient> redAcc = accelerationWrtState(s.getDate().toAbsoluteDate(), s.getFrame(), redPos, Gradient.constant(3, mu.getReal()));
final RealMatrix jacAcc = new Array2DRowRealMatrix(new double[][]{
redAcc.getX().getGradient(),
redAcc.getY().getGradient(),
redAcc.getZ().getGradient()
});
final RealMatrix jacPos = new Array2DRowRealMatrix(new double[][]{
p.getX().getGradient(),
p.getY().getGradient(),
p.getZ().getGradient()
});
final RealMatrix jac = jacAcc.multiply(jacPos);
@SuppressWarnings("unchecked")
final FieldVector3D<T> a = new FieldVector3D<>(
(T) new Gradient(redAcc.getX().getValue(), jac.getRow(0)),
(T) new Gradient(redAcc.getY().getValue(), jac.getRow(1)),
(T) new Gradient(redAcc.getZ().getValue(), jac.getRow(2))
);
return a;
}
// get the position in body frame
final FieldAbsoluteDate<T> date = s.getDate();
final FieldStaticTransform<T> fromBodyFrame = bodyFrame.getStaticTransformTo(s.getFrame(), date);
final FieldStaticTransform<T> toBodyFrame = fromBodyFrame.getInverse();
final FieldVector3D<T> position = toBodyFrame.transformPosition(s.getPosition());
// gradient of the non-central part of the gravity field
return fromBodyFrame.transformVector(new FieldVector3D<>(gradient(date, position, mu)));
}
private <T extends CalculusFieldElement<T>> boolean isGradientStateDerivative2(final FieldSpacecraftState<T> state) {
return (state.getPosition().getX() instanceof Gradient) && state.getDate().getDate().hasZeroField();
}
Cheers,
Alberto