Search in sources :

Example 56 with Vector3D

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

the class BatchLSEstimatorTest method testMultiSatWithParameters.

/**
 * A modified version of the previous test with a selection of propagation drivers to estimate
 *  One common (ยต)
 *  Some specifics for each satellite (Cr and Ca)
 *
 * @throws OrekitException
 */
@Test
public void testMultiSatWithParameters() throws OrekitException {
    // Test: Set the propagator drivers to estimate for each satellite
    final boolean muEstimated = true;
    final boolean crEstimated1 = true;
    final boolean caEstimated1 = true;
    final boolean crEstimated2 = true;
    final boolean caEstimated2 = false;
    // Builder sat 1
    final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder1 = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 1.0, Force.POTENTIAL, Force.SOLAR_RADIATION_PRESSURE);
    // Adding selection of parameters
    String satName = "sat 1";
    for (DelegatingDriver driver : propagatorBuilder1.getPropagationParametersDrivers().getDrivers()) {
        if (driver.getName().equals("central attraction coefficient")) {
            driver.setSelected(muEstimated);
        }
        if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
            driver.setName(driver.getName() + " " + satName);
            driver.setSelected(crEstimated1);
        }
        if (driver.getName().equals(RadiationSensitive.ABSORPTION_COEFFICIENT)) {
            driver.setName(driver.getName() + " " + satName);
            driver.setSelected(caEstimated1);
        }
    }
    // Builder for sat 2
    final Context context2 = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder2 = context2.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 1.0, Force.POTENTIAL, Force.SOLAR_RADIATION_PRESSURE);
    // Adding selection of parameters
    satName = "sat 2";
    for (ParameterDriver driver : propagatorBuilder2.getPropagationParametersDrivers().getDrivers()) {
        if (driver.getName().equals("central attraction coefficient")) {
            driver.setSelected(muEstimated);
        }
        if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
            driver.setName(driver.getName() + " " + satName);
            driver.setSelected(crEstimated2);
        }
        if (driver.getName().equals(RadiationSensitive.ABSORPTION_COEFFICIENT)) {
            driver.setName(driver.getName() + " " + satName);
            driver.setSelected(caEstimated2);
        }
    }
    // Create perfect inter-satellites range measurements
    final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
    final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(), original.getPosition().add(new Vector3D(1000, 2000, 3000)), original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))), context.initialOrbit.getFrame(), context.initialOrbit.getMu());
    final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit, propagatorBuilder2);
    closePropagator.setEphemerisMode();
    closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
    final BoundedPropagator ephemeris = closePropagator.getGeneratedEphemeris();
    Propagator propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder1);
    final List<ObservedMeasurement<?>> r12 = EstimationTestUtils.createMeasurements(propagator1, new InterSatellitesRangeMeasurementCreator(ephemeris), 1.0, 3.0, 300.0);
    // create perfect range measurements for first satellite
    propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder1);
    final List<ObservedMeasurement<?>> r1 = EstimationTestUtils.createMeasurements(propagator1, new RangeMeasurementCreator(context), 1.0, 3.0, 300.0);
    // create orbit estimator
    final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder1, propagatorBuilder2);
    for (final ObservedMeasurement<?> interSat : r12) {
        estimator.addMeasurement(interSat);
    }
    for (final ObservedMeasurement<?> range : r1) {
        estimator.addMeasurement(range);
    }
    estimator.setParametersConvergenceThreshold(1.0e-2);
    estimator.setMaxIterations(10);
    estimator.setMaxEvaluations(20);
    estimator.setObserver(new BatchLSObserver() {

        int lastIter = 0;

        int lastEval = 0;

        /**
         * {@inheritDoc}
         */
        @Override
        public void evaluationPerformed(int iterationsCount, int evaluationscount, Orbit[] orbits, ParameterDriversList estimatedOrbitalParameters, ParameterDriversList estimatedPropagatorParameters, ParameterDriversList estimatedMeasurementsParameters, EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) throws OrekitException {
            if (iterationsCount == lastIter) {
                Assert.assertEquals(lastEval + 1, evaluationscount);
            } else {
                Assert.assertEquals(lastIter + 1, iterationsCount);
            }
            lastIter = iterationsCount;
            lastEval = evaluationscount;
            AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
            for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
                AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
                Assert.assertTrue(current.compareTo(previous) >= 0);
                previous = current;
            }
        }
    });
    List<DelegatingDriver> parameters = estimator.getOrbitalParametersDrivers(true).getDrivers();
    ParameterDriver a0Driver = parameters.get(0);
    Assert.assertEquals("a[0]", a0Driver.getName());
    a0Driver.setValue(a0Driver.getValue() + 1.2);
    a0Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
    ParameterDriver a1Driver = parameters.get(6);
    Assert.assertEquals("a[1]", a1Driver.getName());
    a1Driver.setValue(a1Driver.getValue() - 5.4);
    a1Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
    final Orbit before = new KeplerianOrbit(parameters.get(6).getValue(), parameters.get(7).getValue(), parameters.get(8).getValue(), parameters.get(9).getValue(), parameters.get(10).getValue(), parameters.get(11).getValue(), PositionAngle.TRUE, closeOrbit.getFrame(), closeOrbit.getDate(), closeOrbit.getMu());
    Assert.assertEquals(4.7246, Vector3D.distance(closeOrbit.getPVCoordinates().getPosition(), before.getPVCoordinates().getPosition()), 1.0e-3);
    Assert.assertEquals(0.0010514, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(), before.getPVCoordinates().getVelocity()), 1.0e-6);
    EstimationTestUtils.checkFit(context, estimator, 4, 5, 0.0, 6.0e-06, 0.0, 1.7e-05, 0.0, 4.4e-07, 0.0, 1.7e-10);
    final Orbit determined = new KeplerianOrbit(parameters.get(6).getValue(), parameters.get(7).getValue(), parameters.get(8).getValue(), parameters.get(9).getValue(), parameters.get(10).getValue(), parameters.get(11).getValue(), PositionAngle.TRUE, closeOrbit.getFrame(), closeOrbit.getDate(), closeOrbit.getMu());
    Assert.assertEquals(0.0, Vector3D.distance(closeOrbit.getPVCoordinates().getPosition(), determined.getPVCoordinates().getPosition()), 5.8e-6);
    Assert.assertEquals(0.0, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(), determined.getPVCoordinates().getVelocity()), 3.5e-9);
    // got a default one
    for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
        if (driver.getName().startsWith("a[")) {
            // user-specified reference date
            Assert.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
        } else {
            // default reference date
            Assert.assertEquals(0, driver.getReferenceDate().durationFrom(propagatorBuilder1.getInitialOrbitDate()), 1.0e-15);
        }
    }
}
Also used : CartesianOrbit(org.orekit.orbits.CartesianOrbit) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) AbsoluteDate(org.orekit.time.AbsoluteDate) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) ParameterDriversList(org.orekit.utils.ParameterDriversList) BoundedPropagator(org.orekit.propagation.BoundedPropagator) Propagator(org.orekit.propagation.Propagator) OrekitException(org.orekit.errors.OrekitException) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) DelegatingDriver(org.orekit.utils.ParameterDriversList.DelegatingDriver) BoundedPropagator(org.orekit.propagation.BoundedPropagator) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) EstimationsProvider(org.orekit.estimation.measurements.EstimationsProvider) Context(org.orekit.estimation.Context) Evaluation(org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation) Orbit(org.orekit.orbits.Orbit) CartesianOrbit(org.orekit.orbits.CartesianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) ParameterDriver(org.orekit.utils.ParameterDriver) LevenbergMarquardtOptimizer(org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) InterSatellitesRangeMeasurementCreator(org.orekit.estimation.measurements.InterSatellitesRangeMeasurementCreator) RangeMeasurementCreator(org.orekit.estimation.measurements.RangeMeasurementCreator) InterSatellitesRangeMeasurementCreator(org.orekit.estimation.measurements.InterSatellitesRangeMeasurementCreator) Test(org.junit.Test)

Example 57 with Vector3D

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

the class PVTest method testPVWithCovarianceMatrix.

/**
 * Test the PV constructor with one 6x6 covariance matrix as input.
 * @throws OrekitException
 */
@Test
public void testPVWithCovarianceMatrix() throws OrekitException {
    // Context
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    // Dummy P, V, T
    final Vector3D position = context.initialOrbit.getPVCoordinates().getPosition();
    final Vector3D velocity = context.initialOrbit.getPVCoordinates().getVelocity();
    final AbsoluteDate date = context.initialOrbit.getDate();
    // Initialize standard deviations, weight and corr coeff
    final double[] sigmaPV = { 10., 20., 30., 0.1, 0.2, 0.3 };
    final double baseWeight = 0.5;
    final double[][] corrCoefRef = new double[6][6];
    for (int i = 0; i < 6; i++) {
        for (int j = 0; j < 6; j++) {
            if (i == j) {
                corrCoefRef[i][i] = 1.;
            } else {
                corrCoefRef[i][j] = i + j + 1;
            }
        }
    }
    // Reference covariance matrix
    final double[][] Pref = new double[6][6];
    for (int i = 0; i < 6; i++) {
        for (int j = 0; j < 6; j++) {
            Pref[i][j] = corrCoefRef[i][j] * sigmaPV[i] * sigmaPV[j];
        }
    }
    // Reference propagator numbers
    final int[] propNumRef = { 0, 2 };
    // Create PV measurements
    final PV[] pvs = new PV[2];
    pvs[0] = new PV(date, position, velocity, Pref, baseWeight);
    pvs[1] = new PV(date, position, velocity, Pref, baseWeight, propNumRef[1]);
    // Tolerance
    // tolerance
    final double eps = 1.8e-15;
    // Check data
    for (int k = 0; k < pvs.length; k++) {
        final PV pv = pvs[k];
        // Propagator numbers
        assertEquals(propNumRef[k], pv.getPropagatorsIndices().get(0), eps);
        // Weights
        for (int i = 0; i < 6; i++) {
            assertEquals(baseWeight, pv.getBaseWeight()[i], eps);
        }
        // Sigmas
        for (int i = 0; i < 6; i++) {
            assertEquals(sigmaPV[i], pv.getTheoreticalStandardDeviation()[i], eps);
        }
        // Covariances
        final double[][] P = pv.getCovarianceMatrix();
        // Substract with ref and get the norm
        final double normP = MatrixUtils.createRealMatrix(P).subtract(MatrixUtils.createRealMatrix(Pref)).getNorm();
        assertEquals(0., normP, eps);
        // Correlation coef
        final double[][] corrCoef = pv.getCorrelationCoefficientsMatrix();
        // Substract with ref and get the norm
        final double normCorrCoef = MatrixUtils.createRealMatrix(corrCoef).subtract(MatrixUtils.createRealMatrix(corrCoefRef)).getNorm();
        assertEquals(0., normCorrCoef, eps);
    }
}
Also used : Context(org.orekit.estimation.Context) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) AbsoluteDate(org.orekit.time.AbsoluteDate) Test(org.junit.Test)

Example 58 with Vector3D

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

the class PVTest method testExceptions.

/**
 * Test exceptions raised if the covariance matrix does not have the proper size.
 */
@Test
public void testExceptions() throws OrekitException {
    // Context
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    // Dummy P, V, T
    final Vector3D position = context.initialOrbit.getPVCoordinates().getPosition();
    final Vector3D velocity = context.initialOrbit.getPVCoordinates().getVelocity();
    final AbsoluteDate date = context.initialOrbit.getDate();
    final double weight = 1.;
    // Build with two 3-sized vectors
    try {
        new PV(date, position, velocity, new double[] { 0., 0., 0. }, new double[] { 1. }, weight);
        Assert.fail("An OrekitException should have been thrown");
    } catch (OrekitException e) {
    // An exception should indeed be raised here
    }
    // Build with one 6-sized vector
    try {
        new PV(date, position, velocity, new double[] { 0., 0., 0. }, weight);
        Assert.fail("An OrekitException should have been thrown");
    } catch (OrekitException e) {
    // An exception should indeed be raised here
    }
    // Build with two 3x3 matrices
    try {
        new PV(date, position, velocity, new double[][] { { 0., 0. }, { 0., 0. } }, new double[][] { { 0., 0. }, { 0., 0. } }, weight);
        Assert.fail("An OrekitException should have been thrown");
    } catch (OrekitException e) {
    // An exception should indeed be raised here
    }
    // Build with one 6x6 matrix
    try {
        new PV(date, position, velocity, new double[][] { { 0., 0. }, { 0., 0. } }, weight);
        Assert.fail("An OrekitException should have been thrown");
    } catch (OrekitException e) {
    // An exception should indeed be raised here
    }
}
Also used : Context(org.orekit.estimation.Context) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) OrekitException(org.orekit.errors.OrekitException) AbsoluteDate(org.orekit.time.AbsoluteDate) Test(org.junit.Test)

Example 59 with Vector3D

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

the class RangeAnalytic method theoreticalEvaluationAnalytic.

/**
 * Analytical version of the theoreticalEvaluation function in Range class
 * The derivative structures are not used, an analytical computation is used instead.
 * @param iteration current LS estimator iteration
 * @param evaluation current LS estimator evaluation
 * @param state spacecraft state. At measurement date on first iteration then close to emission date on further iterations
 * @param interpolator Orekit step interpolator
 * @return
 * @throws OrekitException
 */
protected EstimatedMeasurement<Range> theoreticalEvaluationAnalytic(final int iteration, final int evaluation, final SpacecraftState state) throws OrekitException {
    // Station attribute from parent Range class
    final GroundStation groundStation = this.getStation();
    // Station position at signal arrival
    final AbsoluteDate downlinkDate = getDate();
    final Transform topoToInertDownlink = groundStation.getOffsetToInertial(state.getFrame(), downlinkDate);
    final TimeStampedPVCoordinates stationDownlink = topoToInertDownlink.transformPVCoordinates(new TimeStampedPVCoordinates(downlinkDate, PVCoordinates.ZERO));
    // Take propagation time into account
    // (if state has already been set up to pre-compensate propagation delay,
    // we will have offset == downlinkDelay and transitState will be
    // the same as state)
    // Downlink time of flight
    final double tauD = signalTimeOfFlight(state.getPVCoordinates(), stationDownlink.getPosition(), downlinkDate);
    final double delta = downlinkDate.durationFrom(state.getDate());
    final double dt = delta - tauD;
    // Transit state position
    final SpacecraftState transitState = state.shiftedBy(dt);
    final AbsoluteDate transitDate = transitState.getDate();
    final Vector3D transitP = transitState.getPVCoordinates().getPosition();
    // Station position at transit state date
    final Transform topoToInertAtTransitDate = groundStation.getOffsetToInertial(state.getFrame(), transitDate);
    final TimeStampedPVCoordinates stationAtTransitDate = topoToInertAtTransitDate.transformPVCoordinates(new TimeStampedPVCoordinates(transitDate, PVCoordinates.ZERO));
    // Uplink time of flight
    final double tauU = signalTimeOfFlight(stationAtTransitDate, transitP, transitDate);
    final double tau = tauD + tauU;
    // Real date and position of station at signal departure
    final AbsoluteDate uplinkDate = downlinkDate.shiftedBy(-tau);
    final TimeStampedPVCoordinates stationUplink = topoToInertDownlink.shiftedBy(-tau).transformPVCoordinates(new TimeStampedPVCoordinates(uplinkDate, PVCoordinates.ZERO));
    // Prepare the evaluation
    final EstimatedMeasurement<Range> estimated = new EstimatedMeasurement<Range>(this, iteration, evaluation, new SpacecraftState[] { transitState }, new TimeStampedPVCoordinates[] { stationUplink, transitState.getPVCoordinates(), stationDownlink });
    // Set range value
    final double cOver2 = 0.5 * Constants.SPEED_OF_LIGHT;
    estimated.setEstimatedValue(tau * cOver2);
    // Partial derivatives with respect to state
    // The formulas below take into account the fact the measurement is at fixed reception date.
    // When spacecraft position is changed, the downlink delay is changed, and in order
    // to still have the measurement arrive at exactly the same date on ground, we must
    // take the spacecraft-station relative velocity into account.
    final Vector3D v = state.getPVCoordinates().getVelocity();
    final Vector3D qv = stationDownlink.getVelocity();
    final Vector3D downInert = stationDownlink.getPosition().subtract(transitP);
    final double dDown = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT * tauD - Vector3D.dotProduct(downInert, v);
    final Vector3D upInert = transitP.subtract(stationUplink.getPosition());
    // test
    // final double   dUp       = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT * tauU -
    // Vector3D.dotProduct(upInert, qv);
    // test
    final double dUp = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT * tauU - Vector3D.dotProduct(upInert, stationUplink.getVelocity());
    // derivatives of the downlink time of flight
    final double dTauDdPx = -downInert.getX() / dDown;
    final double dTauDdPy = -downInert.getY() / dDown;
    final double dTauDdPz = -downInert.getZ() / dDown;
    // Derivatives of the uplink time of flight
    final double dTauUdPx = 1. / dUp * upInert.dotProduct(Vector3D.PLUS_I.add((qv.subtract(v)).scalarMultiply(dTauDdPx)));
    final double dTauUdPy = 1. / dUp * upInert.dotProduct(Vector3D.PLUS_J.add((qv.subtract(v)).scalarMultiply(dTauDdPy)));
    final double dTauUdPz = 1. / dUp * upInert.dotProduct(Vector3D.PLUS_K.add((qv.subtract(v)).scalarMultiply(dTauDdPz)));
    // derivatives of the range measurement
    final double dRdPx = (dTauDdPx + dTauUdPx) * cOver2;
    final double dRdPy = (dTauDdPy + dTauUdPy) * cOver2;
    final double dRdPz = (dTauDdPz + dTauUdPz) * cOver2;
    estimated.setStateDerivatives(0, new double[] { dRdPx, dRdPy, dRdPz, dRdPx * dt, dRdPy * dt, dRdPz * dt });
    if (groundStation.getEastOffsetDriver().isSelected() || groundStation.getNorthOffsetDriver().isSelected() || groundStation.getZenithOffsetDriver().isSelected()) {
        // Downlink tme of flight derivatives / station position in topocentric frame
        final AngularCoordinates ac = topoToInertDownlink.getAngular().revert();
        // final Rotation rotTopoToInert = ac.getRotation();
        final Vector3D omega = ac.getRotationRate();
        // Inertial frame
        final double dTauDdQIx = downInert.getX() / dDown;
        final double dTauDdQIy = downInert.getY() / dDown;
        final double dTauDdQIz = downInert.getZ() / dDown;
        // Uplink tme of flight derivatives / station position in topocentric frame
        // Inertial frame
        final double dTauUdQIx = 1 / dUp * upInert.dotProduct(Vector3D.MINUS_I.add((qv.subtract(v)).scalarMultiply(dTauDdQIx)).subtract(Vector3D.PLUS_I.crossProduct(omega).scalarMultiply(tau)));
        final double dTauUdQIy = 1 / dUp * upInert.dotProduct(Vector3D.MINUS_J.add((qv.subtract(v)).scalarMultiply(dTauDdQIy)).subtract(Vector3D.PLUS_J.crossProduct(omega).scalarMultiply(tau)));
        final double dTauUdQIz = 1 / dUp * upInert.dotProduct(Vector3D.MINUS_K.add((qv.subtract(v)).scalarMultiply(dTauDdQIz)).subtract(Vector3D.PLUS_K.crossProduct(omega).scalarMultiply(tau)));
        // Range partial derivatives
        // with respect to station position in inertial frame
        final Vector3D dRdQI = new Vector3D((dTauDdQIx + dTauUdQIx) * cOver2, (dTauDdQIy + dTauUdQIy) * cOver2, (dTauDdQIz + dTauUdQIz) * cOver2);
        // convert to topocentric frame, as the station position
        // offset parameter is expressed in this frame
        final Vector3D dRdQT = ac.getRotation().applyTo(dRdQI);
        if (groundStation.getEastOffsetDriver().isSelected()) {
            estimated.setParameterDerivatives(groundStation.getEastOffsetDriver(), dRdQT.getX());
        }
        if (groundStation.getNorthOffsetDriver().isSelected()) {
            estimated.setParameterDerivatives(groundStation.getNorthOffsetDriver(), dRdQT.getY());
        }
        if (groundStation.getZenithOffsetDriver().isSelected()) {
            estimated.setParameterDerivatives(groundStation.getZenithOffsetDriver(), dRdQT.getZ());
        }
    }
    return estimated;
}
Also used : SpacecraftState(org.orekit.propagation.SpacecraftState) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) AngularCoordinates(org.orekit.utils.AngularCoordinates) Transform(org.orekit.frames.Transform) FieldTransform(org.orekit.frames.FieldTransform) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbsoluteDate(org.orekit.time.AbsoluteDate)

Example 60 with Vector3D

use of org.hipparchus.geometry.euclidean.threed.Vector3D 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);
}
Also used : RealMatrix(org.hipparchus.linear.RealMatrix) Array2DRowRealMatrix(org.hipparchus.linear.Array2DRowRealMatrix) Array2DRowRealMatrix(org.hipparchus.linear.Array2DRowRealMatrix) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) Transform(org.orekit.frames.Transform) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D)

Aggregations

Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)750 Test (org.junit.Test)466 AbsoluteDate (org.orekit.time.AbsoluteDate)323 PVCoordinates (org.orekit.utils.PVCoordinates)280 FieldVector3D (org.hipparchus.geometry.euclidean.threed.FieldVector3D)259 TimeStampedPVCoordinates (org.orekit.utils.TimeStampedPVCoordinates)187 SpacecraftState (org.orekit.propagation.SpacecraftState)152 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)124 Rotation (org.hipparchus.geometry.euclidean.threed.Rotation)119 Frame (org.orekit.frames.Frame)115 KeplerianOrbit (org.orekit.orbits.KeplerianOrbit)105 Orbit (org.orekit.orbits.Orbit)100 GeodeticPoint (org.orekit.bodies.GeodeticPoint)84 OrekitException (org.orekit.errors.OrekitException)83 CartesianOrbit (org.orekit.orbits.CartesianOrbit)75 EquinoctialOrbit (org.orekit.orbits.EquinoctialOrbit)68 DateComponents (org.orekit.time.DateComponents)67 Transform (org.orekit.frames.Transform)61 OneAxisEllipsoid (org.orekit.bodies.OneAxisEllipsoid)59 CircularOrbit (org.orekit.orbits.CircularOrbit)59