CR3BPMultipleShooter.java

  1. /* Copyright 2002-2020 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.propagation.numerical.cr3bp;

  18. import java.util.List;
  19. import java.util.Map;

  20. import org.orekit.propagation.SpacecraftState;
  21. import org.orekit.propagation.integration.AdditionalEquations;
  22. import org.orekit.propagation.numerical.NumericalPropagator;
  23. import org.orekit.utils.AbsolutePVCoordinates;
  24. import org.orekit.utils.AbstractMultipleShooting;

  25. /**
  26.  * Multiple shooting method applicable for orbits, either propagation in CR3BP, or in an ephemeris model.
  27.  * @see "TRAJECTORY DESIGN AND ORBIT MAINTENANCE STRATEGIES IN MULTI-BODY DYNAMICAL REGIMES by Thomas A. Pavlak, Purdue University"
  28.  * @author William Desprats
  29.  */
  30. public class CR3BPMultipleShooter extends AbstractMultipleShooting {

  31.     /** Number of patch points. */
  32.     private int npoints;

  33.     /** Simple Constructor.
  34.      * <p> Standard constructor for multiple shooting which can be used with the CR3BP model.</p>
  35.      * @param initialGuessList initial patch points to be corrected.
  36.      * @param propagatorList list of propagators associated to each patch point.
  37.      * @param additionalEquations list of additional equations linked to propagatorList.
  38.      * @param arcDuration initial guess of the duration of each arc.
  39.      * @param tolerance convergence tolerance on the constraint vector
  40.      */
  41.     public CR3BPMultipleShooter(final List<SpacecraftState> initialGuessList, final List<NumericalPropagator> propagatorList,
  42.                                  final List<AdditionalEquations> additionalEquations, final double arcDuration, final double tolerance) {
  43.         super(initialGuessList, propagatorList, additionalEquations, arcDuration, tolerance, "stmEquations");
  44.         this.npoints = initialGuessList.size();
  45.     }

  46.     /** {@inheritDoc} */
  47.     protected SpacecraftState getAugmentedInitialState(final SpacecraftState initialState,
  48.                                                        final AdditionalEquations additionalEquation) {
  49.         return ((STMEquations) additionalEquation).setInitialPhi(initialState);
  50.     }

  51.     /** {@inheritDoc} */
  52.     protected double[][] computeAdditionalJacobianMatrix(final List<SpacecraftState> propagatedSP) {

  53.         final Map<Integer, Double> mapConstraints = getConstraintsMap();

  54.         final boolean isClosedOrbit = isClosedOrbit();

  55.         // Number of additional constraints
  56.         final int n = mapConstraints.size() + (isClosedOrbit ? 6 : 0);

  57.         final int ncolumns = getNumberOfFreeVariables() - 1;

  58.         final double[][] M = new double[n][ncolumns];

  59.         int k = 0;
  60.         if (isClosedOrbit) {
  61.             // The Jacobian matrix has the following form:
  62.             //
  63.             //      [-1  0              0  ...  1  0             ]
  64.             //      [ 0 -1  0           0  ...     1  0          ]
  65.             // C =  [    0 -1  0        0  ...        1  0       ]
  66.             //      [       0 -1  0     0  ...           1  0    ]
  67.             //      [          0  -1 0  0  ...              1  0 ]
  68.             //      [          0  0 -1  0  ...              0  1 ]

  69.             for (int i = 0; i < 6; i++) {
  70.                 M[i][i] = -1;
  71.                 M[i][ncolumns - 6 + i] = 1;
  72.             }
  73.             k = 6;
  74.         }

  75.         for (int index : mapConstraints.keySet()) {
  76.             M[k][index] = 1;
  77.             k++;
  78.         }
  79.         return M;
  80.     }

  81.     /** {@inheritDoc} */
  82.     @Override
  83.     protected double[][] computeEpochJacobianMatrix(final List<SpacecraftState> propagatedSP) {
  84.         final int nFreeEpoch = getNumberOfFreeEpoch();
  85.         // Rows and columns dimensions
  86.         final int ncolumns   = 1 + nFreeEpoch;
  87.         final int nrows      = npoints - 1;
  88.         // Return an empty array
  89.         return new double[nrows][ncolumns];
  90.     }

  91.     /** {@inheritDoc} */
  92.     protected double[] computeAdditionalConstraints(final List<SpacecraftState> propagatedSP) {

  93.         // The additional constraint vector has the following form :

  94.         //           [ xni - x1i ]----
  95.         //           [ yni - x1i ]    |
  96.         // Fadd(X) = [ zni - x1i ] vector's component
  97.         //           [vxni - vx1i] for a closed orbit
  98.         //           [vyni - vy1i]    |
  99.         //           [vzni - vz1i]----
  100.         //           [ y1i - y1d ]---- other constraints (component of
  101.         //           [    ...    ]    | a patch point eaquals to a
  102.         //           [vz2i - vz2d]----  desired value)

  103.         final Map<Integer, Double> mapConstraints = getConstraintsMap();
  104.         final boolean isClosedOrbit = isClosedOrbit();
  105.         // Number of additional constraints
  106.         final int n = mapConstraints.size() + (isClosedOrbit ? 6 : 0);

  107.         final List<SpacecraftState> patchedSpacecraftStates = getPatchedSpacecraftState();

  108.         final double[] fxAdditionnal = new double[n];
  109.         int i = 0;

  110.         if (isClosedOrbit) {

  111.             final AbsolutePVCoordinates apv1i = patchedSpacecraftStates.get(0).getAbsPVA();
  112.             final AbsolutePVCoordinates apvni = patchedSpacecraftStates.get(npoints - 1).getAbsPVA();

  113.             fxAdditionnal[0] = apvni.getPosition().getX() - apv1i.getPosition().getX();
  114.             fxAdditionnal[1] = apvni.getPosition().getY() - apv1i.getPosition().getY();
  115.             fxAdditionnal[2] = apvni.getPosition().getZ() - apv1i.getPosition().getZ();
  116.             fxAdditionnal[3] = apvni.getVelocity().getX() - apv1i.getVelocity().getX();
  117.             fxAdditionnal[4] = apvni.getVelocity().getY() - apv1i.getVelocity().getY();
  118.             fxAdditionnal[5] = apvni.getVelocity().getZ() - apv1i.getVelocity().getZ();

  119.             i = 6;
  120.         }

  121.         // Update additional constraints
  122.         updateAdditionalConstraints(i, fxAdditionnal);
  123.         return fxAdditionnal;
  124.     }

  125. }