Examples with FiniteDifferencesDifferentiator org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator used on opensource projects

Example 16 with FiniteDifferencesDifferentiator

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

the class JB2008Test method testDensityGradient.

@Test
public void testDensityGradient() throws OrekitException {
    final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf);
    final JB2008 atm = new JB2008(new InputParams(), CelestialBodyFactory.getSun(), earth);
    final AbsoluteDate date = InputParams.TC[6];
    // Build the position
    final double alt = 400.;
    final double lat = 60.;
    final double lon = -70.;
    final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(lat), FastMath.toRadians(lon), alt * 1000.);
    final Vector3D pos = earth.transform(point);
    // Run
    DerivativeStructure zero = new DSFactory(1, 1).variable(0, 0.0);
    FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, 10.0);
    DerivativeStructure rhoX = differentiator.differentiate((double x) -> {
        try {
            return atm.getDensity(date, new Vector3D(1, pos, x, Vector3D.PLUS_I), itrf);
        } catch (OrekitException oe) {
            return Double.NaN;
        }
    }).value(zero);
    DerivativeStructure rhoY = differentiator.differentiate((double y) -> {
        try {
            return atm.getDensity(date, new Vector3D(1, pos, y, Vector3D.PLUS_J), itrf);
        } catch (OrekitException oe) {
            return Double.NaN;
        }
    }).value(zero);
    DerivativeStructure rhoZ = differentiator.differentiate((double z) -> {
        try {
            return atm.getDensity(date, new Vector3D(1, pos, z, Vector3D.PLUS_K), itrf);
        } catch (OrekitException oe) {
            return Double.NaN;
        }
    }).value(zero);
    DSFactory factory3 = new DSFactory(3, 1);
    Field<DerivativeStructure> field = factory3.getDerivativeField();
    final DerivativeStructure rhoDS = atm.getDensity(new FieldAbsoluteDate<>(field, date), new FieldVector3D<>(factory3.variable(0, pos.getX()), factory3.variable(1, pos.getY()), factory3.variable(2, pos.getZ())), itrf);
    Assert.assertEquals(rhoX.getValue(), rhoDS.getReal(), rhoX.getValue() * 2.0e-14);
    Assert.assertEquals(rhoY.getValue(), rhoDS.getReal(), rhoY.getValue() * 2.0e-14);
    Assert.assertEquals(rhoZ.getValue(), rhoDS.getReal(), rhoZ.getValue() * 2.0e-14);
    Assert.assertEquals(rhoX.getPartialDerivative(1), rhoDS.getPartialDerivative(1, 0, 0), FastMath.abs(6.0e-10 * rhoX.getPartialDerivative(1)));
    Assert.assertEquals(rhoY.getPartialDerivative(1), rhoDS.getPartialDerivative(0, 1, 0), FastMath.abs(6.0e-10 * rhoY.getPartialDerivative(1)));
    Assert.assertEquals(rhoZ.getPartialDerivative(1), rhoDS.getPartialDerivative(0, 0, 1), FastMath.abs(6.0e-10 * rhoY.getPartialDerivative(1)));
}

Example 17 with FiniteDifferencesDifferentiator

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

the class NRLMSISE00Test method testDensityGradient.

@Test
public void testDensityGradient() throws OrekitException {
    // Build the input params provider
    final InputParams ip = new InputParams();
    // Get Sun
    final PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
    // Get Earth body shape
    final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf);
    // Build the model
    final NRLMSISE00 atm = new NRLMSISE00(ip, sun, earth);
    // Build the date
    final AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 172), new TimeComponents(29000.), TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true));
    // Build the position
    final double alt = 400.;
    final double lat = 60.;
    final double lon = -70.;
    final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(lat), FastMath.toRadians(lon), alt * 1000.);
    final Vector3D pos = earth.transform(point);
    // Run
    DerivativeStructure zero = new DSFactory(1, 1).variable(0, 0.0);
    FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, 10.0);
    DerivativeStructure rhoX = differentiator.differentiate((double x) -> {
        try {
            return atm.getDensity(date, new Vector3D(1, pos, x, Vector3D.PLUS_I), itrf);
        } catch (OrekitException oe) {
            return Double.NaN;
        }
    }).value(zero);
    DerivativeStructure rhoY = differentiator.differentiate((double y) -> {
        try {
            return atm.getDensity(date, new Vector3D(1, pos, y, Vector3D.PLUS_J), itrf);
        } catch (OrekitException oe) {
            return Double.NaN;
        }
    }).value(zero);
    DerivativeStructure rhoZ = differentiator.differentiate((double z) -> {
        try {
            return atm.getDensity(date, new Vector3D(1, pos, z, Vector3D.PLUS_K), itrf);
        } catch (OrekitException oe) {
            return Double.NaN;
        }
    }).value(zero);
    DSFactory factory3 = new DSFactory(3, 1);
    Field<DerivativeStructure> field = factory3.getDerivativeField();
    final DerivativeStructure rhoDS = atm.getDensity(new FieldAbsoluteDate<>(field, date), new FieldVector3D<>(factory3.variable(0, pos.getX()), factory3.variable(1, pos.getY()), factory3.variable(2, pos.getZ())), itrf);
    Assert.assertEquals(rhoX.getValue(), rhoDS.getReal(), rhoX.getValue() * 2.0e-13);
    Assert.assertEquals(rhoY.getValue(), rhoDS.getReal(), rhoY.getValue() * 2.0e-13);
    Assert.assertEquals(rhoZ.getValue(), rhoDS.getReal(), rhoZ.getValue() * 2.0e-13);
    Assert.assertEquals(rhoX.getPartialDerivative(1), rhoDS.getPartialDerivative(1, 0, 0), FastMath.abs(2.0e-10 * rhoX.getPartialDerivative(1)));
    Assert.assertEquals(rhoY.getPartialDerivative(1), rhoDS.getPartialDerivative(0, 1, 0), FastMath.abs(2.0e-10 * rhoY.getPartialDerivative(1)));
    Assert.assertEquals(rhoZ.getPartialDerivative(1), rhoDS.getPartialDerivative(0, 0, 1), FastMath.abs(2.0e-10 * rhoY.getPartialDerivative(1)));
}

Example 18 with FiniteDifferencesDifferentiator

use of org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator 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);
        }
    }
}

Example 19 with FiniteDifferencesDifferentiator

use of org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator 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)));
        }
    }
}

Example 20 with FiniteDifferencesDifferentiator

use of org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator 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();
            }
        }
    };
}