Search in sources :

Example 6 with UnivariateDifferentiableVectorFunction

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)));
        }
    }
}
Also used : ByteArrayInputStream(java.io.ByteArrayInputStream) InputStream(java.io.InputStream) Decimal64(org.hipparchus.util.Decimal64) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) TimeScale(org.orekit.time.TimeScale) UnivariateDifferentiableVectorFunction(org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction) FiniteDifferencesDifferentiator(org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator) Test(org.junit.Test)

Example 7 with UnivariateDifferentiableVectorFunction

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();
            }
        }
    };
}
Also used : SpacecraftState(org.orekit.propagation.SpacecraftState) OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) UnivariateVectorFunction(org.hipparchus.analysis.UnivariateVectorFunction) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) FiniteDifferencesDifferentiator(org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator) UnivariateDifferentiableVectorFunction(org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction)

Example 8 with UnivariateDifferentiableVectorFunction

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);
        }
    }
}
Also used : OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) UnivariateVectorFunction(org.hipparchus.analysis.UnivariateVectorFunction) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) Rotation(org.hipparchus.geometry.euclidean.threed.Rotation) GeodeticPoint(org.orekit.bodies.GeodeticPoint) UnivariateDifferentiableVectorFunction(org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction) AngularCoordinates(org.orekit.utils.AngularCoordinates) OrekitException(org.orekit.errors.OrekitException) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) FiniteDifferencesDifferentiator(org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator)

Aggregations

DerivativeStructure (org.hipparchus.analysis.differentiation.DerivativeStructure)8 UnivariateDifferentiableVectorFunction (org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction)8 DSFactory (org.hipparchus.analysis.differentiation.DSFactory)7 FiniteDifferencesDifferentiator (org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator)6 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)5 UnivariateVectorFunction (org.hipparchus.analysis.UnivariateVectorFunction)4 Test (org.junit.Test)3 GeodeticPoint (org.orekit.bodies.GeodeticPoint)3 AbsoluteDate (org.orekit.time.AbsoluteDate)3 ByteArrayInputStream (java.io.ByteArrayInputStream)2 InputStream (java.io.InputStream)2 HashMap (java.util.HashMap)2 Rotation (org.hipparchus.geometry.euclidean.threed.Rotation)2 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)2 RandomGenerator (org.hipparchus.random.RandomGenerator)2 Well19937a (org.hipparchus.random.Well19937a)2 OneAxisEllipsoid (org.orekit.bodies.OneAxisEllipsoid)2 OrekitExceptionWrapper (org.orekit.errors.OrekitExceptionWrapper)2 Frame (org.orekit.frames.Frame)2 TopocentricFrame (org.orekit.frames.TopocentricFrame)2