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