Examples with Vector3D org.hipparchus.geometry.euclidean.threed.Vector3D used on opensource projects

Example 66 with Vector3D

use of org.hipparchus.geometry.euclidean.threed.Vector3D 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 67 with Vector3D

use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.

the class DTM2000 method getDensity.

/**
 * {@inheritDoc}
 */
@Override
public <T extends RealFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date, final FieldVector3D<T> position, final Frame frame) throws OrekitException {
    // check if data are available :
    final AbsoluteDate dateD = date.toAbsoluteDate();
    if ((dateD.compareTo(inputParams.getMaxDate()) > 0) || (dateD.compareTo(inputParams.getMinDate()) < 0)) {
        throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE, dateD, inputParams.getMinDate(), inputParams.getMaxDate());
    }
    // compute day number in current year
    final Calendar cal = new GregorianCalendar();
    cal.setTime(date.toDate(TimeScalesFactory.getUTC()));
    final int day = cal.get(Calendar.DAY_OF_YEAR);
    // position in ECEF so we only have to do the transform once
    final Frame ecef = earth.getBodyFrame();
    final FieldVector3D<T> pEcef = frame.getTransformTo(ecef, date).transformPosition(position);
    // compute geodetic position
    final FieldGeodeticPoint<T> inBody = earth.transform(pEcef, ecef, date);
    final T alti = inBody.getAltitude();
    final T lon = inBody.getLongitude();
    final T lat = inBody.getLatitude();
    // compute local solar time
    final Vector3D sunPos = sun.getPVCoordinates(dateD, ecef).getPosition();
    final T y = pEcef.getY().multiply(sunPos.getX()).subtract(pEcef.getX().multiply(sunPos.getY()));
    final T x = pEcef.getX().multiply(sunPos.getX()).add(pEcef.getY().multiply(sunPos.getY()));
    final T hl = y.atan2(x).add(FastMath.PI);
    // get current solar activity data and compute
    return getDensity(day, alti, lon, lat, hl, inputParams.getInstantFlux(dateD), inputParams.getMeanFlux(dateD), inputParams.getThreeHourlyKP(dateD), inputParams.get24HoursKp(dateD));
}

Example 68 with Vector3D

use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.

the class NRLMSISE00 method localSolarTime.

/**
 * Get local solar time.
 * @param date current date
 * @param position current position in frame
 * @param frame the frame in which is defined the position
 * @param <T> type of the filed elements
 * @return the local solar time (hour in [0, 24[)
 * @throws OrekitException if position cannot be computed in given frame
 */
private <T extends RealFieldElement<T>> T localSolarTime(final AbsoluteDate date, final FieldVector3D<T> position, final Frame frame) throws OrekitException {
    final Vector3D sunPos = sun.getPVCoordinates(date, frame).getPosition();
    final T y = position.getY().multiply(sunPos.getX()).subtract(position.getX().multiply(sunPos.getY()));
    final T x = position.getX().multiply(sunPos.getX()).add(position.getY().multiply(sunPos.getY()));
    final T hl = y.atan2(x).add(FastMath.PI);
    return hl.multiply(12. / FastMath.PI);
}

Example 69 with Vector3D

use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.

the class HarrisPriester method getDensity.

/**
 * Get the local density.
 * @param sunInEarth position of the Sun in Earth frame (m)
 * @param posInEarth target position in Earth frame (m)
 * @return the local density (kg/m³)
 * @exception OrekitException if altitude is below the model minimal altitude
 */
public double getDensity(final Vector3D sunInEarth, final Vector3D posInEarth) throws OrekitException {
    final double posAlt = getHeight(posInEarth);
    // Check for height boundaries
    if (posAlt < getMinAlt()) {
        throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, posAlt, getMinAlt());
    }
    if (posAlt > getMaxAlt()) {
        return 0.;
    }
    // Diurnal bulge apex direction
    final Vector3D sunDir = sunInEarth.normalize();
    final Vector3D bulDir = new Vector3D(sunDir.getX() * COSLAG - sunDir.getY() * SINLAG, sunDir.getX() * SINLAG + sunDir.getY() * COSLAG, sunDir.getZ());
    // Cosine of angle Psi between the diurnal bulge apex and the satellite
    final double cosPsi = bulDir.normalize().dotProduct(posInEarth.normalize());
    // (1 + cos(Psi))/2 = cos²(Psi/2)
    final double c2Psi2 = (1. + cosPsi) / 2.;
    final double cPsi2 = FastMath.sqrt(c2Psi2);
    final double cosPow = (cPsi2 > MIN_COS) ? c2Psi2 * FastMath.pow(cPsi2, n - 2) : 0.;
    // Search altitude index in density table
    int ia = 0;
    while (ia < tabAltRho.length - 2 && posAlt > tabAltRho[ia + 1][0]) {
        ia++;
    }
    // Fractional satellite height
    final double dH = (tabAltRho[ia][0] - posAlt) / (tabAltRho[ia][0] - tabAltRho[ia + 1][0]);
    // Min exponential density interpolation
    final double rhoMin = tabAltRho[ia][1] * FastMath.pow(tabAltRho[ia + 1][1] / tabAltRho[ia][1], dH);
    if (Precision.equals(cosPow, 0.)) {
        return rhoMin;
    } else {
        // Max exponential density interpolation
        final double rhoMax = tabAltRho[ia][2] * FastMath.pow(tabAltRho[ia + 1][2] / tabAltRho[ia][2], dH);
        return rhoMin + (rhoMax - rhoMin) * cosPow;
    }
}

Example 70 with Vector3D

use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.

the class HarrisPriester method getDensity.

/**
 * Get the local density.
 * @param sunInEarth position of the Sun in Earth frame (m)
 * @param posInEarth target position in Earth frame (m)
 * @return the local density (kg/m³)
 * @param <T> instance of RealFieldElement<T>
 * @exception OrekitException if altitude is below the model minimal altitude
 */
public <T extends RealFieldElement<T>> T getDensity(final Vector3D sunInEarth, final FieldVector3D<T> posInEarth) throws OrekitException {
    final T zero = posInEarth.getX().getField().getZero();
    final T posAlt = getHeight(posInEarth);
    // Check for height boundaries
    if (posAlt.getReal() < getMinAlt()) {
        throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, posAlt, getMinAlt());
    }
    if (posAlt.getReal() > getMaxAlt()) {
        return zero;
    }
    // Diurnal bulge apex direction
    final Vector3D sunDir = sunInEarth.normalize();
    final Vector3D bulDir = new Vector3D(sunDir.getX() * COSLAG - sunDir.getY() * SINLAG, sunDir.getX() * SINLAG + sunDir.getY() * COSLAG, sunDir.getZ());
    // Cosine of angle Psi between the diurnal bulge apex and the satellite
    final T cosPsi = posInEarth.normalize().dotProduct(bulDir.normalize());
    // (1 + cos(Psi))/2 = cos²(Psi/2)
    final T c2Psi2 = cosPsi.add(1.).divide(2);
    final T cPsi2 = c2Psi2.sqrt();
    final T cosPow = (cPsi2.getReal() > MIN_COS) ? c2Psi2.multiply(cPsi2.pow(n - 2)) : zero;
    // Search altitude index in density table
    int ia = 0;
    while (ia < tabAltRho.length - 2 && posAlt.getReal() > tabAltRho[ia + 1][0]) {
        ia++;
    }
    // Fractional satellite height
    final T dH = posAlt.negate().add(tabAltRho[ia][0]).divide(tabAltRho[ia][0] - tabAltRho[ia + 1][0]);
    // Min exponential density interpolation
    final T rhoMin = zero.add(tabAltRho[ia + 1][1] / tabAltRho[ia][1]).pow(dH).multiply(tabAltRho[ia][1]);
    if (Precision.equals(cosPow.getReal(), 0.)) {
        return zero.add(rhoMin);
    } else {
        // Max exponential density interpolation
        final T rhoMax = zero.add(tabAltRho[ia + 1][2] / tabAltRho[ia][2]).pow(dH).multiply(tabAltRho[ia][2]);
        return rhoMin.add(rhoMax.subtract(rhoMin).multiply(cosPow));
    }
}