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.utils;
18
19 import java.util.ArrayList;
20 import java.util.Arrays;
21 import java.util.HashMap;
22 import java.util.List;
23 import java.util.Map;
24
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.geometry.euclidean.threed.Vector3D;
27 import org.hipparchus.linear.LUDecomposition;
28 import org.hipparchus.linear.MatrixUtils;
29 import org.hipparchus.linear.QRDecomposition;
30 import org.hipparchus.linear.RealMatrix;
31 import org.hipparchus.linear.RealVector;
32 import org.orekit.attitudes.Attitude;
33 import org.orekit.attitudes.AttitudeProvider;
34 import org.orekit.propagation.SpacecraftState;
35 import org.orekit.propagation.numerical.NumericalPropagator;
36 import org.orekit.time.AbsoluteDate;
37
38 /**
39 * Multiple shooting method using only constraints on state vectors of patch points (and possibly on epoch and integration time).
40 * @see "TRAJECTORY DESIGN AND ORBIT MAINTENANCE STRATEGIES IN MULTI-BODY DYNAMICAL REGIMES by Thomas A. Pavlak, Purdue University"
41 * @author William Desprats
42 * @author Alberto Fossà
43 * @since 10.2
44 */
45 public abstract class AbstractMultipleShooting implements MultipleShooting {
46
47 /** Patch points along the trajectory. */
48 private final List<SpacecraftState> patchedSpacecraftStates;
49
50 /** List of Propagators. */
51 private final List<NumericalPropagator> propagatorList;
52
53 /** Duration of propagation along each arc. */
54 private final double[] propagationTime;
55
56 /** Free components of patch points. */
57 private final boolean[] freeCompsMap;
58
59 /** Free epochs of patch points. */
60 private final boolean[] freeEpochMap;
61
62 /** Number of free state components. */
63 private int nComps;
64
65 /** Number of free arc duration. */
66 // TODO add possibility to fix integration time span?
67 private final int nDuration;
68
69 /** Number of free epochs. */
70 private int nEpoch;
71
72 /** Scale time for update computation. */
73 private double scaleTime;
74
75 /** Scale length for update computation. */
76 private double scaleLength;
77
78 /** Patch points components which are constrained. */
79 private final Map<Integer, Double> mapConstraints;
80
81 /** True if the dynamical system is autonomous.
82 * In this case epochs and epochs constraints are omitted from the problem formulation. */
83 private final boolean isAutonomous;
84
85 /** Tolerance on the constraint vector. */
86 private final double tolerance;
87
88 /** Maximum number of iterations. */
89 private final int maxIter;
90
91 /** Expected name of the additional equations. */
92 private final String additionalName;
93
94 /** Simple Constructor.
95 * <p> Standard constructor for multiple shooting </p>
96 * @param initialGuessList initial patch points to be corrected
97 * @param propagatorList list of propagators associated to each patch point
98 * @param tolerance convergence tolerance on the constraint vector
99 * @param maxIter maximum number of iterations
100 * @param isAutonomous true if the dynamical system is autonomous (i.e. not dependent on the epoch)
101 * @param additionalName name of the additional equations
102 * @since 11.1
103 */
104 protected AbstractMultipleShooting(final List<SpacecraftState> initialGuessList,
105 final List<NumericalPropagator> propagatorList,
106 final double tolerance, final int maxIter,
107 final boolean isAutonomous, final String additionalName) {
108
109 this.patchedSpacecraftStates = initialGuessList;
110 this.propagatorList = propagatorList;
111 this.isAutonomous = isAutonomous;
112 this.additionalName = additionalName;
113
114 // propagation duration
115 final int propagationNumber = initialGuessList.size() - 1;
116 propagationTime = new double[propagationNumber];
117 for (int i = 0; i < propagationNumber; i++) {
118 propagationTime[i] = initialGuessList.get(i + 1).getDate().durationFrom(initialGuessList.get(i).getDate());
119 }
120
121 // states components freedom
122 this.freeCompsMap = new boolean[6 * initialGuessList.size()];
123 Arrays.fill(freeCompsMap, true);
124
125 // epochs freedom
126 if (isAutonomous) {
127 // epochs omitted from problem formulation
128 this.freeEpochMap = new boolean[0];
129 } else {
130 this.freeEpochMap = new boolean[initialGuessList.size()];
131 Arrays.fill(freeEpochMap, true);
132 }
133
134 // number of free variables
135 this.nComps = 6 * initialGuessList.size();
136 this.nDuration = propagationNumber;
137 this.nEpoch = freeEpochMap.length;
138
139 // convergence criteria
140 this.tolerance = tolerance;
141 this.maxIter = maxIter;
142
143 // scaling
144 this.scaleTime = 1.0;
145 this.scaleLength = 1.0;
146
147 // all additional constraints must be set afterward
148 this.mapConstraints = new HashMap<>();
149 }
150
151 /** Get a patch point.
152 * @param i index of the patch point
153 * @return state of the patch point
154 * @since 11.1
155 */
156 protected SpacecraftState getPatchPoint(final int i) {
157 return patchedSpacecraftStates.get(i);
158 }
159
160 /** Set a component of a patch point to free or not.
161 * @param patchIndex Patch point index (zero-based)
162 * @param componentIndex Index of the component to be constrained (zero-based)
163 * @param isFree constraint value
164 */
165 public void setPatchPointComponentFreedom(final int patchIndex, final int componentIndex, final boolean isFree) {
166 if (freeCompsMap[6 * patchIndex + componentIndex] != isFree) {
167 freeCompsMap[6 * patchIndex + componentIndex] = isFree;
168 nComps += isFree ? 1 : -1;
169 }
170 }
171
172 /** Set the epoch of a patch point to free or not.
173 * @param patchIndex Patch point index (zero-based)
174 * @param isFree constraint value
175 */
176 public void setEpochFreedom(final int patchIndex, final boolean isFree) {
177 if (freeEpochMap[patchIndex] != isFree) {
178 freeEpochMap[patchIndex] = isFree;
179 nEpoch += isFree ? 1 : -1;
180 }
181 }
182
183 /** Set the scale time.
184 * @param scaleTime scale time in seconds
185 */
186 public void setScaleTime(final double scaleTime) {
187 this.scaleTime = scaleTime;
188 }
189
190 /** Set the scale length.
191 * @param scaleLength scale length in meters
192 */
193 public void setScaleLength(final double scaleLength) {
194 this.scaleLength = scaleLength;
195 }
196
197 /** Add a constraint on one component of one patch point.
198 * @param patchIndex Patch point index (zero-based)
199 * @param componentIndex Index of the component which is constrained (zero-based)
200 * @param constraintValue constraint value
201 */
202 public void addConstraint(final int patchIndex, final int componentIndex, final double constraintValue) {
203 mapConstraints.put(patchIndex * 6 + componentIndex, constraintValue);
204 }
205
206 /** {@inheritDoc} */
207 @Override
208 public List<SpacecraftState> compute() {
209
210 int iter = 0; // number of iteration
211 double fxNorm;
212
213 while (iter < maxIter) {
214
215 iter++;
216
217 // propagation (multi-threading see PropagatorsParallelizer)
218 final List<SpacecraftState> propagatedSP = propagatePatchedSpacecraftState();
219
220 // constraints computation
221 final RealVector fx = MatrixUtils.createRealVector(computeConstraint(propagatedSP));
222
223 // convergence check
224 fxNorm = fx.getNorm();
225 // System.out.printf(Locale.US, "Iter: %3d Error: %.16e%n", iter, fxNorm);
226 if (fxNorm < tolerance) {
227 break;
228 }
229
230 // Jacobian matrix computation
231 final RealMatrix M = computeJacobianMatrix(propagatedSP);
232
233 // correction computation using minimum norm approach
234 // (i.e. minimize difference between solutions from successive iterations,
235 // in other word try to stay close to initial guess. This is *not* a least-squares solution)
236 // see equation 3.12 in Pavlak's thesis
237 final RealMatrix MMt = M.multiplyTransposed(M);
238 RealVector dx;
239 try {
240 dx = M.transpose().operate(new LUDecomposition(MMt, 0.0).getSolver().solve(fx));
241 } catch (MathIllegalArgumentException e) {
242 dx = M.transpose().operate(new QRDecomposition(MMt, 0.0).getSolver().solve(fx));
243 }
244
245 // trajectory update
246 updateTrajectory(dx);
247 }
248
249 return patchedSpacecraftStates;
250 }
251
252 /** Compute the Jacobian matrix of the problem.
253 * @param propagatedSP List of propagated SpacecraftStates (patch points)
254 * @return jacobianMatrix Jacobian matrix
255 */
256 private RealMatrix computeJacobianMatrix(final List<SpacecraftState> propagatedSP) {
257
258 // The Jacobian matrix has the following form:
259 //
260 // [ | | ]
261 // [ A | B | C ]
262 // DF(X) = [ | | ]
263 // [-----------------]
264 // [ 0 | D | E ]
265 // [-----------------]
266 // [ F | 0 | 0 ]
267 //
268 // For a problem in which all the components of each patch points are free, A is detailed below:
269 //
270 // [ phi1 -I ]
271 // [ phi2 -I ]
272 // A = [ .... .... ] 6(n-1)x6n
273 // [ .... .... ]
274 // [ phin-1 -I ]
275 //
276 // If the duration of the propagation of each arc is the same:
277 //
278 // [ xdot1f ]
279 // [ xdot2f ] [ -1 ]
280 // B = [ .... ] 6(n-1)x1 and D = [ ... ]
281 // [ .... ] [ -1 ]
282 // [xdotn-1f]
283 //
284 // Otherwise:
285 //
286 // [ xdot1f ]
287 // [ xdot2f ]
288 // B = [ .... ] 6(n-1)x(n-1) and D = -I
289 // [ .... ]
290 // [ xdotn-1f ]
291 //
292 // If the problem is not dependant on the epoch (e.g. CR3BP), the C, D and E matrices are not computed.
293 //
294 // Otherwise:
295 // [ dx1f/dtau1 0 ]
296 // [ dx2f/dtau2 0 ]
297 // C = [ .... 0 ] 6(n-1)xn
298 // [ .... 0 ]
299 // [ dxn-1f/dtaun-1 0 ]
300 //
301 // [ -1 1 0 ]
302 // [ -1 1 0 ]
303 // E = [ .. .. ] n-1xn
304 // [ .. .. ]
305 // [ -1 1 ]
306 //
307 // F is computed according to additional constraints
308 // (for now, closed orbit, or a component of a patch point equals to a specified value)
309
310 final int nArcs = patchedSpacecraftStates.size() - 1;
311 final double scaleVel = scaleLength / scaleTime;
312 final double scaleAcc = scaleVel / scaleTime;
313 final RealMatrix M = MatrixUtils.createRealMatrix(getNumberOfConstraints(), getNumberOfFreeVariables());
314
315 int index = 0; // first column index for matrix A
316 int indexDuration = nComps; // first column index for matrix B
317 int indexEpoch = indexDuration + nDuration; // first column index for matrix C
318 for (int i = 0; i < nArcs; i++) {
319
320 final SpacecraftState finalState = propagatedSP.get(i);
321
322 // PV coordinates and state transition matrix at final time
323 final PVCoordinates pvf = finalState.getPVCoordinates();
324 final double[][] phi = getStateTransitionMatrix(finalState); // already scaled
325
326 // Matrix A
327 for (int j = 0; j < 6; j++) { // Loop on 6 component of the patch point
328 if (freeCompsMap[6 * i + j]) { // If this component is free
329 for (int k = 0; k < 6; k++) { // Loop on the 6 component of the patch point constraint
330 M.setEntry(6 * i + k, index, phi[k][j]);
331 }
332 if (i > 0) {
333 M.setEntry(6 * (i - 1) + j, index, -1.0);
334 }
335 index++;
336 }
337 }
338
339 // Matrix B
340 final double[][] pvfArray = new double[][] {
341 {pvf.getVelocity().getX() / scaleVel},
342 {pvf.getVelocity().getY() / scaleVel},
343 {pvf.getVelocity().getZ() / scaleVel},
344 {pvf.getAcceleration().getX() / scaleAcc},
345 {pvf.getAcceleration().getY() / scaleAcc},
346 {pvf.getAcceleration().getZ() / scaleAcc}};
347
348 M.setSubMatrix(pvfArray, 6 * i, indexDuration);
349 indexDuration++;
350
351 // Matrix C
352 // there is a typo in Pavlak's thesis, equations 3.48-3.49:
353 // the sign in front of the partials of the states with respect to epochs should be plus
354 if (!isAutonomous) {
355 if (freeEpochMap[i]) { // If this component is free
356 final double[] derivatives = finalState.getAdditionalState(additionalName);
357 final double[][] subC = new double[6][1];
358 for (int j = 0; j < 3; j++) {
359 subC[j][0] = derivatives[derivatives.length - 6 + j] / scaleVel;
360 subC[j + 3][0] = derivatives[derivatives.length - 3 + j] / scaleAcc;
361 }
362 M.setSubMatrix(subC, 6 * i, indexEpoch);
363 indexEpoch++;
364 }
365 }
366 }
367
368 // complete Matrix A
369 for (int j = 0; j < 6; j++) { // Loop on 6 component of the patch point
370 if (freeCompsMap[6 * nArcs + j]) { // If this component is free
371 M.setEntry(6 * (nArcs - 1) + j, index, -1.0);
372 index++;
373 }
374 }
375
376 // Matrices D and E
377 if (!isAutonomous) {
378 final double[][] subDE = computeEpochJacobianMatrix(propagatedSP);
379 M.setSubMatrix(subDE, 6 * nArcs, nComps);
380 }
381
382 // Matrix F
383 final double[][] subF = computeAdditionalJacobianMatrix(propagatedSP);
384 if (subF.length > 0) {
385 M.setSubMatrix(subF, isAutonomous ? 6 * nArcs : 7 * nArcs, 0);
386 }
387
388 return M;
389 }
390
391 /** Compute the constraint vector of the problem.
392 * @param propagatedSP propagated SpacecraftState
393 * @return constraint vector
394 */
395 private double[] computeConstraint(final List<SpacecraftState> propagatedSP) {
396
397 // The Constraint vector has the following form :
398 //
399 // [ x1f - x2i ]---
400 // [ x2f - x3i ] |
401 // F(X) = [ ... ] vectors' equality for a continuous trajectory
402 // [ ... ] |
403 // [xn-1f - xni]---
404 // [ d2-(d1+T) ] continuity between epoch
405 // [ ... ] and integration time
406 // [dn-(dn-1+T)]---
407 // [ ... ] additional
408 // [ ... ] constraints
409
410 final int nPoints = patchedSpacecraftStates.size();
411 final double scaleVel = scaleLength / scaleTime;
412 final double[] fx = new double[getNumberOfConstraints()];
413
414 // state continuity
415 for (int i = 0; i < nPoints - 1; i++) {
416 final AbsolutePVCoordinates absPvi = patchedSpacecraftStates.get(i + 1).getAbsPVA();
417 final AbsolutePVCoordinates absPvf = propagatedSP.get(i).getAbsPVA();
418 final double[] deltaPos = absPvf.getPosition().subtract(absPvi.getPosition()).toArray();
419 final double[] deltaVel = absPvf.getVelocity().subtract(absPvi.getVelocity()).toArray();
420 for (int j = 0; j < 3; j++) {
421 fx[6 * i + j] = deltaPos[j] / scaleLength;
422 fx[6 * i + 3 + j] = deltaVel[j] / scaleVel;
423 }
424 }
425
426 int index = 6 * (nPoints - 1);
427
428 // epoch constraints
429 if (!isAutonomous) {
430 for (int i = 0; i < nPoints - 1; i++) {
431 final double deltaEpoch = patchedSpacecraftStates.get(i + 1).getDate()
432 .durationFrom(patchedSpacecraftStates.get(i).getDate());
433 fx[index] = (deltaEpoch - propagationTime[i]) / scaleTime;
434 index++;
435 }
436 }
437
438 // additional constraints
439 final double[] additionalConstraints = computeAdditionalConstraints(propagatedSP);
440 for (double constraint : additionalConstraints) {
441 fx[index] = constraint;
442 index++;
443 }
444
445 return fx;
446 }
447
448 /** Update the trajectory, and the propagation time.
449 * @param dx correction on the initial vector
450 */
451 private void updateTrajectory(final RealVector dx) {
452 // X = [x1, ..., xn, T1, ..., Tn, d1, ..., dn]
453
454 final double scaleVel = scaleLength / scaleTime;
455
456 // Update propagation time
457 int indexDuration = nComps;
458 for (int i = 0; i < nDuration; i++) {
459 propagationTime[i] -= dx.getEntry(indexDuration) * scaleTime;
460 indexDuration++;
461 }
462
463 // Update patch points through SpacecraftStates
464 int index = 0;
465 int indexEpoch = nComps + nDuration;
466
467 for (int i = 0; i < patchedSpacecraftStates.size(); i++) {
468
469 // Get delta in position and velocity
470 final double[] deltaPV = new double[6];
471 for (int j = 0; j < 6; j++) { // Loop on 6 component of the patch point
472 if (freeCompsMap[6 * i + j]) { // If this component is free (is to be updated)
473 deltaPV[j] = dx.getEntry(index);
474 index++;
475 }
476 }
477 final Vector3D deltaP = new Vector3D(deltaPV[0], deltaPV[1], deltaPV[2]).scalarMultiply(scaleLength);
478 final Vector3D deltaV = new Vector3D(deltaPV[3], deltaPV[4], deltaPV[5]).scalarMultiply(scaleVel);
479
480 // Update the PVCoordinates of the patch point
481 final AbsolutePVCoordinates currentAPV = patchedSpacecraftStates.get(i).getAbsPVA();
482 final Vector3D position = currentAPV.getPosition().subtract(deltaP);
483 final Vector3D velocity = currentAPV.getVelocity().subtract(deltaV);
484 final PVCoordinates pv = new PVCoordinates(position, velocity);
485
486 // Update epoch in the AbsolutePVCoordinates
487 AbsoluteDate epoch = currentAPV.getDate();
488 if (!isAutonomous) {
489 if (freeEpochMap[i]) {
490 epoch = epoch.shiftedBy(-dx.getEntry(indexEpoch) * scaleTime);
491 indexEpoch++;
492 }
493 } else {
494 // for autonomous systems we arbitrarily fix the date of the first patch point
495 if (i > 0) {
496 epoch = patchedSpacecraftStates.get(i - 1).getDate().shiftedBy(propagationTime[i - 1]);
497 }
498 }
499 final AbsolutePVCoordinates updatedAPV = new AbsolutePVCoordinates(currentAPV.getFrame(), epoch, pv);
500
501 // Update attitude epoch
502 // Last point does not have an associated propagator. The previous one is then selected.
503 final int iAttitude = i < getPropagatorList().size() ? i : getPropagatorList().size() - 1;
504 final AttitudeProvider attitudeProvider = getPropagatorList().get(iAttitude).getAttitudeProvider();
505 final Attitude attitude = attitudeProvider.getAttitude(updatedAPV, epoch, currentAPV.getFrame());
506
507 // Update the SpacecraftState using previously updated attitude and AbsolutePVCoordinates
508 patchedSpacecraftStates.set(i, new SpacecraftState(updatedAPV, attitude));
509 }
510
511 }
512
513 /** Propagate the patch point states.
514 * @return propagatedSP propagated SpacecraftStates
515 */
516 private List<SpacecraftState> propagatePatchedSpacecraftState() {
517
518 final int nArcs = patchedSpacecraftStates.size() - 1;
519 final ArrayList<SpacecraftState> propagatedSP = new ArrayList<>(nArcs);
520
521 for (int i = 0; i < nArcs; i++) {
522
523 // SpacecraftState initialization
524 final SpacecraftState augmentedInitialState = getAugmentedInitialState(i);
525
526 // Propagator initialization
527 propagatorList.get(i).setInitialState(augmentedInitialState);
528
529 // Propagate trajectory
530 final AbsoluteDate target = augmentedInitialState.getDate().shiftedBy(propagationTime[i]);
531 final SpacecraftState finalState = propagatorList.get(i).propagate(target);
532
533 propagatedSP.add(finalState);
534 }
535
536 return propagatedSP;
537 }
538
539 /** Get the state transition matrix.
540 * @param s current spacecraft state
541 * @return the state transition matrix
542 */
543 private double[][] getStateTransitionMatrix(final SpacecraftState s) {
544 // Additional states
545 final DataDictionary dictionary = s.getAdditionalDataValues();
546 // Initialize state transition matrix
547 final int dim = 6;
548 final double[][] phiM = new double[dim][dim];
549
550 // Loop on entry set
551 for (final DataDictionary.Entry entry : dictionary.getData()) {
552 // Extract entry name
553 final String name = entry.getKey();
554 if (additionalName.equals(name) && entry.getValue() instanceof double[]) {
555 final double[] stm = (double[]) entry.getValue();
556 for (int i = 0; i < 3; i++) {
557 for (int j = 0; j < 3; j++) {
558
559 // partials of final position w.r.t. initial position
560 phiM[i][j] = stm[6 * i + j];
561
562 // partials of final position w.r.t. initial velocity
563 phiM[i][j + 3] = stm[6 * i + j + 3] / scaleTime;
564
565 // partials of final velocity w.r.t. initial position
566 phiM[i + 3][j] = stm[6 * i + j + 18] * scaleTime;
567
568 // partials of final velocity w.r.t. initial velocity
569 phiM[i + 3][j + 3] = stm[6 * i + j + 21];
570 }
571 }
572 }
573 }
574
575 // Return state transition matrix
576 return phiM;
577 }
578
579 /** Compute a part of the Jacobian matrix with derivatives from epoch.
580 * The CR3BP is a time invariant problem. The derivatives w.r.t. epoch are zero.
581 * @param propagatedSP propagatedSP
582 * @return jacobianMatrix Jacobian sub-matrix
583 */
584 protected double[][] computeEpochJacobianMatrix(final List<SpacecraftState> propagatedSP) {
585
586 final int nRows = patchedSpacecraftStates.size() - 1;
587 final double[][] M = new double[nRows][nDuration + nEpoch];
588
589 // The D and E sub-matrices have the following form:
590 //
591 // D = -I
592 //
593 // [-1 +1 0 ]
594 // [ -1 +1 0 ]
595 // F = [ .. .. ]
596 // [ .. .. 0 ]
597 // [ -1 +1 ]
598
599 int index = nDuration;
600 for (int i = 0; i < nRows; i++) {
601
602 // components of D matrix
603 M[i][i] = -1.0;
604
605 // components of E matrix
606 if (freeEpochMap[i]) {
607 M[i][index] = -1.0;
608 index++;
609 }
610 if (freeEpochMap[i + 1]) {
611 M[i][index] = 1.0;
612 }
613 }
614
615 return M;
616 }
617
618 /** Update the array of additional constraints.
619 * @param startIndex start index
620 * @param fxAdditional array of additional constraints
621 */
622 protected void updateAdditionalConstraints(final int startIndex, final double[] fxAdditional) {
623 int iter = startIndex;
624 final double scaleVel = scaleLength / scaleTime;
625 for (final Map.Entry<Integer, Double> entry : getConstraintsMap().entrySet()) {
626 // Extract entry values
627 final int key = entry.getKey();
628 final double value = entry.getValue();
629 final int np = key / 6;
630 final int nc = key % 6;
631 final AbsolutePVCoordinates absPv = getPatchedSpacecraftState().get(np).getAbsPVA();
632 if (nc < 3) {
633 fxAdditional[iter] = (absPv.getPosition().toArray()[nc] - value) / scaleLength;
634 } else {
635 fxAdditional[iter] = (absPv.getVelocity().toArray()[nc - 3] - value) / scaleVel;
636 }
637 iter++;
638 }
639 }
640
641 /** Compute the additional constraints.
642 * @param propagatedSP propagated SpacecraftState
643 * @return fxAdditional additional constraints
644 */
645 protected abstract double[] computeAdditionalConstraints(List<SpacecraftState> propagatedSP);
646
647 /** Compute a part of the Jacobian matrix from additional constraints.
648 * @param propagatedSP propagatedSP
649 * @return jacobianMatrix Jacobian sub-matrix
650 */
651 protected abstract double[][] computeAdditionalJacobianMatrix(List<SpacecraftState> propagatedSP);
652
653 /** Compute the additional state from the additionalEquations.
654 * @param i index of the state
655 * @return augmentedSP SpacecraftState with the additional state within.
656 * @since 11.1
657 */
658 protected abstract SpacecraftState getAugmentedInitialState(int i);
659
660 /** Get the number of free state components.
661 * @return number of free components
662 */
663 protected int getNumberOfFreeComponents() {
664 return nComps;
665 }
666
667 /** Get the total number of constraints.
668 * @return the total number of constraints
669 */
670 protected int getNumberOfConstraints() {
671 final int nArcs = patchedSpacecraftStates.size() - 1;
672 return (isAutonomous ? 6 * nArcs : 7 * nArcs) + mapConstraints.size();
673 }
674
675 /** Get the map of free state components.
676 * @return map of free state components
677 */
678 protected boolean[] getFreeCompsMap() {
679 return freeCompsMap;
680 }
681
682 /** Get the map of patch points components which are constrained.
683 * @return a map of patch points components which are constrained
684 */
685 protected Map<Integer, Double> getConstraintsMap() {
686 return mapConstraints;
687 }
688
689 /** Get the list of patched spacecraft states.
690 * @return a list of patched spacecraft states
691 */
692 protected List<SpacecraftState> getPatchedSpacecraftState() {
693 return patchedSpacecraftStates;
694 }
695
696 /** Get the list of propagators.
697 * @return a list of propagators
698 */
699 private List<NumericalPropagator> getPropagatorList() {
700 return propagatorList;
701 }
702
703 /** Get the number of free variables.
704 * @return the number of free variables
705 */
706 private int getNumberOfFreeVariables() {
707 return nComps + nDuration + nEpoch;
708 }
709
710 }