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