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

Example 61 with DerivativeStructure

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

the class Differentiation method differentiate.

/**
 * Differentiate a scalar function using finite differences.
 * @param function function to differentiate
 * @param driver driver for the parameter
 * @param nbPoints number of points used for finite differences
 * @param step step for finite differences
 * @return scalar function evaluating to the derivative of the original function
 */
public static ParameterFunction differentiate(final ParameterFunction function, final ParameterDriver driver, final int nbPoints, final double step) {
    final UnivariateFunction uf = new UnivariateFunction() {

        /**
         * {@inheritDoc}
         */
        @Override
        public double value(final double normalizedValue) throws OrekitExceptionWrapper {
            try {
                final double saved = driver.getNormalizedValue();
                driver.setNormalizedValue(normalizedValue);
                final double functionValue = function.value(driver);
                driver.setNormalizedValue(saved);
                return functionValue;
            } catch (OrekitException oe) {
                throw new OrekitExceptionWrapper(oe);
            }
        }
    };
    final FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(nbPoints, step);
    final UnivariateDifferentiableFunction differentiated = differentiator.differentiate(uf);
    return new ParameterFunction() {

        /**
         * {@inheritDoc}
         */
        @Override
        public double value(final ParameterDriver parameterDriver) throws OrekitException {
            if (!parameterDriver.getName().equals(driver.getName())) {
                throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, parameterDriver.getName(), driver.getName());
            }
            try {
                final DerivativeStructure dsParam = FACTORY.variable(0, parameterDriver.getNormalizedValue());
                final DerivativeStructure dsValue = differentiated.value(dsParam);
                return dsValue.getPartialDerivative(1);
            } catch (OrekitExceptionWrapper oew) {
                throw oew.getException();
            }
        }
    };
}

Example 62 with DerivativeStructure

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

the class OneAxisEllipsoid method transform.

/**
 * Transform a Cartesian point to a surface-relative point.
 * @param point Cartesian point
 * @param frame frame in which Cartesian point is expressed
 * @param date date of the computation (used for frames conversions)
 * @return point at the same location but as a surface-relative point,
 * using time as the single derivation parameter
 * @exception OrekitException if point cannot be converted to body frame
 */
public FieldGeodeticPoint<DerivativeStructure> transform(final PVCoordinates point, final Frame frame, final AbsoluteDate date) throws OrekitException {
    // transform point to body frame
    final Transform toBody = frame.getTransformTo(bodyFrame, date);
    final PVCoordinates pointInBodyFrame = toBody.transformPVCoordinates(point);
    final FieldVector3D<DerivativeStructure> p = pointInBodyFrame.toDerivativeStructureVector(2);
    final DerivativeStructure pr2 = p.getX().multiply(p.getX()).add(p.getY().multiply(p.getY()));
    final DerivativeStructure pr = pr2.sqrt();
    final DerivativeStructure pz = p.getZ();
    // project point on the ellipsoid surface
    final TimeStampedPVCoordinates groundPoint = projectToGround(new TimeStampedPVCoordinates(date, pointInBodyFrame), bodyFrame);
    final FieldVector3D<DerivativeStructure> gp = groundPoint.toDerivativeStructureVector(2);
    final DerivativeStructure gpr2 = gp.getX().multiply(gp.getX()).add(gp.getY().multiply(gp.getY()));
    final DerivativeStructure gpr = gpr2.sqrt();
    final DerivativeStructure gpz = gp.getZ();
    // relative position of test point with respect to its ellipse sub-point
    final DerivativeStructure dr = pr.subtract(gpr);
    final DerivativeStructure dz = pz.subtract(gpz);
    final double insideIfNegative = g2 * (pr2.getReal() - ae2) + pz.getReal() * pz.getReal();
    return new FieldGeodeticPoint<>(DerivativeStructure.atan2(gpz, gpr.multiply(g2)), DerivativeStructure.atan2(p.getY(), p.getX()), DerivativeStructure.hypot(dr, dz).copySign(insideIfNegative));
}

Example 63 with DerivativeStructure

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

the class BoxAndSolarArraySpacecraft method dragAcceleration.

/**
 * {@inheritDoc}
 */
@Override
public FieldVector3D<DerivativeStructure> dragAcceleration(final AbsoluteDate date, final Frame frame, final Vector3D position, final Rotation rotation, final double mass, final double density, final Vector3D relativeVelocity, final double[] parameters, final String paramName) throws OrekitException {
    final DerivativeStructure dragCoeffDS;
    final DerivativeStructure liftRatioDS;
    final DerivativeStructure oMrDS;
    final Field<DerivativeStructure> field = factory.getDerivativeField();
    if (dragParameterDriver.getName().equals(paramName)) {
        final double liftRatio = liftParameterDriver == null ? 0.0 : parameters[1];
        dragCoeffDS = factory.variable(0, parameters[0]);
        liftRatioDS = factory.constant(liftRatio);
        oMrDS = factory.constant(1 - liftRatio);
    } else if (liftParameterDriver != null && liftParameterDriver.getName().equals(paramName)) {
        dragCoeffDS = factory.constant(parameters[0]);
        liftRatioDS = factory.variable(0, parameters[1]);
        oMrDS = liftRatioDS.negate().add(1);
    } else {
        if (liftParameterDriver == null) {
            throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, paramName, dragParameterDriver.getName());
        } else {
            throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, paramName, dragParameterDriver.getName(), liftParameterDriver.getName());
        }
    }
    // relative velocity in spacecraft frame
    final double vNorm2 = relativeVelocity.getNormSq();
    final double vNorm = FastMath.sqrt(vNorm2);
    final FieldVector3D<DerivativeStructure> vDir = new FieldVector3D<>(field, rotation.applyTo(relativeVelocity.scalarMultiply(1.0 / vNorm)));
    final DerivativeStructure coeff = dragCoeffDS.multiply(0.5 * density * vNorm2 / mass);
    // solar array facet contribution
    final FieldVector3D<DerivativeStructure> frontNormal = new FieldVector3D<>(field, getNormal(date, frame, position, rotation));
    final DerivativeStructure s = coeff.multiply(solarArrayArea).multiply(FieldVector3D.dotProduct(frontNormal, vDir));
    FieldVector3D<DerivativeStructure> acceleration = new FieldVector3D<>(s.abs().multiply(oMrDS), vDir, s.multiply(liftRatioDS).multiply(2), frontNormal);
    // body facets contribution
    for (final Facet facet : facets) {
        final DerivativeStructure dot = FieldVector3D.dotProduct(facet.getNormal(), vDir);
        if (dot.getValue() < 0) {
            // the facet intercepts the incoming flux
            final DerivativeStructure f = coeff.multiply(facet.getArea()).multiply(dot);
            acceleration = new FieldVector3D<>(field.getOne(), acceleration, f.abs().multiply(oMrDS), vDir, f.multiply(liftRatioDS).multiply(2), new FieldVector3D<>(field, facet.getNormal()));
        }
    }
    // convert back to inertial frame
    return new FieldRotation<>(field, rotation).applyInverseTo(acceleration);
}

Example 64 with DerivativeStructure

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

the class BoxAndSolarArraySpacecraft method getNormal.

/**
 * Get solar array normal in spacecraft frame.
 * @param date current date
 * @param frame inertial reference frame for state (both orbit and attitude)
 * @param position position of spacecraft in reference frame
 * @param rotation orientation (attitude) of the spacecraft with respect to reference frame
 * @return solar array normal in spacecraft frame
 * @exception OrekitException if sun direction cannot be computed in best lighting
 * configuration
 */
public synchronized FieldVector3D<DerivativeStructure> getNormal(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldRotation<DerivativeStructure> rotation) throws OrekitException {
    final DerivativeStructure zero = position.getX().getField().getZero();
    if (referenceDate != null) {
        // use a simple rotation at fixed rate
        final DerivativeStructure alpha = zero.add(rotationRate * date.durationFrom(referenceDate));
        return new FieldVector3D<>(alpha.cos(), saX, alpha.sin(), saY);
    }
    // compute orientation for best lighting
    final FieldVector3D<DerivativeStructure> sunInert = position.subtract(sun.getPVCoordinates(date, frame).getPosition()).negate().normalize();
    final FieldVector3D<DerivativeStructure> sunSpacecraft = rotation.applyTo(sunInert);
    final DerivativeStructure d = FieldVector3D.dotProduct(sunSpacecraft, saZ);
    final DerivativeStructure f = d.multiply(d).subtract(1).negate();
    if (f.getValue() < Precision.EPSILON) {
        // we set up an arbitrary normal
        return new FieldVector3D<>(position.getX().getField(), saZ.orthogonal());
    }
    final DerivativeStructure s = f.sqrt().reciprocal();
    return new FieldVector3D<>(s, sunSpacecraft, s.multiply(d).negate(), new FieldVector3D<>(zero.getField(), saZ));
}

Example 65 with DerivativeStructure

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

the class BoxAndSolarArraySpacecraft method radiationPressureAcceleration.

/**
 * {@inheritDoc}
 */
@Override
public FieldVector3D<DerivativeStructure> radiationPressureAcceleration(final AbsoluteDate date, final Frame frame, final Vector3D position, final Rotation rotation, final double mass, final Vector3D flux, final double[] parameters, final String paramName) throws OrekitException {
    if (flux.getNormSq() < Precision.SAFE_MIN) {
        // null illumination (we are probably in umbra)
        return FieldVector3D.getZero(factory.getDerivativeField());
    }
    final DerivativeStructure absorptionCoeffDS;
    final DerivativeStructure specularReflectionCoeffDS;
    if (ABSORPTION_COEFFICIENT.equals(paramName)) {
        absorptionCoeffDS = factory.variable(0, parameters[0]);
        specularReflectionCoeffDS = factory.constant(parameters[1]);
    } else if (REFLECTION_COEFFICIENT.equals(paramName)) {
        absorptionCoeffDS = factory.constant(parameters[0]);
        specularReflectionCoeffDS = factory.variable(0, parameters[1]);
    } else {
        throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, paramName, ABSORPTION_COEFFICIENT + ", " + REFLECTION_COEFFICIENT);
    }
    final DerivativeStructure diffuseReflectionCoeffDS = absorptionCoeffDS.add(specularReflectionCoeffDS).subtract(1).negate();
    // radiation flux in spacecraft frame
    final Vector3D fluxSat = rotation.applyTo(flux);
    // solar array contribution
    Vector3D normal = getNormal(date, frame, position, rotation);
    double dot = Vector3D.dotProduct(normal, fluxSat);
    if (dot > 0) {
        // the solar array is illuminated backward,
        // fix signs to compute contribution correctly
        dot = -dot;
        normal = normal.negate();
    }
    FieldVector3D<DerivativeStructure> force = facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot, specularReflectionCoeffDS, diffuseReflectionCoeffDS);
    // body facets contribution
    for (final Facet bodyFacet : facets) {
        normal = bodyFacet.getNormal();
        dot = Vector3D.dotProduct(normal, fluxSat);
        if (dot < 0) {
            // the facet intercepts the incoming flux
            force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot, specularReflectionCoeffDS, diffuseReflectionCoeffDS));
        }
    }
    // convert to inertial
    return FieldRotation.applyInverseTo(rotation, new FieldVector3D<>(1.0 / mass, force));
}