Examples with ForceModel org.orekit.forces.ForceModel used on opensource projects

Example 11 with ForceModel

use of org.orekit.forces.ForceModel in project Orekit by CS-SI.

the class FieldPropagation method main.

/**
 * Program entry point.
 * @param args program arguments (unused here)
 * @throws IOException
 * @throws OrekitException
 */
public static void main(String[] args) throws IOException, OrekitException {
    // the goal of this example is to make a Montecarlo simulation giving an error on the semiaxis,
    // the inclination and the RAAN. The interest of doing it with Orekit based on the
    // DerivativeStructure is that instead of doing a large number of propagation around the initial
    // point we will do a single propagation of the initial state, and thanks to the Taylor expansion
    // we will see the evolution of the std deviation of the position, which is divided in the
    // CrossTrack, the LongTrack and the Radial error.
    // configure Orekit
    File home = new File(System.getProperty("user.home"));
    File orekitData = new File(home, "orekit-data");
    if (!orekitData.exists()) {
        System.err.format(Locale.US, "Failed to find %s folder%n", orekitData.getAbsolutePath());
        System.err.format(Locale.US, "You need to download %s from the %s page and unzip it in %s for this tutorial to work%n", "orekit-data.zip", "https://www.orekit.org/forge/projects/orekit/files", home.getAbsolutePath());
        System.exit(1);
    }
    DataProvidersManager manager = DataProvidersManager.getInstance();
    manager.addProvider(new DirectoryCrawler(orekitData));
    // output file in user's home directory
    File workingDir = new File(System.getProperty("user.home"));
    File errorFile = new File(workingDir, "error.txt");
    System.out.println("Output file is in : " + errorFile.getAbsolutePath());
    PrintWriter PW = new PrintWriter(errorFile, "UTF-8");
    PW.printf("time \t\tCrossTrackErr \tLongTrackErr  \tRadialErr \tTotalErr%n");
    // setting the parameters of the simulation
    // Order of derivation of the DerivativeStructures
    int params = 3;
    int order = 3;
    DSFactory factory = new DSFactory(params, order);
    // number of samples of the montecarlo simulation
    int montecarlo_size = 100;
    // nominal values of the Orbital parameters
    double a_nominal = 7.278E6;
    double e_nominal = 1e-3;
    double i_nominal = FastMath.toRadians(98.3);
    double pa_nominal = FastMath.PI / 2;
    double raan_nominal = 0.0;
    double ni_nominal = 0.0;
    // mean of the gaussian curve for each of the errors around the nominal values
    // {a, i, RAAN}
    double[] mean = { 0, 0, 0 };
    // standard deviation of the gaussian curve for each of the errors around the nominal values
    // {dA, dI, dRaan}
    double[] dAdIdRaan = { 5, FastMath.toRadians(1e-3), FastMath.toRadians(1e-3) };
    // time of integration
    double final_Dt = 1 * 60 * 60;
    // number of steps per orbit
    double num_step_orbit = 10;
    DerivativeStructure a_0 = factory.variable(0, a_nominal);
    DerivativeStructure e_0 = factory.constant(e_nominal);
    DerivativeStructure i_0 = factory.variable(1, i_nominal);
    DerivativeStructure pa_0 = factory.constant(pa_nominal);
    DerivativeStructure raan_0 = factory.variable(2, raan_nominal);
    DerivativeStructure ni_0 = factory.constant(ni_nominal);
    // sometimes we will need the field of the DerivativeStructure to build new instances
    Field<DerivativeStructure> field = a_0.getField();
    // sometimes we will need the zero of the DerivativeStructure to build new instances
    DerivativeStructure zero = field.getZero();
    // initializing the FieldAbsoluteDate with only the field it will generate the day J2000
    FieldAbsoluteDate<DerivativeStructure> date_0 = new FieldAbsoluteDate<>(field);
    // initialize a basic frame
    Frame frame = FramesFactory.getEME2000();
    // initialize the orbit
    double mu = 3.9860047e14;
    FieldKeplerianOrbit<DerivativeStructure> KO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, pa_0, raan_0, ni_0, PositionAngle.ECCENTRIC, frame, date_0, mu);
    // step of integration (how many times per orbit we take the mesures)
    double int_step = KO.getKeplerianPeriod().getReal() / num_step_orbit;
    // random generator to conduct an
    long number = 23091991;
    RandomGenerator RG = new Well19937a(number);
    GaussianRandomGenerator NGG = new GaussianRandomGenerator(RG);
    UncorrelatedRandomVectorGenerator URVG = new UncorrelatedRandomVectorGenerator(mean, dAdIdRaan, NGG);
    double[][] rand_gen = new double[montecarlo_size][3];
    for (int jj = 0; jj < montecarlo_size; jj++) {
        rand_gen[jj] = URVG.nextVector();
    }
    // 
    FieldSpacecraftState<DerivativeStructure> SS_0 = new FieldSpacecraftState<>(KO);
    // adding force models
    ForceModel fModel_Sun = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
    ForceModel fModel_Moon = new ThirdBodyAttraction(CelestialBodyFactory.getMoon());
    ForceModel fModel_HFAM = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getNormalizedProvider(18, 18));
    // setting an hipparchus field integrator
    OrbitType type = OrbitType.CARTESIAN;
    double[][] tolerance = NumericalPropagator.tolerances(0.001, KO.toOrbit(), type);
    AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(zero.add(60));
    // setting of the field propagator, we used the numerical one in order to add the third body attraction
    // and the holmes featherstone force models
    FieldNumericalPropagator<DerivativeStructure> numProp = new FieldNumericalPropagator<>(field, integrator);
    numProp.setOrbitType(type);
    numProp.setInitialState(SS_0);
    numProp.addForceModel(fModel_Sun);
    numProp.addForceModel(fModel_Moon);
    numProp.addForceModel(fModel_HFAM);
    // with the master mode we will calulcate and print the error on every fixed step on the file error.txt
    // we defined the StepHandler to do that giving him the random number generator,
    // the size of the montecarlo simulation and the initial date
    numProp.setMasterMode(zero.add(int_step), new MyStepHandler<DerivativeStructure>(rand_gen, montecarlo_size, date_0, PW));
    // 
    long START = System.nanoTime();
    FieldSpacecraftState<DerivativeStructure> finalState = numProp.propagate(date_0.shiftedBy(final_Dt));
    long STOP = System.nanoTime();
    System.out.println((STOP - START) / 1E6 + " ms");
    System.out.println(finalState.getDate());
    PW.close();
}

Example 12 with ForceModel

use of org.orekit.forces.ForceModel in project Orekit by CS-SI.

the class PartialDerivativesEquations method computeDerivatives.

/**
 * {@inheritDoc}
 */
public double[] computeDerivatives(final SpacecraftState s, final double[] pDot) throws OrekitException {
    // initialize acceleration Jacobians to zero
    final int paramDim = selected.getNbParams();
    final int dim = 3;
    final double[][] dAccdParam = new double[dim][paramDim];
    final double[][] dAccdPos = new double[dim][dim];
    final double[][] dAccdVel = new double[dim][dim];
    final DSConverter fullConverter = new DSConverter(s, 6, propagator.getAttitudeProvider());
    final DSConverter posOnlyConverter = new DSConverter(s, 3, propagator.getAttitudeProvider());
    // compute acceleration Jacobians, finishing with the largest force: Newtonian attraction
    for (final ForceModel forceModel : propagator.getAllForceModels()) {
        final DSConverter converter = forceModel.dependsOnPositionOnly() ? posOnlyConverter : fullConverter;
        final FieldSpacecraftState<DerivativeStructure> dsState = converter.getState(forceModel);
        final DerivativeStructure[] parameters = converter.getParameters(dsState, forceModel);
        final FieldVector3D<DerivativeStructure> acceleration = forceModel.acceleration(dsState, parameters);
        final double[] derivativesX = acceleration.getX().getAllDerivatives();
        final double[] derivativesY = acceleration.getY().getAllDerivatives();
        final double[] derivativesZ = acceleration.getZ().getAllDerivatives();
        // update Jacobians with respect to state
        addToRow(derivativesX, 0, converter.getFreeStateParameters(), dAccdPos, dAccdVel);
        addToRow(derivativesY, 1, converter.getFreeStateParameters(), dAccdPos, dAccdVel);
        addToRow(derivativesZ, 2, converter.getFreeStateParameters(), dAccdPos, dAccdVel);
        int index = converter.getFreeStateParameters();
        for (ParameterDriver driver : forceModel.getParametersDrivers()) {
            if (driver.isSelected()) {
                final int parameterIndex = map.get(driver);
                ++index;
                dAccdParam[0][parameterIndex] += derivativesX[index];
                dAccdParam[1][parameterIndex] += derivativesY[index];
                dAccdParam[2][parameterIndex] += derivativesZ[index];
            }
        }
    }
    // the variational equations of the complete state Jacobian matrix have the following form:
    // [        |        ]   [                 |                  ]   [     |     ]
    // [  Adot  |  Bdot  ]   [  dVel/dPos = 0  |  dVel/dVel = Id  ]   [  A  |  B  ]
    // [        |        ]   [                 |                  ]   [     |     ]
    // ---------+---------   ------------------+------------------- * ------+------
    // [        |        ]   [                 |                  ]   [     |     ]
    // [  Cdot  |  Ddot  ] = [    dAcc/dPos    |     dAcc/dVel    ]   [  C  |  D  ]
    // [        |        ]   [                 |                  ]   [     |     ]
    // The A, B, C and D sub-matrices and their derivatives (Adot ...) are 3x3 matrices
    // The expanded multiplication above can be rewritten to take into account
    // the fixed values found in the sub-matrices in the left factor. This leads to:
    // [ Adot ] = [ C ]
    // [ Bdot ] = [ D ]
    // [ Cdot ] = [ dAcc/dPos ] * [ A ] + [ dAcc/dVel ] * [ C ]
    // [ Ddot ] = [ dAcc/dPos ] * [ B ] + [ dAcc/dVel ] * [ D ]
    // The following loops compute these expressions taking care of the mapping of the
    // (A, B, C, D) matrices into the single dimension array p and of the mapping of the
    // (Adot, Bdot, Cdot, Ddot) matrices into the single dimension array pDot.
    // copy C and E into Adot and Bdot
    final int stateDim = 6;
    final double[] p = s.getAdditionalState(getName());
    System.arraycopy(p, dim * stateDim, pDot, 0, dim * stateDim);
    // compute Cdot and Ddot
    for (int i = 0; i < dim; ++i) {
        final double[] dAdPi = dAccdPos[i];
        final double[] dAdVi = dAccdVel[i];
        for (int j = 0; j < stateDim; ++j) {
            pDot[(dim + i) * stateDim + j] = dAdPi[0] * p[j] + dAdPi[1] * p[j + stateDim] + dAdPi[2] * p[j + 2 * stateDim] + dAdVi[0] * p[j + 3 * stateDim] + dAdVi[1] * p[j + 4 * stateDim] + dAdVi[2] * p[j + 5 * stateDim];
        }
    }
    for (int k = 0; k < paramDim; ++k) {
        // the variational equations of the parameters Jacobian matrix are computed
        // one column at a time, they have the following form:
        // [      ]   [                 |                  ]   [   ]   [                  ]
        // [ Edot ]   [  dVel/dPos = 0  |  dVel/dVel = Id  ]   [ E ]   [  dVel/dParam = 0 ]
        // [      ]   [                 |                  ]   [   ]   [                  ]
        // --------   ------------------+------------------- * ----- + --------------------
        // [      ]   [                 |                  ]   [   ]   [                  ]
        // [ Fdot ] = [    dAcc/dPos    |     dAcc/dVel    ]   [ F ]   [    dAcc/dParam   ]
        // [      ]   [                 |                  ]   [   ]   [                  ]
        // The E and F sub-columns and their derivatives (Edot, Fdot) are 3 elements columns.
        // The expanded multiplication and addition above can be rewritten to take into
        // account the fixed values found in the sub-matrices in the left factor. This leads to:
        // [ Edot ] = [ F ]
        // [ Fdot ] = [ dAcc/dPos ] * [ E ] + [ dAcc/dVel ] * [ F ] + [ dAcc/dParam ]
        // The following loops compute these expressions taking care of the mapping of the
        // (E, F) columns into the single dimension array p and of the mapping of the
        // (Edot, Fdot) columns into the single dimension array pDot.
        // copy F into Edot
        final int columnTop = stateDim * stateDim + k;
        pDot[columnTop] = p[columnTop + 3 * paramDim];
        pDot[columnTop + paramDim] = p[columnTop + 4 * paramDim];
        pDot[columnTop + 2 * paramDim] = p[columnTop + 5 * paramDim];
        // compute Fdot
        for (int i = 0; i < dim; ++i) {
            final double[] dAdPi = dAccdPos[i];
            final double[] dAdVi = dAccdVel[i];
            pDot[columnTop + (dim + i) * paramDim] = dAccdParam[i][k] + dAdPi[0] * p[columnTop] + dAdPi[1] * p[columnTop + paramDim] + dAdPi[2] * p[columnTop + 2 * paramDim] + dAdVi[0] * p[columnTop + 3 * paramDim] + dAdVi[1] * p[columnTop + 4 * paramDim] + dAdVi[2] * p[columnTop + 5 * paramDim];
        }
    }
    // these equations have no effect on the main state itself
    return null;
}

Example 13 with ForceModel

use of org.orekit.forces.ForceModel in project Orekit by CS-SI.

the class SolidTidesTest method propagate.

private SpacecraftState propagate(Orbit orbit, AbsoluteDate target, ForceModel... forceModels) throws OrekitException {
    double[][] tolerances = NumericalPropagator.tolerances(10, orbit, OrbitType.KEPLERIAN);
    AbstractIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 300, tolerances[0], tolerances[1]);
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    for (ForceModel forceModel : forceModels) {
        propagator.addForceModel(forceModel);
    }
    propagator.setInitialState(new SpacecraftState(orbit));
    return propagator.propagate(target);
}

Example 14 with ForceModel

use of org.orekit.forces.ForceModel in project Orekit by CS-SI.

the class SolidTidesTest method testStateJacobianVs80ImplementationNoPoleTide.

@Test
public void testStateJacobianVs80ImplementationNoPoleTide() throws OrekitException {
    Frame eme2000 = FramesFactory.getEME2000();
    TimeScale utc = TimeScalesFactory.getUTC();
    AbsoluteDate date = new AbsoluteDate(2964, 8, 12, 11, 30, 00.000, utc);
    Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7), FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5), 0, PositionAngle.MEAN, eme2000, date, Constants.EIGEN5C_EARTH_MU);
    Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
    NormalizedSphericalHarmonicsProvider gravityField = GravityFieldFactory.getConstantNormalizedProvider(5, 5);
    ForceModel forceModel = new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(), gravityField.getTideSystem(), false, SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS, IERSConventions.IERS_2010, ut1, CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon());
    Assert.assertTrue(forceModel.dependsOnPositionOnly());
    checkStateJacobianVs80Implementation(new SpacecraftState(orbit), forceModel, new LofOffset(orbit.getFrame(), LOFType.VVLH), 2.0e-15, false);
}

Example 15 with ForceModel

use of org.orekit.forces.ForceModel in project Orekit by CS-SI.

the class SolidTidesTest method testStateJacobianVsFiniteDifferencesNoPoleTide.

@Test
public void testStateJacobianVsFiniteDifferencesNoPoleTide() throws OrekitException {
    Frame eme2000 = FramesFactory.getEME2000();
    TimeScale utc = TimeScalesFactory.getUTC();
    AbsoluteDate date = new AbsoluteDate(2964, 8, 12, 11, 30, 00.000, utc);
    Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7), FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5), 0, PositionAngle.MEAN, eme2000, date, Constants.EIGEN5C_EARTH_MU);
    Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
    NormalizedSphericalHarmonicsProvider gravityField = GravityFieldFactory.getConstantNormalizedProvider(5, 5);
    ForceModel forceModel = new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(), gravityField.getTideSystem(), false, SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS, IERSConventions.IERS_2010, ut1, CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon());
    checkStateJacobianVsFiniteDifferences(new SpacecraftState(orbit), forceModel, Propagator.DEFAULT_LAW, 10.0, 2.0e-10, false);
}