1   /* Copyright 2002-2025 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  
19  import java.util.List;
20  import java.util.Map;
21  
22  import org.orekit.errors.OrekitException;
23  import org.orekit.errors.OrekitMessages;
24  import org.orekit.propagation.SpacecraftState;
25  import org.orekit.propagation.numerical.NumericalPropagator;
26  import org.orekit.utils.AbsolutePVCoordinates;
27  import org.orekit.utils.AbstractMultipleShooting;
28  
29  /**
30   * Multiple shooting method applicable for orbits, either propagation in CR3BP, or in an ephemeris model.
31   * @see "TRAJECTORY DESIGN AND ORBIT MAINTENANCE STRATEGIES IN MULTI-BODY DYNAMICAL REGIMES by Thomas A. Pavlak, Purdue University"
32   * @author William Desprats
33   * @author Alberto Fossà
34   */
35  public class CR3BPMultipleShooter extends AbstractMultipleShooting {
36  
37      /** Name of the derivatives. */
38      private static final String STM = "stmEquations";
39  
40      /** Derivatives linked to the Propagators.
41       * @since 11.1
42       */
43      private final List<STMEquations> stmEquations;
44  
45      /** True if orbit is closed. */
46      private boolean isClosedOrbit;
47  
48      /** Simple Constructor.
49       * <p> Standard constructor for multiple shooting which can be used with the CR3BP model. </p>
50       * @param initialGuessList initial patch points to be corrected
51       * @param propagatorList list of propagators associated to each patch point
52       * @param stmEquations list of additional derivatives providers linked to propagatorList
53       * @param tolerance convergence tolerance on the constraint vector
54       * @param maxIter maximum number of iterations
55       */
56      public CR3BPMultipleShooter(final List<SpacecraftState> initialGuessList,
57                                  final List<NumericalPropagator> propagatorList,
58                                  final List<STMEquations> stmEquations,
59                                  final double tolerance, final int maxIter) {
60          super(initialGuessList, propagatorList, tolerance, maxIter, true, STM);
61          this.stmEquations = stmEquations;
62      }
63  
64      /** {@inheritDoc} */
65      protected SpacecraftState getAugmentedInitialState(final int i) {
66          return stmEquations.get(i).setInitialPhi(getPatchPoint(i));
67      }
68  
69      /** {@inheritDoc} */
70      @Override
71      public void setEpochFreedom(final int patchIndex, final boolean isFree) {
72          throw new OrekitException(OrekitMessages.FUNCTION_NOT_IMPLEMENTED);
73      }
74  
75      /** {@inheritDoc} */
76      @Override
77      public void setScaleLength(final double scaleLength) {
78          throw new OrekitException(OrekitMessages.FUNCTION_NOT_IMPLEMENTED);
79      }
80  
81      /** {@inheritDoc} */
82      @Override
83      public void setScaleTime(final double scaleTime) {
84          throw new OrekitException(OrekitMessages.FUNCTION_NOT_IMPLEMENTED);
85      }
86  
87      /** Set the constraint of a closed orbit or not.
88       *  @param isClosed true if orbit should be closed
89       */
90      public void setClosedOrbitConstraint(final boolean isClosed) {
91          this.isClosedOrbit = isClosed;
92      }
93  
94      /** {@inheritDoc} */
95      @Override
96      protected double[][] computeAdditionalJacobianMatrix(final List<SpacecraftState> propagatedSP) {
97  
98          final Map<Integer, Double> mapConstraints = getConstraintsMap();
99          final boolean[] freeCompsMap              = getFreeCompsMap();
100 
101         // Number of additional constraints
102         final int nRows = mapConstraints.size() + (isClosedOrbit ? 6 : 0);
103         final int nCols = getNumberOfFreeComponents();
104 
105         final double[][] M = new double[nRows][nCols];
106 
107         int j = 0;
108         if (isClosedOrbit) {
109             // The Jacobian matrix has the following form:
110             //
111             //      [-1  0              0  ...  1  0             ]
112             //      [ 0 -1  0           0  ...     1  0          ]
113             // F =  [    0 -1  0        0  ...        1  0       ]
114             //      [       0 -1  0     0  ...           1  0    ]
115             //      [          0 -1  0  0  ...              1  0 ]
116             //      [          0  0 -1  0  ...              0  1 ]
117 
118             int index = 0;
119             for (int i = 0; i < 6; i++) {
120                 if (freeCompsMap[i]) {
121                     M[i][index] = -1.0;
122                     index++;
123                 }
124             }
125             index = nCols - 6;
126             for (int i = 0; i < 6; i++) {
127                 if (freeCompsMap[nCols - 6 + i]) {
128                     M[i][index] = 1.0;
129                     index++;
130                 }
131             }
132             j = 6;
133         }
134 
135         for (int k : mapConstraints.keySet()) {
136             M[j][k] = 1.0;
137             j++;
138         }
139 
140         return M;
141     }
142 
143     /** {@inheritDoc} */
144     @Override
145     protected double[] computeAdditionalConstraints(final List<SpacecraftState> propagatedSP) {
146 
147         // The additional constraint vector has the following form :
148 
149         //           [ xni - x1i ]----
150         //           [ yni - x1i ]    |
151         // Fadd(X) = [ zni - x1i ] vector's component
152         //           [vxni - vx1i] for a closed orbit
153         //           [vyni - vy1i]    |
154         //           [vzni - vz1i]----
155         //           [ y1i - y1d ]---- other constraints (component of
156         //           [    ...    ]    | a patch point equals to a
157         //           [vz2i - vz2d]----  desired value)
158 
159         final Map<Integer, Double> mapConstraints           = getConstraintsMap();
160         final List<SpacecraftState> patchedSpacecraftStates = getPatchedSpacecraftState();
161 
162         final double[] fxAdditional = new double[mapConstraints.size() + (isClosedOrbit ? 6 : 0)];
163         int i = 0;
164 
165         if (isClosedOrbit) {
166 
167             final AbsolutePVCoordinates apv1i = patchedSpacecraftStates.get(0).getAbsPVA();
168             final AbsolutePVCoordinates apvni = patchedSpacecraftStates.get(patchedSpacecraftStates.size() - 1).getAbsPVA();
169 
170             fxAdditional[0] = apvni.getPosition().getX() - apv1i.getPosition().getX();
171             fxAdditional[1] = apvni.getPosition().getY() - apv1i.getPosition().getY();
172             fxAdditional[2] = apvni.getPosition().getZ() - apv1i.getPosition().getZ();
173             fxAdditional[3] = apvni.getVelocity().getX() - apv1i.getVelocity().getX();
174             fxAdditional[4] = apvni.getVelocity().getY() - apv1i.getVelocity().getY();
175             fxAdditional[5] = apvni.getVelocity().getZ() - apv1i.getVelocity().getZ();
176 
177             i = 6;
178         }
179 
180         // Update additional constraints
181         updateAdditionalConstraints(i, fxAdditional);
182         return fxAdditional;
183     }
184 
185     /** {@inheritDoc} */
186     @Override
187     protected int getNumberOfConstraints() {
188         return super.getNumberOfConstraints() + (isClosedOrbit ? 6 : 0);
189     }
190 
191 }