Search in sources :

Example 1 with UnivariateVectorFunction

use of org.hipparchus.analysis.UnivariateVectorFunction in project Orekit by CS-SI.

the class IERSConventionsTest method checkVectorFunctionConsistency.

private void checkVectorFunctionConsistency(final TimeVectorFunction function, final int dimension, final AbsoluteDate date, final double span, final double sampleStep, final double h, final double valueTolerance, final double derivativeTolerance) {
    UnivariateDifferentiableVectorFunction differentiated = new FiniteDifferencesDifferentiator(4, h).differentiate(new UnivariateVectorFunction() {

        @Override
        public double[] value(final double dt) {
            return function.value(date.shiftedBy(dt));
        }
    });
    DSFactory factory = new DSFactory(1, 1);
    FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(date, factory.constant(0.0));
    double maxValueError = 0;
    double maxDerivativeError = 0;
    for (double dt = 0; dt < span; dt += sampleStep) {
        DerivativeStructure dsdt = factory.variable(0, dt);
        DerivativeStructure[] yRef = differentiated.value(dsdt);
        DerivativeStructure[] y = function.value(dsDate.shiftedBy(dsdt));
        Assert.assertEquals(dimension, yRef.length);
        Assert.assertEquals(dimension, y.length);
        for (int i = 0; i < dimension; ++i) {
            maxValueError = FastMath.max(maxValueError, FastMath.abs(yRef[i].getValue() - y[i].getValue()));
            maxDerivativeError = FastMath.max(maxDerivativeError, FastMath.abs(yRef[i].getPartialDerivative(1) - y[i].getPartialDerivative(1)));
        }
    }
    Assert.assertEquals(0, maxValueError, valueTolerance);
    Assert.assertEquals(0, maxDerivativeError, derivativeTolerance);
}
Also used : UnivariateVectorFunction(org.hipparchus.analysis.UnivariateVectorFunction) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) FiniteDifferencesDifferentiator(org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator) UnivariateDifferentiableVectorFunction(org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction)

Example 2 with UnivariateVectorFunction

use of org.hipparchus.analysis.UnivariateVectorFunction in project Orekit by CS-SI.

the class PredefinedIAUPolesTest method testDerivatives.

@Test
public void testDerivatives() {
    final DSFactory factory = new DSFactory(1, 1);
    final AbsoluteDate ref = AbsoluteDate.J2000_EPOCH;
    final FieldAbsoluteDate<DerivativeStructure> refDS = new FieldAbsoluteDate<>(factory.getDerivativeField(), ref);
    FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(8, 60.0);
    for (final IAUPole iaupole : PredefinedIAUPoles.values()) {
        UnivariateDifferentiableVectorFunction dPole = differentiator.differentiate(new UnivariateVectorFunction() {

            @Override
            public double[] value(double t) {
                return iaupole.getPole(ref.shiftedBy(t)).toArray();
            }
        });
        UnivariateDifferentiableFunction dMeridian = differentiator.differentiate(new UnivariateFunction() {

            @Override
            public double value(double t) {
                return iaupole.getPrimeMeridianAngle(ref.shiftedBy(t));
            }
        });
        for (double dt = 0; dt < Constants.JULIAN_YEAR; dt += 3600) {
            final DerivativeStructure dtDS = factory.variable(0, dt);
            final DerivativeStructure[] refPole = dPole.value(dtDS);
            final DerivativeStructure[] fieldPole = iaupole.getPole(refDS.shiftedBy(dtDS)).toArray();
            for (int i = 0; i < 3; ++i) {
                Assert.assertEquals(refPole[i].getValue(), fieldPole[i].getValue(), 2.0e-15);
                Assert.assertEquals(refPole[i].getPartialDerivative(1), fieldPole[i].getPartialDerivative(1), 4.0e-17);
            }
            final DerivativeStructure refMeridian = dMeridian.value(dtDS);
            final DerivativeStructure fieldMeridian = iaupole.getPrimeMeridianAngle(refDS.shiftedBy(dtDS));
            Assert.assertEquals(refMeridian.getValue(), fieldMeridian.getValue(), 4.0e-12);
            Assert.assertEquals(refMeridian.getPartialDerivative(1), fieldMeridian.getPartialDerivative(1), 9.0e-14);
        }
    }
}
Also used : UnivariateVectorFunction(org.hipparchus.analysis.UnivariateVectorFunction) UnivariateFunction(org.hipparchus.analysis.UnivariateFunction) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbsoluteDate(org.orekit.time.AbsoluteDate) UnivariateDifferentiableVectorFunction(org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction) OldIAUPole(org.orekit.bodies.IAUPoleFactory.OldIAUPole) UnivariateDifferentiableFunction(org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) FiniteDifferencesDifferentiator(org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator) Test(org.junit.Test)

Example 3 with UnivariateVectorFunction

use of org.hipparchus.analysis.UnivariateVectorFunction 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 4 with UnivariateVectorFunction

use of org.hipparchus.analysis.UnivariateVectorFunction 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

UnivariateVectorFunction (org.hipparchus.analysis.UnivariateVectorFunction)4 DerivativeStructure (org.hipparchus.analysis.differentiation.DerivativeStructure)4 FiniteDifferencesDifferentiator (org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator)4 UnivariateDifferentiableVectorFunction (org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction)4 DSFactory (org.hipparchus.analysis.differentiation.DSFactory)3 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)3 OrekitExceptionWrapper (org.orekit.errors.OrekitExceptionWrapper)2 UnivariateFunction (org.hipparchus.analysis.UnivariateFunction)1 UnivariateDifferentiableFunction (org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction)1 Rotation (org.hipparchus.geometry.euclidean.threed.Rotation)1 Test (org.junit.Test)1 GeodeticPoint (org.orekit.bodies.GeodeticPoint)1 OldIAUPole (org.orekit.bodies.IAUPoleFactory.OldIAUPole)1 OrekitException (org.orekit.errors.OrekitException)1 SpacecraftState (org.orekit.propagation.SpacecraftState)1 AbsoluteDate (org.orekit.time.AbsoluteDate)1 AngularCoordinates (org.orekit.utils.AngularCoordinates)1