use of org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction in project Orekit by CS-SI.
the class PoissonSeriesParserTest method testDerivativesFromFieldAPI.
@Test
public void testDerivativesFromFieldAPI() throws OrekitException {
Utils.setDataRoot("regular-data");
String directory = "/assets/org/orekit/IERS-conventions/";
PoissonSeriesParser parser = new PoissonSeriesParser(17).withPolynomialPart('t', PolynomialParser.Unit.NO_UNITS).withFirstDelaunay(4).withFirstPlanetary(9).withSinCos(0, 2, 1.0, 3, 1.0);
InputStream xStream = getClass().getResourceAsStream(directory + "2010/tab5.2a.txt");
PoissonSeries xSeries = parser.parse(xStream, "2010/tab5.2a.txt");
InputStream yStream = getClass().getResourceAsStream(directory + "2010/tab5.2b.txt");
PoissonSeries ySeries = parser.parse(yStream, "2010/tab5.2b.txt");
InputStream zStream = getClass().getResourceAsStream(directory + "2010/tab5.2d.txt");
PoissonSeries zSeries = parser.parse(zStream, "2010/tab5.2d.txt");
final PoissonSeries.CompiledSeries compiled = PoissonSeries.compile(xSeries, ySeries, zSeries);
TimeScale ut1 = TimeScalesFactory.getUT1(FramesFactory.getEOPHistory(IERSConventions.IERS_2010, true));
final FundamentalNutationArguments arguments = IERSConventions.IERS_2010.getNutationArguments(ut1);
UnivariateDifferentiableVectorFunction finite = new FiniteDifferencesDifferentiator(4, 0.4).differentiate((double t) -> compiled.value(arguments.evaluateAll(AbsoluteDate.J2000_EPOCH.shiftedBy(t))));
DSFactory factory = new DSFactory(1, 1);
for (double t = 0; t < Constants.JULIAN_DAY; t += 120) {
// computation of derivatives from API
Decimal64[] dAPI = compiled.derivative(arguments.evaluateAll(FieldAbsoluteDate.getJ2000Epoch(Decimal64Field.getInstance()).shiftedBy(t)));
// finite differences computation of derivatives
DerivativeStructure[] d = finite.value(factory.variable(0, t));
Assert.assertEquals(d.length, dAPI.length);
for (int i = 0; i < d.length; ++i) {
Assert.assertEquals(d[i].getPartialDerivative(1), dAPI[i].getReal(), FastMath.abs(2.0e-7 * d[i].getPartialDerivative(1)));
}
}
}
use of org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction in project Orekit by CS-SI.
the class Differentiation method differentiate.
/**
* Differentiate a vector function using finite differences.
* @param function function to differentiate
* @param provider attitude provider to use for modified states
* @param dimension dimension of the vector value of the function
* @param orbitType type used to map the orbit to a one dimensional array
* @param positionAngle type of the position angle used for orbit mapping to array
* @param dP user specified position error, used for step size computation for finite differences
* @param nbPoints number of points used for finite differences
* @return matrix function evaluating to the Jacobian of the original function
*/
public static StateJacobian differentiate(final StateFunction function, final int dimension, final AttitudeProvider provider, final OrbitType orbitType, final PositionAngle positionAngle, final double dP, final int nbPoints) {
return new StateJacobian() {
@Override
public double[][] value(final SpacecraftState state) throws OrekitException {
try {
final double[] tolerances = NumericalPropagator.tolerances(dP, state.getOrbit(), orbitType)[0];
final double[][] jacobian = new double[dimension][6];
for (int j = 0; j < 6; ++j) {
// compute partial derivatives with respect to state component j
final UnivariateVectorFunction componentJ = new StateComponentFunction(j, function, provider, state, orbitType, positionAngle);
final FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(nbPoints, tolerances[j]);
final UnivariateDifferentiableVectorFunction differentiatedJ = differentiator.differentiate(componentJ);
final DerivativeStructure[] c = differentiatedJ.value(FACTORY.variable(0, 0.0));
// populate the j-th column of the Jacobian
for (int i = 0; i < dimension; ++i) {
jacobian[i][j] = c[i].getPartialDerivative(1);
}
}
return jacobian;
} catch (OrekitExceptionWrapper oew) {
throw oew.getException();
}
}
};
}
use of org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction in project Orekit by CS-SI.
the class FramesFactoryTest method doTestDerivatives.
private void doTestDerivatives(AbsoluteDate ref, double duration, double step, boolean forbidInterpolation, double cartesianTolerance, double cartesianDotTolerance, double cartesianDotDotTolerance, double rodriguesTolerance, double rodriguesDotTolerance, double rodriguesDotDotTolerance) throws OrekitException {
final DSFactory factory = new DSFactory(1, 2);
final FieldAbsoluteDate<DerivativeStructure> refDS = new FieldAbsoluteDate<>(factory.getDerivativeField(), ref);
FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(8, 60.0);
for (final Predefined predefined : Predefined.values()) {
final Frame frame = FramesFactory.getFrame(predefined);
final Frame parent = frame.getParent();
if (parent != null) {
UnivariateDifferentiableVectorFunction dCartesian = differentiator.differentiate(new UnivariateVectorFunction() {
@Override
public double[] value(double t) {
try {
return forbidInterpolation ? FramesFactory.getNonInterpolatingTransform(parent, frame, ref.shiftedBy(t)).getTranslation().toArray() : parent.getTransformTo(frame, ref.shiftedBy(t)).getTranslation().toArray();
} catch (OrekitException oe) {
throw new OrekitExceptionWrapper(oe);
}
}
});
UnivariateDifferentiableVectorFunction dOrientation = differentiator.differentiate(new UnivariateVectorFunction() {
double sign = +1.0;
Rotation previous = Rotation.IDENTITY;
@Override
public double[] value(double t) {
try {
AngularCoordinates ac = forbidInterpolation ? FramesFactory.getNonInterpolatingTransform(parent, frame, ref.shiftedBy(t)).getAngular() : parent.getTransformTo(frame, ref.shiftedBy(t)).getAngular();
final double dot = MathArrays.linearCombination(ac.getRotation().getQ0(), previous.getQ0(), ac.getRotation().getQ1(), previous.getQ1(), ac.getRotation().getQ2(), previous.getQ2(), ac.getRotation().getQ3(), previous.getQ3());
sign = FastMath.copySign(1.0, dot * sign);
previous = ac.getRotation();
return ac.getModifiedRodrigues(sign)[0];
} catch (OrekitException oe) {
throw new OrekitExceptionWrapper(oe);
}
}
});
double maxCartesianError = 0;
double maxCartesianDotError = 0;
double maxCartesianDotDotError = 0;
double maxRodriguesError = 0;
double maxRodriguesDotError = 0;
double maxRodriguesDotDotError = 0;
for (double dt = 0; dt < duration; dt += step) {
final DerivativeStructure dtDS = factory.variable(0, dt);
final FieldTransform<DerivativeStructure> tDS = forbidInterpolation ? FramesFactory.getNonInterpolatingTransform(parent, frame, refDS.shiftedBy(dtDS)) : parent.getTransformTo(frame, refDS.shiftedBy(dtDS));
final DerivativeStructure[] refCart = dCartesian.value(dtDS);
final DerivativeStructure[] fieldCart = tDS.getTranslation().toArray();
for (int i = 0; i < 3; ++i) {
maxCartesianError = FastMath.max(maxCartesianError, FastMath.abs(refCart[i].getValue() - fieldCart[i].getValue()));
maxCartesianDotError = FastMath.max(maxCartesianDotError, FastMath.abs(refCart[i].getPartialDerivative(1) - fieldCart[i].getPartialDerivative(1)));
maxCartesianDotDotError = FastMath.max(maxCartesianDotDotError, FastMath.abs(refCart[i].getPartialDerivative(2) - fieldCart[i].getPartialDerivative(2)));
}
final DerivativeStructure[] refOr = dOrientation.value(dtDS);
DerivativeStructure[] fieldOr = tDS.getAngular().getModifiedRodrigues(1.0)[0];
final double dot = refOr[0].linearCombination(refOr, fieldOr).getReal();
if (dot < 0 || Double.isNaN(dot)) {
fieldOr = tDS.getAngular().getModifiedRodrigues(-1.0)[0];
}
for (int i = 0; i < 3; ++i) {
maxRodriguesError = FastMath.max(maxRodriguesError, FastMath.abs(refOr[i].getValue() - fieldOr[i].getValue()));
maxRodriguesDotError = FastMath.max(maxRodriguesDotError, FastMath.abs(refOr[i].getPartialDerivative(1) - fieldOr[i].getPartialDerivative(1)));
maxRodriguesDotDotError = FastMath.max(maxRodriguesDotDotError, FastMath.abs(refOr[i].getPartialDerivative(2) - fieldOr[i].getPartialDerivative(2)));
}
}
Assert.assertEquals(frame.getName(), 0, maxCartesianError, cartesianTolerance);
Assert.assertEquals(frame.getName(), 0, maxCartesianDotError, cartesianDotTolerance);
Assert.assertEquals(frame.getName(), 0, maxCartesianDotDotError, cartesianDotDotTolerance);
Assert.assertEquals(frame.getName(), 0, maxRodriguesError, rodriguesTolerance);
Assert.assertEquals(frame.getName(), 0, maxRodriguesDotError, rodriguesDotTolerance);
Assert.assertEquals(frame.getName(), 0, maxRodriguesDotDotError, rodriguesDotDotTolerance);
}
}
}
Aggregations