1 /* Copyright 2002-2021 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.propagation.SpacecraftState;
23 import org.orekit.propagation.integration.AdditionalEquations;
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 /** Number of patch points. */
36 private int npoints;
37
38 /** Simple Constructor.
39 * <p> Standard constructor for multiple shooting which can be used with the CR3BP model.</p>
40 * @param initialGuessList initial patch points to be corrected.
41 * @param propagatorList list of propagators associated to each patch point.
42 * @param additionalEquations list of additional equations linked to propagatorList.
43 * @param arcDuration initial guess of the duration of each arc.
44 * @param tolerance convergence tolerance on the constraint vector
45 */
46 public CR3BPMultipleShooter(final List<SpacecraftState> initialGuessList, final List<NumericalPropagator> propagatorList,
47 final List<AdditionalEquations> additionalEquations, final double arcDuration, final double tolerance) {
48 super(initialGuessList, propagatorList, additionalEquations, arcDuration, tolerance, "stmEquations");
49 this.npoints = initialGuessList.size();
50 }
51
52 /** {@inheritDoc} */
53 protected SpacecraftState getAugmentedInitialState(final SpacecraftState initialState,
54 final AdditionalEquations additionalEquation) {
55 return ((STMEquations) additionalEquation).setInitialPhi(initialState);
56 }
57
58 /** {@inheritDoc} */
59 protected double[][] computeAdditionalJacobianMatrix(final List<SpacecraftState> propagatedSP) {
60
61 final Map<Integer, Double> mapConstraints = getConstraintsMap();
62
63 final boolean isClosedOrbit = isClosedOrbit();
64
65 // Number of additional constraints
66 final int n = mapConstraints.size() + (isClosedOrbit ? 6 : 0);
67
68 final int ncolumns = getNumberOfFreeVariables() - 1;
69
70 final double[][] M = new double[n][ncolumns];
71
72 int k = 0;
73 if (isClosedOrbit) {
74 // The Jacobian matrix has the following form:
75 //
76 // [-1 0 0 ... 1 0 ]
77 // [ 0 -1 0 0 ... 1 0 ]
78 // C = [ 0 -1 0 0 ... 1 0 ]
79 // [ 0 -1 0 0 ... 1 0 ]
80 // [ 0 -1 0 0 ... 1 0 ]
81 // [ 0 0 -1 0 ... 0 1 ]
82
83 for (int i = 0; i < 6; i++) {
84 M[i][i] = -1;
85 M[i][ncolumns - 6 + i] = 1;
86 }
87 k = 6;
88 }
89
90 for (int index : mapConstraints.keySet()) {
91 M[k][index] = 1;
92 k++;
93 }
94 return M;
95 }
96
97 /** {@inheritDoc} */
98 @Override
99 protected double[][] computeEpochJacobianMatrix(final List<SpacecraftState> propagatedSP) {
100 final int nFreeEpoch = getNumberOfFreeEpoch();
101 // Rows and columns dimensions
102 final int ncolumns = 1 + nFreeEpoch;
103 final int nrows = npoints - 1;
104 // Return an empty array
105 return new double[nrows][ncolumns];
106 }
107
108 /** {@inheritDoc} */
109 protected double[] computeAdditionalConstraints(final List<SpacecraftState> propagatedSP) {
110
111 // The additional constraint vector has the following form :
112
113 // [ xni - x1i ]----
114 // [ yni - x1i ] |
115 // Fadd(X) = [ zni - x1i ] vector's component
116 // [vxni - vx1i] for a closed orbit
117 // [vyni - vy1i] |
118 // [vzni - vz1i]----
119 // [ y1i - y1d ]---- other constraints (component of
120 // [ ... ] | a patch point eaquals to a
121 // [vz2i - vz2d]---- desired value)
122
123 final Map<Integer, Double> mapConstraints = getConstraintsMap();
124 final boolean isClosedOrbit = isClosedOrbit();
125 // Number of additional constraints
126 final int n = mapConstraints.size() + (isClosedOrbit ? 6 : 0);
127
128 final List<SpacecraftState> patchedSpacecraftStates = getPatchedSpacecraftState();
129
130 final double[] fxAdditionnal = new double[n];
131 int i = 0;
132
133 if (isClosedOrbit) {
134
135 final AbsolutePVCoordinates apv1i = patchedSpacecraftStates.get(0).getAbsPVA();
136 final AbsolutePVCoordinates apvni = patchedSpacecraftStates.get(npoints - 1).getAbsPVA();
137
138 fxAdditionnal[0] = apvni.getPosition().getX() - apv1i.getPosition().getX();
139 fxAdditionnal[1] = apvni.getPosition().getY() - apv1i.getPosition().getY();
140 fxAdditionnal[2] = apvni.getPosition().getZ() - apv1i.getPosition().getZ();
141 fxAdditionnal[3] = apvni.getVelocity().getX() - apv1i.getVelocity().getX();
142 fxAdditionnal[4] = apvni.getVelocity().getY() - apv1i.getVelocity().getY();
143 fxAdditionnal[5] = apvni.getVelocity().getZ() - apv1i.getVelocity().getZ();
144
145 i = 6;
146 }
147
148 // Update additional constraints
149 updateAdditionalConstraints(i, fxAdditionnal);
150 return fxAdditionnal;
151 }
152
153 }