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

Example 6 with DSFactory

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

the class HolmesFeatherstoneAttractionModelTest method RealFieldTest.

/**
 *Testing if the propagation between the FieldPropagation and the propagation
 * is equivalent.
 * Also testing if propagating X+dX with the propagation is equivalent to
 * propagation X with the FieldPropagation and then applying the taylor
 * expansion of dX to the result.
 */
@Test
public void RealFieldTest() throws OrekitException {
    DSFactory factory = new DSFactory(6, 4);
    DerivativeStructure a_0 = factory.variable(0, 7201009.7124401);
    DerivativeStructure e_0 = factory.variable(1, 1e-3);
    DerivativeStructure i_0 = factory.variable(2, 98.7 * FastMath.PI / 180);
    DerivativeStructure R_0 = factory.variable(3, 15.0 * 22.5 * FastMath.PI / 180);
    DerivativeStructure O_0 = factory.variable(4, 93.0 * FastMath.PI / 180);
    DerivativeStructure n_0 = factory.variable(5, 0.1);
    Field<DerivativeStructure> field = a_0.getField();
    DerivativeStructure zero = field.getZero();
    FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
    Frame EME = FramesFactory.getEME2000();
    FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0, PositionAngle.MEAN, EME, J2000, Constants.EIGEN5C_EARTH_MU);
    FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
    SpacecraftState iSR = initialState.toSpacecraftState();
    OrbitType type = OrbitType.EQUINOCTIAL;
    double[][] tolerance = NumericalPropagator.tolerances(10.0, FKO.toOrbit(), type);
    AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(zero.add(60));
    AdaptiveStepsizeIntegrator RIntegrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
    RIntegrator.setInitialStepSize(60);
    FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
    FNP.setOrbitType(type);
    FNP.setInitialState(initialState);
    NumericalPropagator NP = new NumericalPropagator(RIntegrator);
    NP.setOrbitType(type);
    NP.setInitialState(iSR);
    double[][] c = new double[3][1];
    c[0][0] = 0.0;
    c[2][0] = normalizedC20;
    double[][] s = new double[3][1];
    NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(6378136.460, mu, TideSystem.UNKNOWN, c, s);
    HolmesFeatherstoneAttractionModel forceModel = new HolmesFeatherstoneAttractionModel(itrf, provider);
    FNP.addForceModel(forceModel);
    NP.addForceModel(forceModel);
    FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1005.);
    FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
    SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
    FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
    PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
    Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getX(), finPVC_R.getPosition().getX(), FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
    Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getY(), finPVC_R.getPosition().getY(), FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
    Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getZ(), finPVC_R.getPosition().getZ(), FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
    long number = 23091991;
    RandomGenerator RG = new Well19937a(number);
    GaussianRandomGenerator NGG = new GaussianRandomGenerator(RG);
    UncorrelatedRandomVectorGenerator URVG = new UncorrelatedRandomVectorGenerator(new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, new double[] { 1e1, 0.001, 0.001, 0.001, 0.001, 0.001 }, NGG);
    double a_R = a_0.getReal();
    double e_R = e_0.getReal();
    double i_R = i_0.getReal();
    double R_R = R_0.getReal();
    double O_R = O_0.getReal();
    double n_R = n_0.getReal();
    for (int ii = 0; ii < 1; ii++) {
        double[] rand_next = URVG.nextVector();
        double a_shift = a_R + rand_next[0];
        double e_shift = e_R + rand_next[1];
        double i_shift = i_R + rand_next[2];
        double R_shift = R_R + rand_next[3];
        double O_shift = O_R + rand_next[4];
        double n_shift = n_R + rand_next[5];
        KeplerianOrbit shiftedOrb = new KeplerianOrbit(a_shift, e_shift, i_shift, R_shift, O_shift, n_shift, PositionAngle.MEAN, EME, J2000.toAbsoluteDate(), Constants.EIGEN5C_EARTH_MU);
        SpacecraftState shift_iSR = new SpacecraftState(shiftedOrb);
        NumericalPropagator shift_NP = new NumericalPropagator(RIntegrator);
        shift_NP.setOrbitType(type);
        shift_NP.setInitialState(shift_iSR);
        shift_NP.addForceModel(forceModel);
        SpacecraftState finalState_shift = shift_NP.propagate(target.toAbsoluteDate());
        PVCoordinates finPVC_shift = finalState_shift.getPVCoordinates();
        // position check
        FieldVector3D<DerivativeStructure> pos_DS = finPVC_DS.getPosition();
        double x_DS = pos_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double y_DS = pos_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double z_DS = pos_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double x = finPVC_shift.getPosition().getX();
        double y = finPVC_shift.getPosition().getY();
        double z = finPVC_shift.getPosition().getZ();
        Assert.assertEquals(x_DS, x, FastMath.abs(x - pos_DS.getX().getReal()) * 1e-8);
        Assert.assertEquals(y_DS, y, FastMath.abs(y - pos_DS.getY().getReal()) * 1e-8);
        Assert.assertEquals(z_DS, z, FastMath.abs(z - pos_DS.getZ().getReal()) * 1e-8);
        // velocity check
        FieldVector3D<DerivativeStructure> vel_DS = finPVC_DS.getVelocity();
        double vx_DS = vel_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vy_DS = vel_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vz_DS = vel_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vx = finPVC_shift.getVelocity().getX();
        double vy = finPVC_shift.getVelocity().getY();
        double vz = finPVC_shift.getVelocity().getZ();
        Assert.assertEquals(vx_DS, vx, FastMath.abs(vx) * 1e-9);
        Assert.assertEquals(vy_DS, vy, FastMath.abs(vy) * 1e-9);
        Assert.assertEquals(vz_DS, vz, FastMath.abs(vz) * 1e-9);
        // acceleration check
        FieldVector3D<DerivativeStructure> acc_DS = finPVC_DS.getAcceleration();
        double ax_DS = acc_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double ay_DS = acc_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double az_DS = acc_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double ax = finPVC_shift.getAcceleration().getX();
        double ay = finPVC_shift.getAcceleration().getY();
        double az = finPVC_shift.getAcceleration().getZ();
        Assert.assertEquals(ax_DS, ax, FastMath.abs(ax) * 1e-9);
        Assert.assertEquals(ay_DS, ay, FastMath.abs(ay) * 1e-9);
        Assert.assertEquals(az_DS, az, FastMath.abs(az) * 1e-9);
    }
}

Example 7 with DSFactory

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

the class AbstractForceModelTest method toDS.

protected FieldSpacecraftState<DerivativeStructure> toDS(final SpacecraftState state, final AttitudeProvider attitudeProvider) throws OrekitException {
    final Vector3D p = state.getPVCoordinates().getPosition();
    final Vector3D v = state.getPVCoordinates().getVelocity();
    final Vector3D a = state.getPVCoordinates().getAcceleration();
    DSFactory factory = new DSFactory(6, 1);
    Field<DerivativeStructure> field = factory.getDerivativeField();
    final FieldAbsoluteDate<DerivativeStructure> fDate = new FieldAbsoluteDate<>(field, state.getDate());
    final TimeStampedFieldPVCoordinates<DerivativeStructure> fPVA = new TimeStampedFieldPVCoordinates<>(fDate, new FieldVector3D<>(factory.variable(0, p.getX()), factory.variable(1, p.getY()), factory.variable(2, p.getZ())), new FieldVector3D<>(factory.variable(3, v.getX()), factory.variable(4, v.getY()), factory.variable(5, v.getZ())), new FieldVector3D<>(factory.constant(a.getX()), factory.constant(a.getY()), factory.constant(a.getZ())));
    final FieldCartesianOrbit<DerivativeStructure> orbit = new FieldCartesianOrbit<>(fPVA, state.getFrame(), state.getMu());
    final FieldAttitude<DerivativeStructure> attitude = attitudeProvider.getAttitude(orbit, orbit.getDate(), orbit.getFrame());
    return new FieldSpacecraftState<>(orbit, attitude, field.getZero().add(state.getMass()));
}

Example 8 with DSFactory

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

the class PoissonSeriesParserTest method testDerivativesFromDoubleAPI.

@Test
public void testDerivativesFromDoubleAPI() 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
        double[] dAPI = compiled.derivative(arguments.evaluateAll(AbsoluteDate.J2000_EPOCH.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], FastMath.abs(2.0e-7 * d[i].getPartialDerivative(1)));
        }
    }
}

Example 9 with DSFactory

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

the class PoissonSeriesParserTest method testDerivativesAsField.

@Test
public void testDerivativesAsField() 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);
    PoissonSeries xSeries = parser.parse(getClass().getResourceAsStream(directory + "2010/tab5.2a.txt"), "2010/tab5.2a.txt");
    PoissonSeries ySeries = parser.parse(getClass().getResourceAsStream(directory + "2010/tab5.2b.txt"), "2010/tab5.2b.txt");
    PoissonSeries zSeries = parser.parse(getClass().getResourceAsStream(directory + "2010/tab5.2d.txt"), "2010/tab5.2d.txt");
    TimeScale ut1 = TimeScalesFactory.getUT1(FramesFactory.getEOPHistory(IERSConventions.IERS_2010, true));
    FundamentalNutationArguments arguments = IERSConventions.IERS_2010.getNutationArguments(ut1);
    Coordinate xCoordinate = new Coordinate(xSeries, arguments);
    Coordinate yCoordinate = new Coordinate(ySeries, arguments);
    Coordinate zCoordinate = new Coordinate(zSeries, arguments);
    UnivariateDifferentiableFunction dx = new FiniteDifferencesDifferentiator(4, 0.4).differentiate(xCoordinate);
    UnivariateDifferentiableFunction dy = new FiniteDifferencesDifferentiator(4, 0.4).differentiate(yCoordinate);
    UnivariateDifferentiableFunction dz = new FiniteDifferencesDifferentiator(4, 0.4).differentiate(zCoordinate);
    DSFactory factory = new DSFactory(1, 1);
    FieldAbsoluteDate<DerivativeStructure> ds2000 = FieldAbsoluteDate.getJ2000Epoch(factory.getDerivativeField());
    for (double t = 0; t < Constants.JULIAN_DAY; t += 120) {
        final FieldAbsoluteDate<DerivativeStructure> date = ds2000.shiftedBy(factory.variable(0, t));
        // direct computation of derivatives
        FieldBodiesElements<DerivativeStructure> elements = arguments.evaluateAll(date);
        Assert.assertEquals(0.0, elements.getDate().durationFrom(date).getValue(), 1.0e-15);
        DerivativeStructure xDirect = xSeries.value(elements);
        DerivativeStructure yDirect = ySeries.value(elements);
        DerivativeStructure zDirect = zSeries.value(elements);
        // finite differences computation of derivatives
        DerivativeStructure zero = factory.variable(0, 0.0);
        xCoordinate.setDate(date.toAbsoluteDate());
        DerivativeStructure xFinite = dx.value(zero);
        yCoordinate.setDate(date.toAbsoluteDate());
        DerivativeStructure yFinite = dy.value(zero);
        zCoordinate.setDate(date.toAbsoluteDate());
        DerivativeStructure zFinite = dz.value(zero);
        Assert.assertEquals(xFinite.getValue(), xDirect.getValue(), FastMath.abs(7.0e-15 * xFinite.getValue()));
        Assert.assertEquals(xFinite.getPartialDerivative(1), xDirect.getPartialDerivative(1), FastMath.abs(2.0e-07 * xFinite.getPartialDerivative(1)));
        Assert.assertEquals(yFinite.getValue(), yDirect.getValue(), FastMath.abs(7.0e-15 * yFinite.getValue()));
        Assert.assertEquals(yFinite.getPartialDerivative(1), yDirect.getPartialDerivative(1), FastMath.abs(2.0e-07 * yFinite.getPartialDerivative(1)));
        Assert.assertEquals(zFinite.getValue(), zDirect.getValue(), FastMath.abs(7.0e-15 * zFinite.getValue()));
        Assert.assertEquals(zFinite.getPartialDerivative(1), zDirect.getPartialDerivative(1), FastMath.abs(2.0e-07 * zFinite.getPartialDerivative(1)));
    }
}

Example 10 with DSFactory

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

the class DragForce method getDensityWrtStateUsingFiniteDifferences.

/**
 * Compute density and its derivatives.
 * Using finite differences for the derivatives.
 * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
 * <p>
 * From a theoretical point of view, this method computes the same values
 * as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
 * specific case of {@link DerivativeStructure} with respect to state, so
 * it is less general. However, it is *much* faster in this important case.
 * <p>
 * <p>
 * The derivatives should be computed with respect to position. The input
 * parameters already take into account the free parameters (6, 7 or 8 depending
 * on derivation with respect to drag coefficient and lift ratio being considered or not)
 * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
 * with respect to position. Free parameters at indices 3, 4 and 5 correspond
 * to derivatives with respect to velocity (these derivatives will remain zero
 * as the atmospheric density does not depend on velocity). Free parameter
 * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
 * and/or lift ratio (one of these or both).
 * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
 * </p>
 * @param date current date
 * @param frame inertial reference frame for state (both orbit and attitude)
 * @param position position of spacecraft in inertial frame
 * @param <T> type of the elements
 * @return the density and its derivatives
 * @exception OrekitException if derivatives cannot be computed
 * @since 9.0
 */
private <T extends RealFieldElement<T>> T getDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date, final Frame frame, final FieldVector3D<T> position) throws OrekitException {
    // Retrieve derivation properties for parameter T
    // It is implied here that T is a DerivativeStructure
    // With order 1 and 6, 7 or 8 free parameters
    // This is all checked before in method isStateDerivatives
    final DSFactory factory = ((DerivativeStructure) position.getX()).getFactory();
    // Build a DerivativeStructure using only derivatives with respect to position
    final DSFactory factory3 = new DSFactory(3, 1);
    final FieldVector3D<DerivativeStructure> position3 = new FieldVector3D<>(factory3.variable(0, position.getX().getReal()), factory3.variable(1, position.getY().getReal()), factory3.variable(2, position.getZ().getReal()));
    // Get atmosphere properties in atmosphere own frame
    final Frame atmFrame = atmosphere.getFrame();
    final Transform toBody = frame.getTransformTo(atmFrame, date);
    final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position3);
    final Vector3D posBody = posBodyDS.toVector3D();
    // Estimate density model by finite differences and composition
    // Using a delta of 1m
    final double delta = 1.0;
    final double x = posBody.getX();
    final double y = posBody.getY();
    final double z = posBody.getZ();
    final double rho0 = atmosphere.getDensity(date, posBody, atmFrame);
    final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y, z), atmFrame) - rho0) / delta;
    final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x, y + delta, z), atmFrame) - rho0) / delta;
    final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x, y, z + delta), atmFrame) - rho0) / delta;
    final double[] dXdQ = posBodyDS.getX().getAllDerivatives();
    final double[] dYdQ = posBodyDS.getY().getAllDerivatives();
    final double[] dZdQ = posBodyDS.getZ().getAllDerivatives();
    // Density with derivatives:
    // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
    // - Others are set to 0.
    final int p = factory.getCompiler().getFreeParameters();
    final double[] rhoAll = new double[p + 1];
    rhoAll[0] = rho0;
    for (int i = 1; i < 4; ++i) {
        rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
    }
    @SuppressWarnings("unchecked") final T rho = (T) (factory.build(rhoAll));
    return rho;
}