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.estimation.sequential;
18  
19  import org.hipparchus.filtering.kalman.ProcessEstimate;
20  import org.hipparchus.filtering.kalman.unscented.UnscentedEvolution;
21  import org.hipparchus.filtering.kalman.unscented.UnscentedKalmanFilter;
22  import org.hipparchus.filtering.kalman.unscented.UnscentedProcess;
23  import org.hipparchus.linear.ArrayRealVector;
24  import org.hipparchus.linear.MatrixUtils;
25  import org.hipparchus.linear.RealMatrix;
26  import org.hipparchus.linear.RealVector;
27  import org.hipparchus.util.FastMath;
28  import org.orekit.estimation.measurements.EstimatedMeasurement;
29  import org.orekit.estimation.measurements.EstimatedMeasurementBase;
30  import org.orekit.estimation.measurements.ObservedMeasurement;
31  import org.orekit.orbits.Orbit;
32  import org.orekit.orbits.OrbitType;
33  import org.orekit.orbits.PositionAngleType;
34  import org.orekit.propagation.PropagationType;
35  import org.orekit.propagation.SpacecraftState;
36  import org.orekit.propagation.conversion.DSSTPropagatorBuilder;
37  import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
38  import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
39  import org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms;
40  import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
41  import org.orekit.time.AbsoluteDate;
42  import org.orekit.time.ChronologicalComparator;
43  import org.orekit.utils.ParameterDriver;
44  import org.orekit.utils.ParameterDriversList;
45  import org.orekit.utils.ParameterDriversList.DelegatingDriver;
46  
47  import java.util.ArrayList;
48  import java.util.Comparator;
49  import java.util.List;
50  
51  /** Class defining the process model dynamics to use with a {@link SemiAnalyticalUnscentedKalmanEstimator}.
52   * @author Gaëtan Pierre
53   * @author Bryan Cazabonne
54   * @since 11.3
55   */
56  public class SemiAnalyticalUnscentedKalmanModel implements KalmanEstimation, UnscentedProcess<MeasurementDecorator>, SemiAnalyticalProcess {
57  
58      /** Initial builder for propagator. */
59      private final DSSTPropagatorBuilder builder;
60  
61      /** Estimated orbital parameters. */
62      private final ParameterDriversList estimatedOrbitalParameters;
63  
64      /** Estimated propagation parameters. */
65      private final ParameterDriversList estimatedPropagationParameters;
66  
67      /** Estimated measurements parameters. */
68      private final ParameterDriversList estimatedMeasurementsParameters;
69  
70      /** Provider for covariance matrice. */
71      private final CovarianceMatrixProvider covarianceMatrixProvider;
72  
73      /** Process noise matrix provider for measurement parameters. */
74      private final CovarianceMatrixProvider measurementProcessNoiseMatrix;
75  
76      /** Position angle type used during orbit determination. */
77      private final PositionAngleType angleType;
78  
79      /** Orbit type used during orbit determination. */
80      private final OrbitType orbitType;
81  
82      /** Current corrected estimate. */
83      private ProcessEstimate correctedEstimate;
84  
85      /** Observer to retrieve current estimation info. */
86      private KalmanObserver observer;
87  
88      /** Current number of measurement. */
89      private int currentMeasurementNumber;
90  
91      /** Current date. */
92      private AbsoluteDate currentDate;
93  
94      /** Nominal mean spacecraft state. */
95      private SpacecraftState nominalMeanSpacecraftState;
96  
97      /** Previous nominal mean spacecraft state. */
98      private SpacecraftState previousNominalMeanSpacecraftState;
99  
100     /** Predicted osculating spacecraft state. */
101     private SpacecraftState predictedSpacecraftState;
102 
103     /** Corrected mean spacecraft state. */
104     private SpacecraftState correctedSpacecraftState;
105 
106     /** Predicted measurement. */
107     private EstimatedMeasurement<?> predictedMeasurement;
108 
109     /** Corrected measurement. */
110     private EstimatedMeasurement<?> correctedMeasurement;
111 
112     /** Predicted mean element filter correction. */
113     private RealVector predictedFilterCorrection;
114 
115     /** Corrected mean element filter correction. */
116     private RealVector correctedFilterCorrection;
117 
118     /** Propagators for the reference trajectories, up to current date. */
119     private final DSSTPropagator dsstPropagator;
120 
121     /** Short period terms for the nominal mean spacecraft state. */
122     private RealVector shortPeriodicTerms;
123 
124     /** Unscented Kalman process model constructor (package private).
125      * @param propagatorBuilder propagators builders used to evaluate the orbits.
126      * @param covarianceMatrixProvider provider for covariance matrix
127      * @param estimatedMeasurementParameters measurement parameters to estimate
128      * @param measurementProcessNoiseMatrix provider for measurement process noise matrix
129      */
130     protected SemiAnalyticalUnscentedKalmanModel(final DSSTPropagatorBuilder propagatorBuilder,
131                                                  final CovarianceMatrixProvider covarianceMatrixProvider,
132                                                  final ParameterDriversList estimatedMeasurementParameters,
133                                                  final CovarianceMatrixProvider measurementProcessNoiseMatrix) {
134 
135         this.builder                         = propagatorBuilder;
136         this.angleType                       = propagatorBuilder.getPositionAngleType();
137         this.orbitType                       = propagatorBuilder.getOrbitType();
138         this.estimatedMeasurementsParameters = estimatedMeasurementParameters;
139         this.currentMeasurementNumber        = 0;
140         this.currentDate                     = propagatorBuilder.getInitialOrbitDate();
141         this.covarianceMatrixProvider        = covarianceMatrixProvider;
142         this.measurementProcessNoiseMatrix   = measurementProcessNoiseMatrix;
143 
144         // Number of estimated parameters
145         int columns = 0;
146 
147         // Set estimated orbital parameters
148         this.estimatedOrbitalParameters = new ParameterDriversList();
149         for (final ParameterDriver driver : propagatorBuilder.getOrbitalParametersDrivers().getDrivers()) {
150 
151             // Verify if the driver reference date has been set
152             if (driver.getReferenceDate() == null) {
153                 driver.setReferenceDate(currentDate);
154             }
155 
156             // Verify if the driver is selected
157             if (driver.isSelected()) {
158                 estimatedOrbitalParameters.add(driver);
159                 columns++;
160             }
161 
162         }
163 
164         // Set estimated propagation parameters
165         this.estimatedPropagationParameters = new ParameterDriversList();
166         final List<String> estimatedPropagationParametersNames = new ArrayList<>();
167         for (final ParameterDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
168 
169             // Verify if the driver reference date has been set
170             if (driver.getReferenceDate() == null) {
171                 driver.setReferenceDate(currentDate);
172             }
173 
174             // Verify if the driver is selected
175             if (driver.isSelected()) {
176                 estimatedPropagationParameters.add(driver);
177                 final String driverName = driver.getName();
178                 // Add the driver name if it has not been added yet
179                 if (!estimatedPropagationParametersNames.contains(driverName)) {
180                     estimatedPropagationParametersNames.add(driverName);
181                     ++columns;
182                 }
183             }
184 
185         }
186         estimatedPropagationParametersNames.sort(Comparator.naturalOrder());
187 
188         // Set the estimated measurement parameters
189         for (final ParameterDriver parameter : estimatedMeasurementsParameters.getDrivers()) {
190             if (parameter.getReferenceDate() == null) {
191                 parameter.setReferenceDate(currentDate);
192             }
193             ++columns;
194         }
195 
196         // Number of estimated measurement parameters
197         final int nbMeas = estimatedMeasurementParameters.getNbParams();
198 
199         // Number of estimated dynamic parameters (orbital + propagation)
200         final int nbDyn  = estimatedOrbitalParameters.getNbParams() +
201                            estimatedPropagationParameters.getNbParams();
202 
203         // Build the reference propagator
204         this.dsstPropagator = getEstimatedPropagator();
205         final SpacecraftState meanState = dsstPropagator.initialIsOsculating() ?
206                          DSSTPropagator.computeMeanState(dsstPropagator.getInitialState(), dsstPropagator.getAttitudeProvider(), dsstPropagator.getAllForceModels()) :
207                          dsstPropagator.getInitialState();
208         this.nominalMeanSpacecraftState         = meanState;
209         this.predictedSpacecraftState           = meanState;
210         this.correctedSpacecraftState           = predictedSpacecraftState;
211         this.previousNominalMeanSpacecraftState = nominalMeanSpacecraftState;
212 
213         // Initialize the estimated mean element filter correction
214         this.predictedFilterCorrection = MatrixUtils.createRealVector(columns);
215         this.correctedFilterCorrection = predictedFilterCorrection;
216 
217         // Covariance matrix
218         final RealMatrix noiseK = MatrixUtils.createRealMatrix(nbDyn + nbMeas, nbDyn + nbMeas);
219         final RealMatrix noiseP = covarianceMatrixProvider.getInitialCovarianceMatrix(nominalMeanSpacecraftState);
220         noiseK.setSubMatrix(noiseP.getData(), 0, 0);
221         if (measurementProcessNoiseMatrix != null) {
222             final RealMatrix noiseM = measurementProcessNoiseMatrix.getInitialCovarianceMatrix(nominalMeanSpacecraftState);
223             noiseK.setSubMatrix(noiseM.getData(), nbDyn, nbDyn);
224         }
225 
226         KalmanEstimatorUtil.checkDimension(noiseK.getRowDimension(),
227                                            propagatorBuilder.getOrbitalParametersDrivers(),
228                                            propagatorBuilder.getPropagationParametersDrivers(),
229                                            estimatedMeasurementsParameters);
230 
231         // Initialize corrected estimate
232         this.correctedEstimate = new ProcessEstimate(0.0, correctedFilterCorrection, noiseK);
233 
234     }
235 
236     /** {@inheritDoc} */
237     @Override
238     public KalmanObserver getObserver() {
239         return observer;
240     }
241 
242     /** Set the observer.
243      * @param observer the observer
244      */
245     public void setObserver(final KalmanObserver observer) {
246         this.observer = observer;
247     }
248 
249     /** Get the current corrected estimate.
250      * <p>
251      * For the Unscented Semi-analytical Kalman Filter
252      * it corresponds to the corrected filter correction.
253      * In other words, it doesn't represent an orbital state.
254      * </p>
255      * @return current corrected estimate
256      */
257     public ProcessEstimate getEstimate() {
258         return correctedEstimate;
259     }
260 
261     /** Process measurements.
262      * @param observedMeasurements the list of measurements to process
263      * @param filter Unscented Kalman Filter
264      * @return estimated propagator
265      */
266     public DSSTPropagator processMeasurements(final List<ObservedMeasurement<?>> observedMeasurements,
267                                               final UnscentedKalmanFilter<MeasurementDecorator> filter) {
268 
269         // Sort the measurement
270         observedMeasurements.sort(new ChronologicalComparator());
271         final AbsoluteDate tStart             = observedMeasurements.get(0).getDate();
272         final AbsoluteDate tEnd               = observedMeasurements.get(observedMeasurements.size() - 1).getDate();
273         final double       overshootTimeRange = FastMath.nextAfter(tEnd.durationFrom(tStart),
274                                                 Double.POSITIVE_INFINITY);
275 
276         // Initialize step handler and set it to a parallelized propagator
277         final SemiAnalyticalMeasurementHandler  stepHandler = new SemiAnalyticalMeasurementHandler(this, filter, observedMeasurements, builder.getInitialOrbitDate(), true);
278         dsstPropagator.getMultiplexer().add(stepHandler);
279         dsstPropagator.propagate(tStart, tStart.shiftedBy(overshootTimeRange));
280 
281         // Return the last estimated propagator
282         return getEstimatedPropagator();
283 
284     }
285 
286     /** Get the propagator estimated with the values set in the propagator builder.
287      * @return propagator based on the current values in the builder
288      */
289     public DSSTPropagator getEstimatedPropagator() {
290         // Return propagator built with current instantiation of the propagator builder
291         return (DSSTPropagator) builder.buildPropagator();
292     }
293 
294     /** {@inheritDoc} */
295     @Override
296     public UnscentedEvolution getEvolution(final double previousTime, final RealVector[] sigmaPoints,
297                                            final MeasurementDecorator measurement) {
298 
299         // Set a reference date for all measurements parameters that lack one (including the not estimated ones)
300         final ObservedMeasurement<?> observedMeasurement = measurement.getObservedMeasurement();
301         for (final ParameterDriver driver : observedMeasurement.getParametersDrivers()) {
302             if (driver.getReferenceDate() == null) {
303                 driver.setReferenceDate(builder.getInitialOrbitDate());
304             }
305         }
306 
307         // Increment measurement number
308         ++currentMeasurementNumber;
309 
310         // Update the current date
311         currentDate = measurement.getObservedMeasurement().getDate();
312 
313         // STM for the prediction of the filter correction
314         final RealMatrix stm = getStm();
315 
316         // Predicted states
317         final RealVector[] predictedStates = new RealVector[sigmaPoints.length];
318         for (int k = 0; k < sigmaPoints.length; ++k) {
319             // Predict filter correction for the current sigma point
320             final RealVector predicted = stm.operate(sigmaPoints[k]);
321             predictedStates[k] = predicted;
322         }
323 
324         // Return
325         return new UnscentedEvolution(measurement.getTime(), predictedStates);
326 
327     }
328 
329     /** {@inheritDoc} */
330     @Override
331     public RealMatrix getProcessNoiseMatrix(final double previousTime, final RealVector predictedState,
332                                             final MeasurementDecorator measurement) {
333 
334         // Number of estimated measurement parameters
335         final int nbMeas = getNumberSelectedMeasurementDrivers();
336 
337         // Number of estimated dynamic parameters (orbital + propagation)
338         final int nbDyn  = getNumberSelectedOrbitalDrivers() + getNumberSelectedPropagationDrivers();
339 
340         // Covariance matrix
341         final RealMatrix noiseK = MatrixUtils.createRealMatrix(nbDyn + nbMeas, nbDyn + nbMeas);
342         final RealMatrix noiseP = covarianceMatrixProvider.getProcessNoiseMatrix(previousNominalMeanSpacecraftState, nominalMeanSpacecraftState);
343         noiseK.setSubMatrix(noiseP.getData(), 0, 0);
344         if (measurementProcessNoiseMatrix != null) {
345             final RealMatrix noiseM = measurementProcessNoiseMatrix.getProcessNoiseMatrix(previousNominalMeanSpacecraftState, nominalMeanSpacecraftState);
346             noiseK.setSubMatrix(noiseM.getData(), nbDyn, nbDyn);
347         }
348 
349         // Verify dimension
350         KalmanEstimatorUtil.checkDimension(noiseK.getRowDimension(),
351                 builder.getOrbitalParametersDrivers(),
352                 builder.getPropagationParametersDrivers(),
353                 estimatedMeasurementsParameters);
354 
355         return noiseK;
356     }
357 
358     /** {@inheritDoc} */
359     @Override
360     public RealVector[] getPredictedMeasurements(final RealVector[] predictedSigmaPoints, final MeasurementDecorator measurement) {
361 
362         // Observed measurement
363         final ObservedMeasurement<?> observedMeasurement = measurement.getObservedMeasurement();
364 
365         // Initialize arrays of predicted states and measurements
366         final RealVector[] predictedMeasurements = new RealVector[predictedSigmaPoints.length];
367 
368         // Loop on sigma points
369         for (int k = 0; k < predictedSigmaPoints.length; ++k) {
370 
371             // Calculate the predicted osculating elements for the current mean state
372             final RealVector osculating = computeOsculatingElements(predictedSigmaPoints[k], nominalMeanSpacecraftState, shortPeriodicTerms);
373             final Orbit osculatingOrbit = orbitType.mapArrayToOrbit(osculating.toArray(), null, angleType,
374                                                                     currentDate, builder.getMu(), builder.getFrame());
375 
376             // Then, estimate the measurement
377             final EstimatedMeasurement<?> estimated = estimateMeasurement(observedMeasurement, currentMeasurementNumber,
378                 new SpacecraftState[] { new SpacecraftState(osculatingOrbit) });
379             predictedMeasurements[k] = new ArrayRealVector(estimated.getEstimatedValue());
380 
381         }
382 
383         // Return
384         return predictedMeasurements;
385 
386     }
387 
388     /** {@inheritDoc} */
389     @Override
390     public RealVector getInnovation(final MeasurementDecorator measurement, final RealVector predictedMeas,
391                                     final RealVector predictedState, final RealMatrix innovationCovarianceMatrix) {
392 
393         // Predicted filter correction
394         predictedFilterCorrection = predictedState;
395 
396         // Predicted measurement
397         final RealVector osculating = computeOsculatingElements(predictedFilterCorrection, nominalMeanSpacecraftState, shortPeriodicTerms);
398         final Orbit osculatingOrbit = orbitType.mapArrayToOrbit(osculating.toArray(), null, angleType,
399                                                                 currentDate, builder.getMu(), builder.getFrame());
400         predictedSpacecraftState = new SpacecraftState(osculatingOrbit);
401         predictedMeasurement = estimateMeasurement(measurement.getObservedMeasurement(), currentMeasurementNumber,
402             getPredictedSpacecraftStates());
403         predictedMeasurement.setEstimatedValue(predictedMeas.toArray());
404 
405         // Apply the dynamic outlier filter, if it exists
406         KalmanEstimatorUtil.applyDynamicOutlierFilter(predictedMeasurement, innovationCovarianceMatrix);
407 
408         // Compute the innovation vector (not normalized for unscented Kalman filter)
409         return KalmanEstimatorUtil.computeInnovationVector(predictedMeasurement);
410 
411     }
412 
413 
414     /** {@inheritDoc} */
415     @Override
416     public void finalizeEstimation(final ObservedMeasurement<?> observedMeasurement,
417                                    final ProcessEstimate estimate) {
418         // Update the process estimate
419         correctedEstimate = estimate;
420         // Corrected filter correction
421         correctedFilterCorrection = estimate.getState();
422 
423         // Update the previous nominal mean spacecraft state
424         previousNominalMeanSpacecraftState = nominalMeanSpacecraftState;
425 
426         // Update the previous nominal mean spacecraft state
427         // Calculate the corrected osculating elements
428         final RealVector osculating = computeOsculatingElements(correctedFilterCorrection, nominalMeanSpacecraftState, shortPeriodicTerms);
429         final Orbit osculatingOrbit = orbitType.mapArrayToOrbit(osculating.toArray(), null, builder.getPositionAngleType(),
430                                                                 currentDate, builder.getMu(), builder.getFrame());
431 
432         // Compute the corrected measurements
433         correctedSpacecraftState = new SpacecraftState(osculatingOrbit);
434         correctedMeasurement = estimateMeasurement(observedMeasurement, currentMeasurementNumber,
435             getCorrectedSpacecraftStates());
436 
437         // Call the observer if the user add one
438         if (observer != null) {
439             observer.evaluationPerformed(this);
440         }
441     }
442 
443     /**
444      * Estimate measurement (without derivatives).
445      * @param <T> measurement type
446      * @param observedMeasurement observed measurement
447      * @param measurementNumber measurement number
448      * @param spacecraftStates states
449      * @return estimated measurements
450      * @since 12.1
451      */
452     private static <T extends ObservedMeasurement<T>> EstimatedMeasurement<?> estimateMeasurement(final ObservedMeasurement<T> observedMeasurement,
453                                                                                                   final int measurementNumber,
454                                                                                                   final SpacecraftState[] spacecraftStates) {
455         final EstimatedMeasurementBase<T> estimatedMeasurementBase = observedMeasurement.
456                 estimateWithoutDerivatives(measurementNumber, measurementNumber,
457                         KalmanEstimatorUtil.filterRelevant(observedMeasurement, spacecraftStates));
458         final EstimatedMeasurement<T> estimatedMeasurement = new EstimatedMeasurement<>(estimatedMeasurementBase.getObservedMeasurement(),
459                 estimatedMeasurementBase.getIteration(), estimatedMeasurementBase.getCount(),
460                 estimatedMeasurementBase.getStates(), estimatedMeasurementBase.getParticipants());
461         estimatedMeasurement.setEstimatedValue(estimatedMeasurementBase.getEstimatedValue());
462         return estimatedMeasurement;
463     }
464 
465     /** Get the state transition matrix used to predict the filter correction.
466      * <p>
467      * The state transition matrix is not computed by the DSST propagator.
468      * It is analytically calculated considering Keplerian contribution only
469      * </p>
470      * @return the state transition matrix used to predict the filter correction
471      */
472     private RealMatrix getStm() {
473 
474         // initialize the STM
475         final int nbDym  = getNumberSelectedOrbitalDrivers() + getNumberSelectedPropagationDrivers();
476         final int nbMeas = getNumberSelectedMeasurementDrivers();
477         final RealMatrix stm = MatrixUtils.createRealIdentityMatrix(nbDym + nbMeas);
478 
479         // State transition matrix using Keplerian contribution only
480         final double mu  = builder.getMu();
481         final double sma = previousNominalMeanSpacecraftState.getOrbit().getA();
482         final double dt  = currentDate.durationFrom(previousNominalMeanSpacecraftState.getDate());
483         final double contribution = -1.5 * dt * FastMath.sqrt(mu / FastMath.pow(sma, 5));
484         stm.setEntry(5, 0, contribution);
485 
486         // Return
487         return stm;
488 
489     }
490 
491     /** {@inheritDoc} */
492     @Override
493     public void finalizeOperationsObservationGrid() {
494         // Update parameters
495         updateParameters();
496     }
497 
498     /** {@inheritDoc} */
499     @Override
500     public ParameterDriversList getEstimatedOrbitalParameters() {
501         return estimatedOrbitalParameters;
502     }
503 
504     /** {@inheritDoc} */
505     @Override
506     public ParameterDriversList getEstimatedPropagationParameters() {
507         return estimatedPropagationParameters;
508     }
509 
510     /** {@inheritDoc} */
511     @Override
512     public ParameterDriversList getEstimatedMeasurementsParameters() {
513         return estimatedMeasurementsParameters;
514     }
515 
516     /** {@inheritDoc}
517      * <p>
518      * Predicted state is osculating.
519      * </p>
520      */
521     @Override
522     public SpacecraftState[] getPredictedSpacecraftStates() {
523         return new SpacecraftState[] {predictedSpacecraftState};
524     }
525 
526     /** {@inheritDoc}
527      * <p>
528      * Corrected state is osculating.
529      * </p>
530      */
531     @Override
532     public SpacecraftState[] getCorrectedSpacecraftStates() {
533         return new SpacecraftState[] {correctedSpacecraftState};
534     }
535 
536     /** {@inheritDoc} */
537     @Override
538     public RealVector getPhysicalEstimatedState() {
539         // Method {@link ParameterDriver#getValue()} is used to get
540         // the physical values of the state.
541         // The scales'array is used to get the size of the state vector
542         final RealVector physicalEstimatedState = new ArrayRealVector(getEstimate().getState().getDimension());
543         int i = 0;
544         for (final DelegatingDriver driver : getEstimatedOrbitalParameters().getDrivers()) {
545             physicalEstimatedState.setEntry(i++, driver.getValue());
546         }
547         for (final DelegatingDriver driver : getEstimatedPropagationParameters().getDrivers()) {
548             physicalEstimatedState.setEntry(i++, driver.getValue());
549         }
550         for (final DelegatingDriver driver : getEstimatedMeasurementsParameters().getDrivers()) {
551             physicalEstimatedState.setEntry(i++, driver.getValue());
552         }
553 
554         return physicalEstimatedState;
555     }
556 
557     /** {@inheritDoc} */
558     @Override
559     public RealMatrix getPhysicalEstimatedCovarianceMatrix() {
560         return correctedEstimate.getCovariance();
561     }
562 
563     /** {@inheritDoc} */
564     @Override
565     public RealMatrix getPhysicalStateTransitionMatrix() {
566         return null;
567     }
568 
569     /** {@inheritDoc} */
570     @Override
571     public RealMatrix getPhysicalMeasurementJacobian() {
572         return null;
573     }
574 
575     /** {@inheritDoc} */
576     @Override
577     public RealMatrix getPhysicalInnovationCovarianceMatrix() {
578         return correctedEstimate.getInnovationCovariance();
579     }
580 
581     /** {@inheritDoc} */
582     @Override
583     public RealMatrix getPhysicalKalmanGain() {
584         return correctedEstimate.getKalmanGain();
585     }
586 
587     /** {@inheritDoc} */
588     @Override
589     public int getCurrentMeasurementNumber() {
590         return currentMeasurementNumber;
591     }
592 
593     /** {@inheritDoc} */
594     @Override
595     public AbsoluteDate getCurrentDate() {
596         return currentDate;
597     }
598 
599     /** {@inheritDoc} */
600     @Override
601     public EstimatedMeasurement<?> getPredictedMeasurement() {
602         return predictedMeasurement;
603     }
604 
605     /** {@inheritDoc} */
606     @Override
607     public EstimatedMeasurement<?> getCorrectedMeasurement() {
608         return correctedMeasurement;
609     }
610 
611     /** {@inheritDoc} */
612     @Override
613     public void updateNominalSpacecraftState(final SpacecraftState nominal) {
614         this.nominalMeanSpacecraftState = nominal;
615         // Short period terms
616         shortPeriodicTerms = new ArrayRealVector(dsstPropagator.getShortPeriodTermsValue(nominalMeanSpacecraftState));
617         // Update the builder with the nominal mean elements orbit
618         builder.resetOrbit(nominal.getOrbit(), PropagationType.MEAN);
619     }
620 
621     /** {@inheritDoc} */
622     @Override
623     public void updateShortPeriods(final SpacecraftState state) {
624         // Loop on DSST force models
625         for (final DSSTForceModel model : dsstPropagator.getAllForceModels()) {
626             model.updateShortPeriodTerms(model.getParameters(), state);
627         }
628     }
629 
630     /** {@inheritDoc} */
631     @Override
632     public void initializeShortPeriodicTerms(final SpacecraftState meanState) {
633         final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<>();
634         for (final DSSTForceModel force :  builder.getAllForceModels()) {
635             shortPeriodTerms.addAll(force.initializeShortPeriodTerms(new AuxiliaryElements(meanState.getOrbit(), 1), PropagationType.OSCULATING, force.getParameters()));
636         }
637         dsstPropagator.setShortPeriodTerms(shortPeriodTerms);
638     }
639 
640     /** Compute the predicted osculating elements.
641      * @param filterCorrection physical kalman filter correction
642      * @param meanState mean spacecraft state
643      * @param shortPeriodTerms short period terms for the given mean state
644      * @return the predicted osculating element
645      */
646     private RealVector computeOsculatingElements(final RealVector filterCorrection,
647                                                  final SpacecraftState meanState,
648                                                  final RealVector shortPeriodTerms) {
649 
650         // Convert the input predicted mean state to a SpacecraftState
651         final RealVector stateVector = toRealVector(meanState);
652 
653         // Return
654         return stateVector.add(filterCorrection).add(shortPeriodTerms);
655 
656     }
657 
658     /** Convert a SpacecraftState to a RealVector.
659      * @param state the input SpacecraftState
660      * @return the corresponding RealVector
661      */
662     private RealVector toRealVector(final SpacecraftState state) {
663 
664         // Convert orbit to array
665         final double[] stateArray = new double[6];
666         orbitType.mapOrbitToArray(state.getOrbit(), angleType, stateArray, null);
667 
668         // Return the RealVector
669         return new ArrayRealVector(stateArray);
670 
671     }
672 
673     /** Get the number of estimated orbital parameters.
674      * @return the number of estimated orbital parameters
675      */
676     public int getNumberSelectedOrbitalDrivers() {
677         return estimatedOrbitalParameters.getNbParams();
678     }
679 
680     /** Get the number of estimated propagation parameters.
681      * @return the number of estimated propagation parameters
682      */
683     public int getNumberSelectedPropagationDrivers() {
684         return estimatedPropagationParameters.getNbParams();
685     }
686 
687     /** Get the number of estimated measurement parameters.
688      * @return the number of estimated measurement parameters
689      */
690     public int getNumberSelectedMeasurementDrivers() {
691         return estimatedMeasurementsParameters.getNbParams();
692     }
693 
694     /** Update the estimated parameters after the correction phase of the filter.
695      * The min/max allowed values are handled by the parameter themselves.
696      */
697     private void updateParameters() {
698         final RealVector correctedState = correctedEstimate.getState();
699         int i = 0;
700         for (final DelegatingDriver driver : getEstimatedOrbitalParameters().getDrivers()) {
701             // let the parameter handle min/max clipping
702             driver.setValue(driver.getValue() + correctedState.getEntry(i++));
703         }
704         for (final DelegatingDriver driver : getEstimatedPropagationParameters().getDrivers()) {
705             // let the parameter handle min/max clipping
706             driver.setValue(driver.getValue() + correctedState.getEntry(i++));
707         }
708         for (final DelegatingDriver driver : getEstimatedMeasurementsParameters().getDrivers()) {
709             // let the parameter handle min/max clipping
710             driver.setValue(driver.getValue() + correctedState.getEntry(i++));
711         }
712     }
713 
714 }