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.estimation.sequential;
18  
19  import java.util.ArrayList;
20  import java.util.Arrays;
21  import java.util.Comparator;
22  import java.util.HashMap;
23  import java.util.List;
24  import java.util.Map;
25  
26  import org.hipparchus.filtering.kalman.ProcessEstimate;
27  import org.hipparchus.filtering.kalman.extended.NonLinearEvolution;
28  import org.hipparchus.filtering.kalman.extended.NonLinearProcess;
29  import org.hipparchus.linear.Array2DRowRealMatrix;
30  import org.hipparchus.linear.ArrayRealVector;
31  import org.hipparchus.linear.MatrixUtils;
32  import org.hipparchus.linear.RealMatrix;
33  import org.hipparchus.linear.RealVector;
34  import org.hipparchus.util.FastMath;
35  import org.orekit.errors.OrekitException;
36  import org.orekit.errors.OrekitMessages;
37  import org.orekit.estimation.measurements.EstimatedMeasurement;
38  import org.orekit.estimation.measurements.EstimationModifier;
39  import org.orekit.estimation.measurements.ObservableSatellite;
40  import org.orekit.estimation.measurements.ObservedMeasurement;
41  import org.orekit.estimation.measurements.modifiers.DynamicOutlierFilter;
42  import org.orekit.orbits.Orbit;
43  import org.orekit.propagation.MatricesHarvester;
44  import org.orekit.propagation.PropagationType;
45  import org.orekit.propagation.Propagator;
46  import org.orekit.propagation.SpacecraftState;
47  import org.orekit.propagation.conversion.OrbitDeterminationPropagatorBuilder;
48  import org.orekit.propagation.integration.AbstractJacobiansMapper;
49  import org.orekit.time.AbsoluteDate;
50  import org.orekit.utils.ParameterDriver;
51  import org.orekit.utils.ParameterDriversList;
52  import org.orekit.utils.ParameterDriversList.DelegatingDriver;
53  
54  /** Abstract class defining the process model dynamics to use with a {@link KalmanEstimator}.
55   * @author Romain Gerbaud
56   * @author Maxime Journot
57   * @author Bryan Cazabonne
58   * @author Thomas Paulet
59   * @since 11.0
60   */
61  public abstract class AbstractKalmanModel implements KalmanEstimation, NonLinearProcess<MeasurementDecorator> {
62  
63      /** Builders for propagators. */
64      private final List<OrbitDeterminationPropagatorBuilder> builders;
65  
66      /** Estimated orbital parameters. */
67      private final ParameterDriversList allEstimatedOrbitalParameters;
68  
69      /** Estimated propagation drivers. */
70      private final ParameterDriversList allEstimatedPropagationParameters;
71  
72      /** Per-builder estimated orbita parameters drivers.
73       * @since 11.1
74       */
75      private final ParameterDriversList[] estimatedOrbitalParameters;
76  
77      /** Per-builder estimated propagation drivers. */
78      private final ParameterDriversList[] estimatedPropagationParameters;
79  
80      /** Estimated measurements parameters. */
81      private final ParameterDriversList estimatedMeasurementsParameters;
82  
83      /** Start columns for each estimated orbit. */
84      private final int[] orbitsStartColumns;
85  
86      /** End columns for each estimated orbit. */
87      private final int[] orbitsEndColumns;
88  
89      /** Map for propagation parameters columns. */
90      private final Map<String, Integer> propagationParameterColumns;
91  
92      /** Map for measurements parameters columns. */
93      private final Map<String, Integer> measurementParameterColumns;
94  
95      /** Providers for covariance matrices. */
96      private final List<CovarianceMatrixProvider> covarianceMatricesProviders;
97  
98      /** Process noise matrix provider for measurement parameters. */
99      private final CovarianceMatrixProvider measurementProcessNoiseMatrix;
100 
101     /** Indirection arrays to extract the noise components for estimated parameters. */
102     private final int[][] covarianceIndirection;
103 
104     /** Scaling factors. */
105     private final double[] scale;
106 
107     /** Harvesters for extracting Jacobians from integrated states. */
108     private MatricesHarvester[] harvesters;
109 
110     /** Propagators for the reference trajectories, up to current date. */
111     private Propagator[] referenceTrajectories;
112 
113     /** Current corrected estimate. */
114     private ProcessEstimate correctedEstimate;
115 
116     /** Current number of measurement. */
117     private int currentMeasurementNumber;
118 
119     /** Reference date. */
120     private AbsoluteDate referenceDate;
121 
122     /** Current date. */
123     private AbsoluteDate currentDate;
124 
125     /** Predicted spacecraft states. */
126     private SpacecraftState[] predictedSpacecraftStates;
127 
128     /** Corrected spacecraft states. */
129     private SpacecraftState[] correctedSpacecraftStates;
130 
131     /** Predicted measurement. */
132     private EstimatedMeasurement<?> predictedMeasurement;
133 
134     /** Corrected measurement. */
135     private EstimatedMeasurement<?> correctedMeasurement;
136 
137     /** Type of the orbit used for the propagation.*/
138     private PropagationType propagationType;
139 
140     /** Type of the elements used to define the orbital state.*/
141     private PropagationType stateType;
142 
143     /** Kalman process model constructor (package private).
144      * This constructor is used whenever state type and propagation type do not matter.
145      * It is used for {@link KalmanModel} and {@link TLEKalmanModel}.
146      * @param propagatorBuilders propagators builders used to evaluate the orbits.
147      * @param covarianceMatricesProviders providers for covariance matrices
148      * @param estimatedMeasurementParameters measurement parameters to estimate
149      * @param measurementProcessNoiseMatrix provider for measurement process noise matrix
150      * @param harvesters harvesters for extracting Jacobians from integrated states
151      */
152     protected AbstractKalmanModel(final List<OrbitDeterminationPropagatorBuilder> propagatorBuilders,
153                                   final List<CovarianceMatrixProvider> covarianceMatricesProviders,
154                                   final ParameterDriversList estimatedMeasurementParameters,
155                                   final CovarianceMatrixProvider measurementProcessNoiseMatrix,
156                                   final MatricesHarvester[] harvesters) {
157         this(propagatorBuilders, covarianceMatricesProviders, estimatedMeasurementParameters,
158              measurementProcessNoiseMatrix, harvesters, PropagationType.MEAN, PropagationType.MEAN);
159     }
160 
161     /** Kalman process model constructor (package private).
162      * This constructor is used whenever propagation type and/or state type are to be specified.
163      * It is used for {@link DSSTKalmanModel}.
164      * @param propagatorBuilders propagators builders used to evaluate the orbits.
165      * @param covarianceMatricesProviders providers for covariance matrices
166      * @param estimatedMeasurementParameters measurement parameters to estimate
167      * @param measurementProcessNoiseMatrix provider for measurement process noise matrix
168      * @param harvesters harvesters for extracting Jacobians from integrated states
169      * @param propagationType type of the orbit used for the propagation (mean or osculating), applicable only for DSST
170      * @param stateType type of the elements used to define the orbital state (mean or osculating), applicable only for DSST
171      */
172     protected AbstractKalmanModel(final List<OrbitDeterminationPropagatorBuilder> propagatorBuilders,
173                                   final List<CovarianceMatrixProvider> covarianceMatricesProviders,
174                                   final ParameterDriversList estimatedMeasurementParameters,
175                                   final CovarianceMatrixProvider measurementProcessNoiseMatrix,
176                                   final MatricesHarvester[] harvesters,
177                                   final PropagationType propagationType,
178                                   final PropagationType stateType) {
179 
180         this.builders                        = propagatorBuilders;
181         this.estimatedMeasurementsParameters = estimatedMeasurementParameters;
182         this.measurementParameterColumns     = new HashMap<>(estimatedMeasurementsParameters.getDrivers().size());
183         this.currentMeasurementNumber        = 0;
184         this.referenceDate                   = propagatorBuilders.get(0).getInitialOrbitDate();
185         this.currentDate                     = referenceDate;
186         this.propagationType                 = propagationType;
187         this.stateType                       = stateType;
188 
189         final Map<String, Integer> orbitalParameterColumns = new HashMap<>(6 * builders.size());
190         orbitsStartColumns      = new int[builders.size()];
191         orbitsEndColumns        = new int[builders.size()];
192         int columns = 0;
193         allEstimatedOrbitalParameters = new ParameterDriversList();
194         estimatedOrbitalParameters    = new ParameterDriversList[builders.size()];
195         for (int k = 0; k < builders.size(); ++k) {
196             estimatedOrbitalParameters[k] = new ParameterDriversList();
197             orbitsStartColumns[k] = columns;
198             final String suffix = propagatorBuilders.size() > 1 ? "[" + k + "]" : null;
199             for (final ParameterDriver driver : builders.get(k).getOrbitalParametersDrivers().getDrivers()) {
200                 if (driver.getReferenceDate() == null) {
201                     driver.setReferenceDate(currentDate);
202                 }
203                 if (suffix != null && !driver.getName().endsWith(suffix)) {
204                     // we add suffix only conditionally because the method may already have been called
205                     // and suffixes may have already been appended
206                     driver.setName(driver.getName() + suffix);
207                 }
208                 if (driver.isSelected()) {
209                     allEstimatedOrbitalParameters.add(driver);
210                     estimatedOrbitalParameters[k].add(driver);
211                     orbitalParameterColumns.put(driver.getName(), columns++);
212                 }
213             }
214             orbitsEndColumns[k] = columns;
215         }
216 
217         // Gather all the propagation drivers names in a list
218         allEstimatedPropagationParameters = new ParameterDriversList();
219         estimatedPropagationParameters    = new ParameterDriversList[builders.size()];
220         final List<String> estimatedPropagationParametersNames = new ArrayList<>();
221         for (int k = 0; k < builders.size(); ++k) {
222             estimatedPropagationParameters[k] = new ParameterDriversList();
223             for (final ParameterDriver driver : builders.get(k).getPropagationParametersDrivers().getDrivers()) {
224                 if (driver.getReferenceDate() == null) {
225                     driver.setReferenceDate(currentDate);
226                 }
227                 if (driver.isSelected()) {
228                     allEstimatedPropagationParameters.add(driver);
229                     estimatedPropagationParameters[k].add(driver);
230                     final String driverName = driver.getName();
231                     // Add the driver name if it has not been added yet
232                     if (!estimatedPropagationParametersNames.contains(driverName)) {
233                         estimatedPropagationParametersNames.add(driverName);
234                     }
235                 }
236             }
237         }
238         estimatedPropagationParametersNames.sort(Comparator.naturalOrder());
239 
240         // Populate the map of propagation drivers' columns and update the total number of columns
241         propagationParameterColumns = new HashMap<>(estimatedPropagationParametersNames.size());
242         for (final String driverName : estimatedPropagationParametersNames) {
243             propagationParameterColumns.put(driverName, columns);
244             ++columns;
245         }
246 
247         // Populate the map of measurement drivers' columns and update the total number of columns
248         for (final ParameterDriver parameter : estimatedMeasurementsParameters.getDrivers()) {
249             if (parameter.getReferenceDate() == null) {
250                 parameter.setReferenceDate(currentDate);
251             }
252             measurementParameterColumns.put(parameter.getName(), columns);
253             ++columns;
254         }
255 
256         // Store providers for process noise matrices
257         this.covarianceMatricesProviders = covarianceMatricesProviders;
258         this.measurementProcessNoiseMatrix = measurementProcessNoiseMatrix;
259         this.covarianceIndirection       = new int[covarianceMatricesProviders.size()][columns];
260         for (int k = 0; k < covarianceIndirection.length; ++k) {
261             final ParameterDriversList orbitDrivers      = builders.get(k).getOrbitalParametersDrivers();
262             final ParameterDriversList parametersDrivers = builders.get(k).getPropagationParametersDrivers();
263             Arrays.fill(covarianceIndirection[k], -1);
264             int i = 0;
265             for (final ParameterDriver driver : orbitDrivers.getDrivers()) {
266                 final Integer c = orbitalParameterColumns.get(driver.getName());
267                 covarianceIndirection[k][i++] = (c == null) ? -1 : c.intValue();
268             }
269             for (final ParameterDriver driver : parametersDrivers.getDrivers()) {
270                 final Integer c = propagationParameterColumns.get(driver.getName());
271                 if (c != null) {
272                     covarianceIndirection[k][i++] = c.intValue();
273                 }
274             }
275             for (final ParameterDriver driver : estimatedMeasurementParameters.getDrivers()) {
276                 final Integer c = measurementParameterColumns.get(driver.getName());
277                 if (c != null) {
278                     covarianceIndirection[k][i++] = c.intValue();
279                 }
280             }
281         }
282 
283         // Compute the scale factors
284         this.scale = new double[columns];
285         int index = 0;
286         for (final ParameterDriver driver : allEstimatedOrbitalParameters.getDrivers()) {
287             scale[index++] = driver.getScale();
288         }
289         for (final ParameterDriver driver : allEstimatedPropagationParameters.getDrivers()) {
290             scale[index++] = driver.getScale();
291         }
292         for (final ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
293             scale[index++] = driver.getScale();
294         }
295 
296         // Build the reference propagators and add their partial derivatives equations implementation
297         this.harvesters = harvesters.clone();
298         updateReferenceTrajectories(getEstimatedPropagators(), propagationType, stateType);
299         this.predictedSpacecraftStates = new SpacecraftState[referenceTrajectories.length];
300         for (int i = 0; i < predictedSpacecraftStates.length; ++i) {
301             predictedSpacecraftStates[i] = referenceTrajectories[i].getInitialState();
302         };
303         this.correctedSpacecraftStates = predictedSpacecraftStates.clone();
304 
305         // Initialize the estimated normalized state and fill its values
306         final RealVector correctedState      = MatrixUtils.createRealVector(columns);
307 
308         int p = 0;
309         for (final ParameterDriver driver : allEstimatedOrbitalParameters.getDrivers()) {
310             correctedState.setEntry(p++, driver.getNormalizedValue());
311         }
312         for (final ParameterDriver driver : allEstimatedPropagationParameters.getDrivers()) {
313             correctedState.setEntry(p++, driver.getNormalizedValue());
314         }
315         for (final ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
316             correctedState.setEntry(p++, driver.getNormalizedValue());
317         }
318 
319         // Set up initial covariance
320         final RealMatrix physicalProcessNoise = MatrixUtils.createRealMatrix(columns, columns);
321         for (int k = 0; k < covarianceMatricesProviders.size(); ++k) {
322 
323             // Number of estimated measurement parameters
324             final int nbMeas = estimatedMeasurementParameters.getNbParams();
325 
326             // Number of estimated dynamic parameters (orbital + propagation)
327             final int nbDyn  = orbitsEndColumns[k] - orbitsStartColumns[k] +
328                                estimatedPropagationParameters[k].getNbParams();
329 
330             // Covariance matrix
331             final RealMatrix noiseK = MatrixUtils.createRealMatrix(nbDyn + nbMeas, nbDyn + nbMeas);
332             final RealMatrix noiseP = covarianceMatricesProviders.get(k).
333                                       getInitialCovarianceMatrix(correctedSpacecraftStates[k]);
334             noiseK.setSubMatrix(noiseP.getData(), 0, 0);
335             if (measurementProcessNoiseMatrix != null) {
336                 final RealMatrix noiseM = measurementProcessNoiseMatrix.
337                                           getInitialCovarianceMatrix(correctedSpacecraftStates[k]);
338                 noiseK.setSubMatrix(noiseM.getData(), nbDyn, nbDyn);
339             }
340 
341             checkDimension(noiseK.getRowDimension(),
342                            builders.get(k).getOrbitalParametersDrivers(),
343                            builders.get(k).getPropagationParametersDrivers(),
344                            estimatedMeasurementsParameters);
345 
346             final int[] indK = covarianceIndirection[k];
347             for (int i = 0; i < indK.length; ++i) {
348                 if (indK[i] >= 0) {
349                     for (int j = 0; j < indK.length; ++j) {
350                         if (indK[j] >= 0) {
351                             physicalProcessNoise.setEntry(indK[i], indK[j], noiseK.getEntry(i, j));
352                         }
353                     }
354                 }
355             }
356 
357         }
358         final RealMatrix correctedCovariance = normalizeCovarianceMatrix(physicalProcessNoise);
359 
360         correctedEstimate = new ProcessEstimate(0.0, correctedState, correctedCovariance);
361 
362     }
363 
364     /** Update the reference trajectories using the propagators as input.
365      * @param propagators The new propagators to use
366      * @param pType propagationType type of the orbit used for the propagation (mean or osculating)
367      * @param sType type of the elements used to define the orbital state (mean or osculating)
368      */
369     protected abstract void updateReferenceTrajectories(Propagator[] propagators,
370                                                         PropagationType pType,
371                                                         PropagationType sType);
372 
373     /** Not used anymore.
374      * @param mapper Jacobian mapper to calculate short period perturbations
375      * @param state state used to calculate short period perturbations
376      * @deprecated as of 11.1, not used anymore
377      */
378     protected void analyticalDerivativeComputations(final AbstractJacobiansMapper mapper, final SpacecraftState state) {
379         // nothing by default
380     }
381 
382     /** Check dimension.
383      * @param dimension dimension to check
384      * @param orbitalParameters orbital parameters
385      * @param propagationParameters propagation parameters
386      * @param measurementParameters measurements parameters
387      */
388     private void checkDimension(final int dimension,
389                                 final ParameterDriversList orbitalParameters,
390                                 final ParameterDriversList propagationParameters,
391                                 final ParameterDriversList measurementParameters) {
392 
393         // count parameters, taking care of counting all orbital parameters
394         // regardless of them being estimated or not
395         int requiredDimension = orbitalParameters.getNbParams();
396         for (final ParameterDriver driver : propagationParameters.getDrivers()) {
397             if (driver.isSelected()) {
398                 ++requiredDimension;
399             }
400         }
401         for (final ParameterDriver driver : measurementParameters.getDrivers()) {
402             if (driver.isSelected()) {
403                 ++requiredDimension;
404             }
405         }
406 
407         if (dimension != requiredDimension) {
408             // there is a problem, set up an explicit error message
409             final StringBuilder builder = new StringBuilder();
410             for (final ParameterDriver driver : orbitalParameters.getDrivers()) {
411                 if (builder.length() > 0) {
412                     builder.append(", ");
413                 }
414                 builder.append(driver.getName());
415             }
416             for (final ParameterDriver driver : propagationParameters.getDrivers()) {
417                 if (driver.isSelected()) {
418                     builder.append(driver.getName());
419                 }
420             }
421             for (final ParameterDriver driver : measurementParameters.getDrivers()) {
422                 if (driver.isSelected()) {
423                     builder.append(driver.getName());
424                 }
425             }
426             throw new OrekitException(OrekitMessages.DIMENSION_INCONSISTENT_WITH_PARAMETERS,
427                                       dimension, builder.toString());
428         }
429 
430     }
431 
432     /** {@inheritDoc} */
433     @Override
434     public RealMatrix getPhysicalStateTransitionMatrix() {
435         //  Un-normalize the state transition matrix (ฯ†) from Hipparchus and return it.
436         // ฯ† is an mxm matrix where m = nbOrb + nbPropag + nbMeas
437         // For each element [i,j] of normalized ฯ† (ฯ†n), the corresponding physical value is:
438         // ฯ†[i,j] = ฯ†n[i,j] * scale[i] / scale[j]
439 
440         // Normalized matrix
441         final RealMatrix normalizedSTM = correctedEstimate.getStateTransitionMatrix();
442 
443         if (normalizedSTM == null) {
444             return null;
445         } else {
446             // Initialize physical matrix
447             final int nbParams = normalizedSTM.getRowDimension();
448             final RealMatrix physicalSTM = MatrixUtils.createRealMatrix(nbParams, nbParams);
449 
450             // Un-normalize the matrix
451             for (int i = 0; i < nbParams; ++i) {
452                 for (int j = 0; j < nbParams; ++j) {
453                     physicalSTM.setEntry(i, j,
454                                          normalizedSTM.getEntry(i, j) * scale[i] / scale[j]);
455                 }
456             }
457             return physicalSTM;
458         }
459     }
460 
461     /** {@inheritDoc} */
462     @Override
463     public RealMatrix getPhysicalMeasurementJacobian() {
464         // Un-normalize the measurement matrix (H) from Hipparchus and return it.
465         // H is an nxm matrix where:
466         //  - m = nbOrb + nbPropag + nbMeas is the number of estimated parameters
467         //  - n is the size of the measurement being processed by the filter
468         // For each element [i,j] of normalized H (Hn) the corresponding physical value is:
469         // H[i,j] = Hn[i,j] * ฯƒ[i] / scale[j]
470 
471         // Normalized matrix
472         final RealMatrix normalizedH = correctedEstimate.getMeasurementJacobian();
473 
474         if (normalizedH == null) {
475             return null;
476         } else {
477             // Get current measurement sigmas
478             final double[] sigmas = correctedMeasurement.getObservedMeasurement().getTheoreticalStandardDeviation();
479 
480             // Initialize physical matrix
481             final int nbLine = normalizedH.getRowDimension();
482             final int nbCol  = normalizedH.getColumnDimension();
483             final RealMatrix physicalH = MatrixUtils.createRealMatrix(nbLine, nbCol);
484 
485             // Un-normalize the matrix
486             for (int i = 0; i < nbLine; ++i) {
487                 for (int j = 0; j < nbCol; ++j) {
488                     physicalH.setEntry(i, j, normalizedH.getEntry(i, j) * sigmas[i] / scale[j]);
489                 }
490             }
491             return physicalH;
492         }
493     }
494 
495     /** {@inheritDoc} */
496     @Override
497     public RealMatrix getPhysicalInnovationCovarianceMatrix() {
498         // Un-normalize the innovation covariance matrix (S) from Hipparchus and return it.
499         // S is an nxn matrix where n is the size of the measurement being processed by the filter
500         // For each element [i,j] of normalized S (Sn) the corresponding physical value is:
501         // S[i,j] = Sn[i,j] * ฯƒ[i] * ฯƒ[j]
502 
503         // Normalized matrix
504         final RealMatrix normalizedS = correctedEstimate.getInnovationCovariance();
505 
506         if (normalizedS == null) {
507             return null;
508         } else {
509             // Get current measurement sigmas
510             final double[] sigmas = correctedMeasurement.getObservedMeasurement().getTheoreticalStandardDeviation();
511 
512             // Initialize physical matrix
513             final int nbMeas = sigmas.length;
514             final RealMatrix physicalS = MatrixUtils.createRealMatrix(nbMeas, nbMeas);
515 
516             // Un-normalize the matrix
517             for (int i = 0; i < nbMeas; ++i) {
518                 for (int j = 0; j < nbMeas; ++j) {
519                     physicalS.setEntry(i, j, normalizedS.getEntry(i, j) * sigmas[i] *   sigmas[j]);
520                 }
521             }
522             return physicalS;
523         }
524     }
525 
526     /** {@inheritDoc} */
527     @Override
528     public RealMatrix getPhysicalKalmanGain() {
529         // Un-normalize the Kalman gain (K) from Hipparchus and return it.
530         // K is an mxn matrix where:
531         //  - m = nbOrb + nbPropag + nbMeas is the number of estimated parameters
532         //  - n is the size of the measurement being processed by the filter
533         // For each element [i,j] of normalized K (Kn) the corresponding physical value is:
534         // K[i,j] = Kn[i,j] * scale[i] / ฯƒ[j]
535 
536         // Normalized matrix
537         final RealMatrix normalizedK = correctedEstimate.getKalmanGain();
538 
539         if (normalizedK == null) {
540             return null;
541         } else {
542             // Get current measurement sigmas
543             final double[] sigmas = correctedMeasurement.getObservedMeasurement().getTheoreticalStandardDeviation();
544 
545             // Initialize physical matrix
546             final int nbLine = normalizedK.getRowDimension();
547             final int nbCol  = normalizedK.getColumnDimension();
548             final RealMatrix physicalK = MatrixUtils.createRealMatrix(nbLine, nbCol);
549 
550             // Un-normalize the matrix
551             for (int i = 0; i < nbLine; ++i) {
552                 for (int j = 0; j < nbCol; ++j) {
553                     physicalK.setEntry(i, j, normalizedK.getEntry(i, j) * scale[i] / sigmas[j]);
554                 }
555             }
556             return physicalK;
557         }
558     }
559 
560     /** {@inheritDoc} */
561     @Override
562     public SpacecraftState[] getPredictedSpacecraftStates() {
563         return predictedSpacecraftStates.clone();
564     }
565 
566     /** {@inheritDoc} */
567     @Override
568     public SpacecraftState[] getCorrectedSpacecraftStates() {
569         return correctedSpacecraftStates.clone();
570     }
571 
572     /** {@inheritDoc} */
573     @Override
574     public int getCurrentMeasurementNumber() {
575         return currentMeasurementNumber;
576     }
577 
578     /** {@inheritDoc} */
579     @Override
580     public AbsoluteDate getCurrentDate() {
581         return currentDate;
582     }
583 
584     /** {@inheritDoc} */
585     @Override
586     public EstimatedMeasurement<?> getPredictedMeasurement() {
587         return predictedMeasurement;
588     }
589 
590     /** {@inheritDoc} */
591     @Override
592     public EstimatedMeasurement<?> getCorrectedMeasurement() {
593         return correctedMeasurement;
594     }
595 
596     /** {@inheritDoc} */
597     @Override
598     public RealVector getPhysicalEstimatedState() {
599         // Method {@link ParameterDriver#getValue()} is used to get
600         // the physical values of the state.
601         // The scales'array is used to get the size of the state vector
602         final RealVector physicalEstimatedState = new ArrayRealVector(scale.length);
603         int i = 0;
604         for (final DelegatingDriver driver : getEstimatedOrbitalParameters().getDrivers()) {
605             physicalEstimatedState.setEntry(i++, driver.getValue());
606         }
607         for (final DelegatingDriver driver : getEstimatedPropagationParameters().getDrivers()) {
608             physicalEstimatedState.setEntry(i++, driver.getValue());
609         }
610         for (final DelegatingDriver driver : getEstimatedMeasurementsParameters().getDrivers()) {
611             physicalEstimatedState.setEntry(i++, driver.getValue());
612         }
613 
614         return physicalEstimatedState;
615     }
616 
617     /** {@inheritDoc} */
618     @Override
619     public RealMatrix getPhysicalEstimatedCovarianceMatrix() {
620         // Un-normalize the estimated covariance matrix (P) from Hipparchus and return it.
621         // The covariance P is an mxm matrix where m = nbOrb + nbPropag + nbMeas
622         // For each element [i,j] of P the corresponding normalized value is:
623         // Pn[i,j] = P[i,j] / (scale[i]*scale[j])
624         // Consequently: P[i,j] = Pn[i,j] * scale[i] * scale[j]
625 
626         // Normalized covariance matrix
627         final RealMatrix normalizedP = correctedEstimate.getCovariance();
628 
629         // Initialize physical covariance matrix
630         final int nbParams = normalizedP.getRowDimension();
631         final RealMatrix physicalP = MatrixUtils.createRealMatrix(nbParams, nbParams);
632 
633         // Un-normalize the covairance matrix
634         for (int i = 0; i < nbParams; ++i) {
635             for (int j = 0; j < nbParams; ++j) {
636                 physicalP.setEntry(i, j, normalizedP.getEntry(i, j) * scale[i] * scale[j]);
637             }
638         }
639         return physicalP;
640     }
641 
642     /** {@inheritDoc} */
643     @Override
644     public ParameterDriversList getEstimatedOrbitalParameters() {
645         return allEstimatedOrbitalParameters;
646     }
647 
648     /** {@inheritDoc} */
649     @Override
650     public ParameterDriversList getEstimatedPropagationParameters() {
651         return allEstimatedPropagationParameters;
652     }
653 
654     /** {@inheritDoc} */
655     @Override
656     public ParameterDriversList getEstimatedMeasurementsParameters() {
657         return estimatedMeasurementsParameters;
658     }
659 
660     /** Get the current corrected estimate.
661      * @return current corrected estimate
662      */
663     public ProcessEstimate getEstimate() {
664         return correctedEstimate;
665     }
666 
667     /** Get the normalized error state transition matrix (STM) from previous point to current point.
668      * The STM contains the partial derivatives of current state with respect to previous state.
669      * The  STM is an mxm matrix where m is the size of the state vector.
670      * m = nbOrb + nbPropag + nbMeas
671      * @return the normalized error state transition matrix
672      */
673     private RealMatrix getErrorStateTransitionMatrix() {
674 
675         /* The state transition matrix is obtained as follows, with:
676          *  - Y  : Current state vector
677          *  - Y0 : Initial state vector
678          *  - Pp : Current propagation parameter
679          *  - Pp0: Initial propagation parameter
680          *  - Mp : Current measurement parameter
681          *  - Mp0: Initial measurement parameter
682          *
683          *       |        |         |         |   |        |        |   .    |
684          *       | dY/dY0 | dY/dPp  | dY/dMp  |   | dY/dY0 | dY/dPp | ..0..  |
685          *       |        |         |         |   |        |        |   .    |
686          *       |--------|---------|---------|   |--------|--------|--------|
687          *       |        |         |         |   |   .    | 1 0 0..|   .    |
688          * STM = | dP/dY0 | dP/dPp0 | dP/dMp  | = | ..0..  | 0 1 0..| ..0..  |
689          *       |        |         |         |   |   .    | 0 0 1..|   .    |
690          *       |--------|---------|---------|   |--------|--------|--------|
691          *       |        |         |         |   |   .    |   .    | 1 0 0..|
692          *       | dM/dY0 | dM/dPp0 | dM/dMp0 |   | ..0..  | ..0..  | 0 1 0..|
693          *       |        |         |         |   |   .    |   .    | 0 0 1..|
694          */
695 
696         // Initialize to the proper size identity matrix
697         final RealMatrix stm = MatrixUtils.createRealIdentityMatrix(correctedEstimate.getState().getDimension());
698 
699         // loop over all orbits
700         for (int k = 0; k < predictedSpacecraftStates.length; ++k) {
701 
702             // Indexes
703             final int[] indK = covarianceIndirection[k];
704 
705             // Reset reference (for example compute short periodic terms in DSST)
706             harvesters[k].setReferenceState(predictedSpacecraftStates[k]);
707 
708             // Derivatives of the state vector with respect to initial state vector
709             final int nbOrbParams = estimatedOrbitalParameters[k].getNbParams();
710             if (nbOrbParams > 0) {
711                 final RealMatrix dYdY0 = harvesters[k].getStateTransitionMatrix(predictedSpacecraftStates[k]);
712 
713                 // Fill upper left corner (dY/dY0)
714                 for (int i = 0; i < dYdY0.getRowDimension(); ++i) {
715                     for (int j = 0; j < nbOrbParams; ++j) {
716                         stm.setEntry(indK[i], indK[j], dYdY0.getEntry(i, j));
717                     }
718                 }
719             }
720 
721             // Derivatives of the state vector with respect to propagation parameters
722             final int nbParams = estimatedPropagationParameters[k].getNbParams();
723             if (nbParams > 0) {
724                 final RealMatrix dYdPp = harvesters[k].getParametersJacobian(predictedSpacecraftStates[k]);
725 
726                 // Fill 1st row, 2nd column (dY/dPp)
727                 for (int i = 0; i < dYdPp.getRowDimension(); ++i) {
728                     for (int j = 0; j < nbParams; ++j) {
729                         stm.setEntry(indK[i], indK[j + 6], dYdPp.getEntry(i, j));
730                     }
731                 }
732 
733             }
734 
735         }
736 
737         // Normalization of the STM
738         // normalized(STM)ij = STMij*Sj/Si
739         for (int i = 0; i < scale.length; i++) {
740             for (int j = 0; j < scale.length; j++ ) {
741                 stm.setEntry(i, j, stm.getEntry(i, j) * scale[j] / scale[i]);
742             }
743         }
744 
745         // Return the error state transition matrix
746         return stm;
747 
748     }
749 
750     /** Get the normalized measurement matrix H.
751      * H contains the partial derivatives of the measurement with respect to the state.
752      * H is an nxm matrix where n is the size of the measurement vector and m the size of the state vector.
753      * @return the normalized measurement matrix H
754      */
755     private RealMatrix getMeasurementMatrix() {
756 
757         // Observed measurement characteristics
758         final SpacecraftState[]      evaluationStates    = predictedMeasurement.getStates();
759         final ObservedMeasurement<?> observedMeasurement = predictedMeasurement.getObservedMeasurement();
760         final double[] sigma  = observedMeasurement.getTheoreticalStandardDeviation();
761 
762         // Initialize measurement matrix H: nxm
763         // n: Number of measurements in current measurement
764         // m: State vector size
765         final RealMatrix measurementMatrix = MatrixUtils.
766                         createRealMatrix(observedMeasurement.getDimension(),
767                                          correctedEstimate.getState().getDimension());
768 
769         // loop over all orbits involved in the measurement
770         for (int k = 0; k < evaluationStates.length; ++k) {
771             final int p = observedMeasurement.getSatellites().get(k).getPropagatorIndex();
772 
773             // Predicted orbit
774             final Orbit predictedOrbit = evaluationStates[k].getOrbit();
775 
776             // Measurement matrix's columns related to orbital parameters
777             // ----------------------------------------------------------
778 
779             // Partial derivatives of the current Cartesian coordinates with respect to current orbital state
780             final double[][] aCY = new double[6][6];
781             predictedOrbit.getJacobianWrtParameters(builders.get(p).getPositionAngle(), aCY);   //dC/dY
782             final RealMatrix dCdY = new Array2DRowRealMatrix(aCY, false);
783 
784             // Jacobian of the measurement with respect to current Cartesian coordinates
785             final RealMatrix dMdC = new Array2DRowRealMatrix(predictedMeasurement.getStateDerivatives(k), false);
786 
787             // Jacobian of the measurement with respect to current orbital state
788             final RealMatrix dMdY = dMdC.multiply(dCdY);
789 
790             // Fill the normalized measurement matrix's columns related to estimated orbital parameters
791             for (int i = 0; i < dMdY.getRowDimension(); ++i) {
792                 int jOrb = orbitsStartColumns[p];
793                 for (int j = 0; j < dMdY.getColumnDimension(); ++j) {
794                     final ParameterDriver driver = builders.get(p).getOrbitalParametersDrivers().getDrivers().get(j);
795                     if (driver.isSelected()) {
796                         measurementMatrix.setEntry(i, jOrb++,
797                                                    dMdY.getEntry(i, j) / sigma[i] * driver.getScale());
798                     }
799                 }
800             }
801 
802             // Normalized measurement matrix's columns related to propagation parameters
803             // --------------------------------------------------------------
804 
805             // Jacobian of the measurement with respect to propagation parameters
806             final int nbParams = estimatedPropagationParameters[p].getNbParams();
807             if (nbParams > 0) {
808                 final RealMatrix dYdPp = harvesters[p].getParametersJacobian(evaluationStates[k]);
809                 final RealMatrix dMdPp = dMdY.multiply(dYdPp);
810                 for (int i = 0; i < dMdPp.getRowDimension(); ++i) {
811                     for (int j = 0; j < nbParams; ++j) {
812                         final ParameterDriver delegating = estimatedPropagationParameters[p].getDrivers().get(j);
813                         measurementMatrix.setEntry(i, propagationParameterColumns.get(delegating.getName()),
814                                                    dMdPp.getEntry(i, j) / sigma[i] * delegating.getScale());
815                     }
816                 }
817             }
818 
819             // Normalized measurement matrix's columns related to measurement parameters
820             // --------------------------------------------------------------
821 
822             // Jacobian of the measurement with respect to measurement parameters
823             // Gather the measurement parameters linked to current measurement
824             for (final ParameterDriver driver : observedMeasurement.getParametersDrivers()) {
825                 if (driver.isSelected()) {
826                     // Derivatives of current measurement w/r to selected measurement parameter
827                     final double[] aMPm = predictedMeasurement.getParameterDerivatives(driver);
828 
829                     // Check that the measurement parameter is managed by the filter
830                     if (measurementParameterColumns.get(driver.getName()) != null) {
831                         // Column of the driver in the measurement matrix
832                         final int driverColumn = measurementParameterColumns.get(driver.getName());
833 
834                         // Fill the corresponding indexes of the measurement matrix
835                         for (int i = 0; i < aMPm.length; ++i) {
836                             measurementMatrix.setEntry(i, driverColumn,
837                                                        aMPm[i] / sigma[i] * driver.getScale());
838                         }
839                     }
840                 }
841             }
842         }
843 
844         // Return the normalized measurement matrix
845         return measurementMatrix;
846 
847     }
848 
849     /** Normalize a covariance matrix.
850      * The covariance P is an mxm matrix where m = nbOrb + nbPropag + nbMeas
851      * For each element [i,j] of P the corresponding normalized value is:
852      * Pn[i,j] = P[i,j] / (scale[i]*scale[j])
853      * @param physicalCovarianceMatrix The "physical" covariance matrix in input
854      * @return the normalized covariance matrix
855      */
856     private RealMatrix normalizeCovarianceMatrix(final RealMatrix physicalCovarianceMatrix) {
857 
858         // Initialize output matrix
859         final int nbParams = physicalCovarianceMatrix.getRowDimension();
860         final RealMatrix normalizedCovarianceMatrix = MatrixUtils.createRealMatrix(nbParams, nbParams);
861 
862         // Normalize the state matrix
863         for (int i = 0; i < nbParams; ++i) {
864             for (int j = 0; j < nbParams; ++j) {
865                 normalizedCovarianceMatrix.setEntry(i, j,
866                                                     physicalCovarianceMatrix.getEntry(i, j) /
867                                                     (scale[i] * scale[j]));
868             }
869         }
870         return normalizedCovarianceMatrix;
871     }
872 
873     /** Set and apply a dynamic outlier filter on a measurement.<p>
874      * Loop on the modifiers to see if a dynamic outlier filter needs to be applied.<p>
875      * Compute the sigma array using the matrix in input and set the filter.<p>
876      * Apply the filter by calling the modify method on the estimated measurement.<p>
877      * Reset the filter.
878      * @param measurement measurement to filter
879      * @param innovationCovarianceMatrix So called innovation covariance matrix S, with:<p>
880      *        S = H.Ppred.Ht + R<p>
881      *        Where:<p>
882      *         - H is the normalized measurement matrix (Ht its transpose)<p>
883      *         - Ppred is the normalized predicted covariance matrix<p>
884      *         - R is the normalized measurement noise matrix
885      * @param <T> the type of measurement
886      */
887     private <T extends ObservedMeasurement<T>> void applyDynamicOutlierFilter(final EstimatedMeasurement<T> measurement,
888                                                                               final RealMatrix innovationCovarianceMatrix) {
889 
890         // Observed measurement associated to the predicted measurement
891         final ObservedMeasurement<T> observedMeasurement = measurement.getObservedMeasurement();
892 
893         // Check if a dynamic filter was added to the measurement
894         // If so, update its sigma value and apply it
895         for (EstimationModifier<T> modifier : observedMeasurement.getModifiers()) {
896             if (modifier instanceof DynamicOutlierFilter<?>) {
897                 final DynamicOutlierFilter<T> dynamicOutlierFilter = (DynamicOutlierFilter<T>) modifier;
898 
899                 // Initialize the values of the sigma array used in the dynamic filter
900                 final double[] sigmaDynamic     = new double[innovationCovarianceMatrix.getColumnDimension()];
901                 final double[] sigmaMeasurement = observedMeasurement.getTheoreticalStandardDeviation();
902 
903                 // Set the sigma value for each element of the measurement
904                 // Here we do use the value suggested by David A. Vallado (see [1]ยง10.6):
905                 // sigmaDynamic[i] = sqrt(diag(S))*sigma[i]
906                 // With S = H.Ppred.Ht + R
907                 // Where:
908                 //  - S is the measurement error matrix in input
909                 //  - H is the normalized measurement matrix (Ht its transpose)
910                 //  - Ppred is the normalized predicted covariance matrix
911                 //  - R is the normalized measurement noise matrix
912                 //  - sigma[i] is the theoretical standard deviation of the ith component of the measurement.
913                 //    It is used here to un-normalize the value before it is filtered
914                 for (int i = 0; i < sigmaDynamic.length; i++) {
915                     sigmaDynamic[i] = FastMath.sqrt(innovationCovarianceMatrix.getEntry(i, i)) * sigmaMeasurement[i];
916                 }
917                 dynamicOutlierFilter.setSigma(sigmaDynamic);
918 
919                 // Apply the modifier on the estimated measurement
920                 modifier.modify(measurement);
921 
922                 // Re-initialize the value of the filter for the next measurement of the same type
923                 dynamicOutlierFilter.setSigma(null);
924             }
925         }
926     }
927 
928     /** {@inheritDoc} */
929     @Override
930     public NonLinearEvolution getEvolution(final double previousTime, final RealVector previousState,
931                                            final MeasurementDecorator measurement) {
932 
933         // Set a reference date for all measurements parameters that lack one (including the not estimated ones)
934         final ObservedMeasurement<?> observedMeasurement = measurement.getObservedMeasurement();
935         for (final ParameterDriver driver : observedMeasurement.getParametersDrivers()) {
936             if (driver.getReferenceDate() == null) {
937                 driver.setReferenceDate(builders.get(0).getInitialOrbitDate());
938             }
939         }
940 
941         ++currentMeasurementNumber;
942         currentDate = measurement.getObservedMeasurement().getDate();
943 
944         // Note:
945         // - n = size of the current measurement
946         //  Example:
947         //   * 1 for Range, RangeRate and TurnAroundRange
948         //   * 2 for Angular (Azimuth/Elevation or Right-ascension/Declination)
949         //   * 6 for Position/Velocity
950         // - m = size of the state vector. n = nbOrb + nbPropag + nbMeas
951 
952         // Predict the state vector (mx1)
953         final RealVector predictedState = predictState(observedMeasurement.getDate());
954 
955         // Get the error state transition matrix (mxm)
956         final RealMatrix stateTransitionMatrix = getErrorStateTransitionMatrix();
957 
958         // Predict the measurement based on predicted spacecraft state
959         // Compute the innovations (i.e. residuals of the predicted measurement)
960         // ------------------------------------------------------------
961 
962         // Predicted measurement
963         // Note: here the "iteration/evaluation" formalism from the batch LS method
964         // is twisted to fit the need of the Kalman filter.
965         // The number of "iterations" is actually the number of measurements processed by the filter
966         // so far. We use this to be able to apply the OutlierFilter modifiers on the predicted measurement.
967         predictedMeasurement = observedMeasurement.estimate(currentMeasurementNumber,
968                                                             currentMeasurementNumber,
969                                                             filterRelevant(observedMeasurement, predictedSpacecraftStates));
970 
971         // Normalized measurement matrix (nxm)
972         final RealMatrix measurementMatrix = getMeasurementMatrix();
973 
974         // compute process noise matrix
975         final RealMatrix physicalProcessNoise = MatrixUtils.createRealMatrix(previousState.getDimension(),
976                                                                              previousState.getDimension());
977         for (int k = 0; k < covarianceMatricesProviders.size(); ++k) {
978 
979             // Number of estimated measurement parameters
980             final int nbMeas = estimatedMeasurementsParameters.getNbParams();
981 
982             // Number of estimated dynamic parameters (orbital + propagation)
983             final int nbDyn  = orbitsEndColumns[k] - orbitsStartColumns[k] +
984                                estimatedPropagationParameters[k].getNbParams();
985 
986             // Covariance matrix
987             final RealMatrix noiseK = MatrixUtils.createRealMatrix(nbDyn + nbMeas, nbDyn + nbMeas);
988             final RealMatrix noiseP = covarianceMatricesProviders.get(k).
989                                       getProcessNoiseMatrix(correctedSpacecraftStates[k],
990                                                             predictedSpacecraftStates[k]);
991             noiseK.setSubMatrix(noiseP.getData(), 0, 0);
992             if (measurementProcessNoiseMatrix != null) {
993                 final RealMatrix noiseM = measurementProcessNoiseMatrix.
994                                           getProcessNoiseMatrix(correctedSpacecraftStates[k],
995                                                                 predictedSpacecraftStates[k]);
996                 noiseK.setSubMatrix(noiseM.getData(), nbDyn, nbDyn);
997             }
998 
999             checkDimension(noiseK.getRowDimension(),
1000                            builders.get(k).getOrbitalParametersDrivers(),
1001                            builders.get(k).getPropagationParametersDrivers(),
1002                            estimatedMeasurementsParameters);
1003 
1004             final int[] indK = covarianceIndirection[k];
1005             for (int i = 0; i < indK.length; ++i) {
1006                 if (indK[i] >= 0) {
1007                     for (int j = 0; j < indK.length; ++j) {
1008                         if (indK[j] >= 0) {
1009                             physicalProcessNoise.setEntry(indK[i], indK[j], noiseK.getEntry(i, j));
1010                         }
1011                     }
1012                 }
1013             }
1014 
1015         }
1016         final RealMatrix normalizedProcessNoise = normalizeCovarianceMatrix(physicalProcessNoise);
1017 
1018         return new NonLinearEvolution(measurement.getTime(), predictedState,
1019                                       stateTransitionMatrix, normalizedProcessNoise, measurementMatrix);
1020 
1021     }
1022 
1023 
1024     /** {@inheritDoc} */
1025     @Override
1026     public RealVector getInnovation(final MeasurementDecorator measurement, final NonLinearEvolution evolution,
1027                                     final RealMatrix innovationCovarianceMatrix) {
1028 
1029         // Apply the dynamic outlier filter, if it exists
1030         applyDynamicOutlierFilter(predictedMeasurement, innovationCovarianceMatrix);
1031         if (predictedMeasurement.getStatus() == EstimatedMeasurement.Status.REJECTED)  {
1032             // set innovation to null to notify filter measurement is rejected
1033             return null;
1034         } else {
1035             // Normalized innovation of the measurement (Nx1)
1036             final double[] observed  = predictedMeasurement.getObservedMeasurement().getObservedValue();
1037             final double[] estimated = predictedMeasurement.getEstimatedValue();
1038             final double[] sigma     = predictedMeasurement.getObservedMeasurement().getTheoreticalStandardDeviation();
1039             final double[] residuals = new double[observed.length];
1040 
1041             for (int i = 0; i < observed.length; i++) {
1042                 residuals[i] = (observed[i] - estimated[i]) / sigma[i];
1043             }
1044             return MatrixUtils.createRealVector(residuals);
1045         }
1046     }
1047 
1048     /** Finalize estimation.
1049      * @param observedMeasurement measurement that has just been processed
1050      * @param estimate corrected estimate
1051      */
1052     public void finalizeEstimation(final ObservedMeasurement<?> observedMeasurement,
1053                                    final ProcessEstimate estimate) {
1054         // Update the parameters with the estimated state
1055         // The min/max values of the parameters are handled by the ParameterDriver implementation
1056         correctedEstimate = estimate;
1057         updateParameters();
1058 
1059         // Get the estimated propagator (mirroring parameter update in the builder)
1060         // and the estimated spacecraft state
1061         final Propagator[] estimatedPropagators = getEstimatedPropagators();
1062         for (int k = 0; k < estimatedPropagators.length; ++k) {
1063             correctedSpacecraftStates[k] = estimatedPropagators[k].getInitialState();
1064         }
1065 
1066         // Compute the estimated measurement using estimated spacecraft state
1067         correctedMeasurement = observedMeasurement.estimate(currentMeasurementNumber,
1068                                                             currentMeasurementNumber,
1069                                                             filterRelevant(observedMeasurement, correctedSpacecraftStates));
1070         // Update the trajectory
1071         // ---------------------
1072         updateReferenceTrajectories(estimatedPropagators, propagationType, stateType);
1073 
1074     }
1075 
1076     /** Filter relevant states for a measurement.
1077      * @param observedMeasurement measurement to consider
1078      * @param allStates all states
1079      * @return array containing only the states relevant to the measurement
1080      * @since 10.1
1081      */
1082     private SpacecraftState[] filterRelevant(final ObservedMeasurement<?> observedMeasurement, final SpacecraftState[] allStates) {
1083         final List<ObservableSatellite> satellites = observedMeasurement.getSatellites();
1084         final SpacecraftState[] relevantStates = new SpacecraftState[satellites.size()];
1085         for (int i = 0; i < relevantStates.length; ++i) {
1086             relevantStates[i] = allStates[satellites.get(i).getPropagatorIndex()];
1087         }
1088         return relevantStates;
1089     }
1090 
1091     /** Set the predicted normalized state vector.
1092      * The predicted/propagated orbit is used to update the state vector
1093      * @param date prediction date
1094      * @return predicted state
1095      */
1096     private RealVector predictState(final AbsoluteDate date) {
1097 
1098         // Predicted state is initialized to previous estimated state
1099         final RealVector predictedState = correctedEstimate.getState().copy();
1100 
1101         // Orbital parameters counter
1102         int jOrb = 0;
1103 
1104         for (int k = 0; k < predictedSpacecraftStates.length; ++k) {
1105 
1106             // Propagate the reference trajectory to measurement date
1107             predictedSpacecraftStates[k] = referenceTrajectories[k].propagate(date);
1108 
1109             // Update the builder with the predicted orbit
1110             // This updates the orbital drivers with the values of the predicted orbit
1111             builders.get(k).resetOrbit(predictedSpacecraftStates[k].getOrbit());
1112 
1113             // The orbital parameters in the state vector are replaced with their predicted values
1114             // The propagation & measurement parameters are not changed by the prediction (i.e. the propagation)
1115             // As the propagator builder was previously updated with the predicted orbit,
1116             // the selected orbital drivers are already up to date with the prediction
1117             for (DelegatingDriver orbitalDriver : builders.get(k).getOrbitalParametersDrivers().getDrivers()) {
1118                 if (orbitalDriver.isSelected()) {
1119                     predictedState.setEntry(jOrb++, orbitalDriver.getNormalizedValue());
1120                 }
1121             }
1122 
1123         }
1124 
1125         return predictedState;
1126 
1127     }
1128 
1129     /** Update the estimated parameters after the correction phase of the filter.
1130      * The min/max allowed values are handled by the parameter themselves.
1131      */
1132     private void updateParameters() {
1133         final RealVector correctedState = correctedEstimate.getState();
1134         int i = 0;
1135         for (final DelegatingDriver driver : getEstimatedOrbitalParameters().getDrivers()) {
1136             // let the parameter handle min/max clipping
1137             driver.setNormalizedValue(correctedState.getEntry(i));
1138             correctedState.setEntry(i++, driver.getNormalizedValue());
1139         }
1140         for (final DelegatingDriver driver : getEstimatedPropagationParameters().getDrivers()) {
1141             // let the parameter handle min/max clipping
1142             driver.setNormalizedValue(correctedState.getEntry(i));
1143             correctedState.setEntry(i++, driver.getNormalizedValue());
1144         }
1145         for (final DelegatingDriver driver : getEstimatedMeasurementsParameters().getDrivers()) {
1146             // let the parameter handle min/max clipping
1147             driver.setNormalizedValue(correctedState.getEntry(i));
1148             correctedState.setEntry(i++, driver.getNormalizedValue());
1149         }
1150     }
1151 
1152     /** Getter for the propagators.
1153      * @return the propagators
1154      */
1155     public List<OrbitDeterminationPropagatorBuilder> getBuilders() {
1156         return builders;
1157     }
1158 
1159     /** Getter for the reference trajectories.
1160      * @return the referencetrajectories
1161      */
1162     public Propagator[] getReferenceTrajectories() {
1163         return referenceTrajectories.clone();
1164     }
1165 
1166     /** Setter for the reference trajectories.
1167      * @param referenceTrajectories the reference trajectories to be setted
1168      */
1169     public void setReferenceTrajectories(final Propagator[] referenceTrajectories) {
1170         this.referenceTrajectories = referenceTrajectories.clone();
1171     }
1172 
1173     /** Getter for the jacobian mappers.
1174      * @return the jacobian mappers
1175      * @deprecated as of 11.1, not used anymore
1176      */
1177     @Deprecated
1178     public AbstractJacobiansMapper[] getMappers() {
1179         return null;
1180     }
1181 
1182     /** Setter for the jacobian mappers.
1183      * @param mappers the jacobian mappers to set
1184      * @deprecated as of 11.1, replaced by {@link #setHarvesters(MatricesHarvester[])}
1185      */
1186     @Deprecated
1187     public void setMappers(final AbstractJacobiansMapper[] mappers) {
1188         setHarvesters(mappers);
1189     }
1190 
1191     /** Setter for the jacobian harvesters.
1192      * @param harvesters the jacobian harvesters to set
1193      * @since 11.1
1194      */
1195     public void setHarvesters(final MatricesHarvester[] harvesters) {
1196         this.harvesters = harvesters.clone();
1197     }
1198 
1199     /** Get the propagators estimated with the values set in the propagators builders.
1200      * @return propagators based on the current values in the builder
1201      */
1202     public Propagator[] getEstimatedPropagators() {
1203         // Return propagators built with current instantiation of the propagator builders
1204         final Propagator[] propagators = new Propagator[getBuilders().size()];
1205         for (int k = 0; k < getBuilders().size(); ++k) {
1206             propagators[k] = getBuilders().get(k).buildPropagator(getBuilders().get(k).getSelectedNormalizedParameters());
1207         }
1208         return propagators;
1209     }
1210 
1211 }