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

Example 11 with DerivativeStructure

use of org.hipparchus.analysis.differentiation.DerivativeStructure 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 12 with DerivativeStructure

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

the class HolmesFeatherstoneAttractionModel method isStateDerivative.

/**
 * Check if a field state corresponds to derivatives with respect to state.
 * @param state state to check
 * @param <T> type of the filed elements
 * @return true if state corresponds to derivatives with respect to state
 * @since 9.0
 */
private <T extends RealFieldElement<T>> boolean isStateDerivative(final FieldSpacecraftState<T> state) {
    try {
        final DerivativeStructure dsMass = (DerivativeStructure) state.getMass();
        final int o = dsMass.getOrder();
        final int p = dsMass.getFreeParameters();
        if (o != 1 || (p < 3)) {
            return false;
        }
        @SuppressWarnings("unchecked") final FieldPVCoordinates<DerivativeStructure> pv = (FieldPVCoordinates<DerivativeStructure>) state.getPVCoordinates();
        return isVariable(pv.getPosition().getX(), 0) && isVariable(pv.getPosition().getY(), 1) && isVariable(pv.getPosition().getZ(), 2);
    } catch (ClassCastException cce) {
        return false;
    }
}

Example 13 with DerivativeStructure

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

the class HolmesFeatherstoneAttractionModel method accelerationWrtState.

/**
 * Compute acceleration derivatives with respect to state parameters.
 * <p>
 * From a theoretical point of view, this method computes the same values
 * as {@link #acceleration(FieldSpacecraftState, RealFieldElement[])} 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 or 7 depending
 * on derivation with respect to mass 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 acceleration due to gravity does not depend on velocity). Free parameter
 * at index 6 (if present) corresponds to to derivatives with respect to mass
 * (this derivative will remain zero as acceleration due to gravity does not
 * depend on mass).
 * </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 mu central attraction coefficient to use
 * @return acceleration with all derivatives specified by the input parameters
 * own derivatives
 * @exception OrekitException if derivatives cannot be computed
 * @since 6.0
 */
private FieldVector3D<DerivativeStructure> accelerationWrtState(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final DerivativeStructure mu) throws OrekitException {
    // get the position in body frame
    final Transform fromBodyFrame = bodyFrame.getTransformTo(frame, date);
    final Transform toBodyFrame = fromBodyFrame.getInverse();
    final Vector3D positionBody = toBodyFrame.transformPosition(position.toVector3D());
    // compute gradient and Hessian
    final GradientHessian gh = gradientHessian(date, positionBody, mu.getReal());
    // gradient of the non-central part of the gravity field
    final double[] gInertial = fromBodyFrame.transformVector(new Vector3D(gh.getGradient())).toArray();
    // Hessian of the non-central part of the gravity field
    final RealMatrix hBody = new Array2DRowRealMatrix(gh.getHessian(), false);
    final RealMatrix rot = new Array2DRowRealMatrix(toBodyFrame.getRotation().getMatrix());
    final RealMatrix hInertial = rot.transpose().multiply(hBody).multiply(rot);
    // distribute all partial derivatives in a compact acceleration vector
    final double[] derivatives = new double[1 + position.getX().getFreeParameters()];
    final DerivativeStructure[] accDer = new DerivativeStructure[3];
    for (int i = 0; i < 3; ++i) {
        // first element is value of acceleration (i.e. gradient of field)
        derivatives[0] = gInertial[i];
        // next three elements are one row of the Jacobian of acceleration (i.e. Hessian of field)
        derivatives[1] = hInertial.getEntry(i, 0);
        derivatives[2] = hInertial.getEntry(i, 1);
        derivatives[3] = hInertial.getEntry(i, 2);
        // next element is derivative with respect to parameter mu
        if (derivatives.length > 4 && isVariable(mu, 3)) {
            derivatives[4] = gInertial[i] / mu.getReal();
        }
        accDer[i] = position.getX().getFactory().build(derivatives);
    }
    return new FieldVector3D<>(accDer);
}

Example 14 with DerivativeStructure

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

Example 15 with DerivativeStructure

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

the class IERSConventionsTest method checkScalarFunctionConsistency.

private void checkScalarFunctionConsistency(final TimeScalarFunction function, final AbsoluteDate date, final double span, final double sampleStep, final double h, final double valueTolerance, final double derivativeTolerance) {
    UnivariateDifferentiableFunction differentiated = new FiniteDifferencesDifferentiator(4, h).differentiate(new UnivariateFunction() {

        @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));
        maxValueError = FastMath.max(maxValueError, FastMath.abs(yRef.getValue() - y.getValue()));
        maxDerivativeError = FastMath.max(maxDerivativeError, FastMath.abs(yRef.getPartialDerivative(1) - y.getPartialDerivative(1)));
    }
    Assert.assertEquals(0, maxValueError, valueTolerance);
    Assert.assertEquals(0, maxDerivativeError, derivativeTolerance);
}