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.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.linear.LUDecomposition;
21  import org.hipparchus.linear.MatrixUtils;
22  import org.hipparchus.linear.RealMatrix;
23  import org.hipparchus.linear.RealVector;
24  import org.junit.jupiter.api.Assertions;
25  import org.junit.jupiter.api.BeforeEach;
26  import org.junit.jupiter.api.Test;
27  import org.orekit.estimation.Context;
28  import org.orekit.estimation.EstimationTestUtils;
29  import org.orekit.estimation.Force;
30  import org.orekit.estimation.measurements.EstimatedMeasurement;
31  import org.orekit.estimation.measurements.GroundStation;
32  import org.orekit.estimation.measurements.ObservableSatellite;
33  import org.orekit.estimation.measurements.ObservedMeasurement;
34  import org.orekit.estimation.measurements.PV;
35  import org.orekit.estimation.measurements.Range;
36  import org.orekit.estimation.measurements.modifiers.Bias;
37  import org.orekit.forces.maneuvers.ConstantThrustManeuver;
38  import org.orekit.forces.radiation.RadiationSensitive;
39  import org.orekit.orbits.Orbit;
40  import org.orekit.orbits.OrbitType;
41  import org.orekit.orbits.PositionAngleType;
42  import org.orekit.propagation.MatricesHarvester;
43  import org.orekit.propagation.Propagator;
44  import org.orekit.propagation.SpacecraftState;
45  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
46  import org.orekit.time.AbsoluteDate;
47  import org.orekit.utils.PVCoordinates;
48  import org.orekit.utils.ParameterDriver;
49  import org.orekit.utils.ParameterDriversList;
50  import org.orekit.utils.TimeSpanMap.Span;
51  
52  import java.util.ArrayList;
53  import java.util.Collections;
54  import java.util.List;
55  
56  /** Test class for Kalman model.
57   * This class is deeply entangled with KalmanEstimator class. Thus it is difficult to test as a stand-alone.
58   * Here we simply test the functions from KalmanEstimation interface that return the "physical" values of the
59   * different state vectors and matrices used in a Kalman filter:
60   * state transition, measurement, kalman gain matrices etc.
61   * @author Maxime Journot
62   */
63  public class KalmanModelTest {
64  
65      /** Orbit type for propagation. */
66      private final OrbitType orbitType = OrbitType.CARTESIAN;
67  
68      /** Position angle for propagation. */
69      private final PositionAngleType positionAngleType = PositionAngleType.TRUE;
70  
71      /** Initial orbit. */
72      private Orbit orbit0;
73  
74      /** Propagator builder. */
75      private NumericalPropagatorBuilder propagatorBuilder;
76  
77      /** Covariance matrix provider. */
78      private CovarianceMatrixProvider covMatrixProvider;
79  
80      /** Estimated measurement parameters list. */
81      private ParameterDriversList estimatedMeasurementsParameters;
82  
83      /** Kalman extended estimator containing models. */
84      private KalmanEstimator kalman;
85  
86      /** Kalman observer. */
87      private ModelLogger modelLogger;
88  
89      /** State size. */
90      private int M;
91  
92      /** PV at t0. */
93      private PV pv;
94  
95      /** Range after t0. */
96      private Range range;
97  
98      /** Driver for satellite range bias. */
99      private ParameterDriver satRangeBiasDriver;
100 
101     /** Driver for SRP coefficient. */
102     private ParameterDriver srpCoefDriver;
103 
104     /** Tolerance for the test. */
105     private final double tol = 2e-16;
106 
107     /** Setup an eccentric orbit with Keplerian gravity force and isotropic solar radiation pressure model.
108      * Add a satellite range bias to the range measurements and select its estimation.
109      * Select estimation of radiation pressure reflection coefficient.
110      * Select estimation of all orbital parameters.
111      * Create a Kalman filter from this.
112      *
113      * Create one perfect PV measurement at t0
114      * Create one range measurement at t0 + 10s, modified by the satellite range bias mentionned above.
115      */
116     @BeforeEach
117     public void setup() {
118         // Create context
119         final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
120 
121         // Initial orbit and date
122         this.orbit0 = context.initialOrbit;
123         ObservableSatellite sat = new ObservableSatellite(0);
124 
125         // Create propagator builder
126         this.propagatorBuilder = context.createBuilder(orbitType, positionAngleType, true,
127                                                        1.0e-6, 60.0, 10., Force.SOLAR_RADIATION_PRESSURE);
128 
129         // Create PV at t0
130         final AbsoluteDate date0 = context.initialOrbit.getDate();
131         this.pv = new PV(date0,
132                              context.initialOrbit.getPosition(),
133                              context.initialOrbit.getPVCoordinates().getVelocity(),
134                              new double[] {1., 2., 3., 1e-3, 2e-3, 3e-3}, 1.,
135                              sat);
136 
137         // Create one 0m range measurement at t0 + 10s
138         final AbsoluteDate date  = date0.shiftedBy(10.);
139         final GroundStation station = context.stations.get(0);
140         this.range = new Range(station, true, date, 18616150., 10., 1., sat);
141         // Exact range value is 1.8616150246470984E7 m
142 
143         // Add sat range bias to PV and select it
144         final Bias<Range> satRangeBias = new Bias<Range>(new String[] {"sat range bias"},
145                                                          new double[] {100.},
146                                                          new double[] {10.},
147                                                          new double[] {0.},
148                                                          new double[] {100.});
149         this.satRangeBiasDriver = satRangeBias.getParametersDrivers().get(0);
150         satRangeBiasDriver.setSelected(true);
151         satRangeBiasDriver.setReferenceDate(date);
152         range.addModifier(satRangeBias);
153         for (ParameterDriver driver : range.getParametersDrivers()) {
154             driver.setReferenceDate(date);
155         }
156 
157         // Gather list of meas parameters (only sat range bias here)
158         this.estimatedMeasurementsParameters = new ParameterDriversList();
159         for (final ParameterDriver driver : range.getParametersDrivers()) {
160             if (driver.isSelected()) {
161                 estimatedMeasurementsParameters.add(driver);
162             }
163         }
164         // Select SRP coefficient
165         this.srpCoefDriver = propagatorBuilder.getPropagationParametersDrivers().
166                         findByName(RadiationSensitive.REFLECTION_COEFFICIENT);
167         srpCoefDriver.setReferenceDate(date);
168         srpCoefDriver.setSelected(true);
169 
170         // Create a covariance matrix using the scales of the estimated parameters
171         final double[] scales = getParametersScale(propagatorBuilder, estimatedMeasurementsParameters);
172         this.M = scales.length;
173         this.covMatrixProvider = setInitialCovarianceMatrix(scales);
174 
175         // Initialize Kalman
176         final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
177         kalmanBuilder.addPropagationConfiguration(propagatorBuilder, covMatrixProvider);
178         kalmanBuilder.estimatedMeasurementsParameters(estimatedMeasurementsParameters, null);
179         this.kalman = kalmanBuilder.build();
180         this.modelLogger = new ModelLogger();
181         kalman.setObserver(modelLogger);
182     }
183 
184     /** Test of the physical matrices and vectors returned by the methods from KalmanEstimation interface.
185      *  First, we perform a check before any measurement is added. Most of the matrices should be null.
186      *  Then we process a perfect PV at t0 in the Kalman and check the matrices.
187      *  Finally we process a range measurement after t0 in the Kalman and check the matrices.
188      */
189     @Test
190     public void ModelPhysicalOutputsTest() {
191 
192         // Check model at t0 before any measurement is added
193         // -------------------------------------------------
194         checkModelAtT0();
195 
196 
197         // Check model after PV measurement at t0 is added
198         // -----------------------------------------------
199 
200         // Constant process noise covariance matrix Q
201         final RealMatrix Q = covMatrixProvider.getProcessNoiseMatrix(new SpacecraftState(orbit0),
202                                                                      new SpacecraftState(orbit0));
203 
204         // Initial covariance matrix
205         final RealMatrix P0 = covMatrixProvider.getInitialCovarianceMatrix(new SpacecraftState(orbit0));
206 
207         // Physical predicted covariance matrix at t0
208         // State transition matrix is the identity matrix at t0
209         RealMatrix Ppred = P0.add(Q);
210 
211         // Predicted orbit is equal to initial orbit at t0
212         Orbit orbitPred = orbit0;
213 
214         // Expected measurement matrix for a PV measurement is the 6-sized identity matrix for cartesian orbital parameters
215         // + zeros for other estimated parameters
216         RealMatrix expH = MatrixUtils.createRealMatrix(6, M);
217         for (int i = 0; i < 6; i++) {
218             expH.setEntry(i, i, 1.);
219         }
220 
221         // Expected state transition matrix
222         // State transition matrix is the identity matrix at t0
223         RealMatrix expPhi = MatrixUtils.createRealIdentityMatrix(M);
224         // Add PV measurement and check model afterwards
225         checkModelAfterMeasurementAdded(1, pv, Ppred, orbitPred, expPhi, expH);
226 
227 
228         // Check model after range measurement after t0 is added
229         // -----------------------------------------------------
230 
231         // Get the estimated propagator from Kalman filter and propagate it to
232         // range measurement date
233         final Propagator propagator = propagatorBuilder.buildPropagator();
234 
235         // Set derivatives computation for the propagator
236         final String equationName = KalmanEstimator.class.getName() + "-derivatives-";
237         final MatricesHarvester harvester = propagator.setupMatricesComputation(equationName, null, null);
238 
239         // Propagate to range date and get predicted orbit
240         final SpacecraftState scPred = propagator.propagate(range.getDate());
241         orbitPred = scPred.getOrbit();
242 
243         // Expected state transition matrix
244         expPhi = MatrixUtils.createRealIdentityMatrix(M);
245 
246         // Derivatives of the state vector with respect to initial state vector
247         final double[][] dYdY0 =  harvester.getStateTransitionMatrix(scPred).getData();
248         expPhi.setSubMatrix(dYdY0, 0, 0);
249 
250         // Derivatives of SRP coef with respect to state
251         final double[][] dYdPp  = harvester.getParametersJacobian(scPred).getData();
252         expPhi.setSubMatrix(dYdPp, 0, 6);
253 
254         // Estimated cov matrix from last measurement
255         RealMatrix Pest = kalman.getPhysicalEstimatedCovarianceMatrix();
256 
257         // Predicted covariance matrix for this measurement
258         Ppred = expPhi.multiply(Pest.multiplyTransposed(expPhi)).add(Q);
259 
260         // Expected measurement matrix
261         expH = MatrixUtils.createRealMatrix(1, M);
262         // State part
263         EstimatedMeasurement<Range> rangeEstimated = range.estimate(0, 0, new SpacecraftState[] {scPred});
264         final RealMatrix dMdY = MatrixUtils.createRealMatrix(rangeEstimated.getStateDerivatives(0));
265         expH.setSubMatrix(dMdY.getData(), 0, 0);
266         // SRP part
267         final SpacecraftState scTransition = scPred.shiftedBy(-rangeEstimated.getTimeOffset());
268         final double[][] dYdPpTransition = harvester.getParametersJacobian(scTransition).getData();
269         final RealMatrix dMdCr = dMdY.multiply(MatrixUtils.createRealMatrix(dYdPpTransition));
270         expH.setEntry(0, 6, dMdCr.getEntry(0, 0));
271         // Sat range bias part
272         expH.setEntry(0, 7, rangeEstimated.getParameterDerivatives(satRangeBiasDriver, new AbsoluteDate())[0]);
273         
274         // Add range measurement and check model afterwards
275         checkModelAfterMeasurementAdded(2, range, Ppred, orbitPred, expPhi, expH);
276     }
277 
278     /** Check the model physical outputs at t0 before any measurement is added. */
279     private void checkModelAtT0() {
280 
281         // Instantiate a Model from attributes
282         final KalmanModel model = new KalmanModel(Collections.singletonList(propagatorBuilder),
283                                                   Collections.singletonList(covMatrixProvider),
284                                                   estimatedMeasurementsParameters,
285                                                   null);
286 
287         // Evaluate at t0
288         // --------------
289 
290         // Time
291         Assertions.assertEquals(0., model.getEstimate().getTime(), 0.);
292         Assertions.assertEquals(0., model.getCurrentDate().durationFrom(orbit0.getDate()), 0.);
293 
294         // Measurement number
295         Assertions.assertEquals(0, model.getCurrentMeasurementNumber());
296 
297         // Normalized state - is zeros
298         final RealVector stateN = model.getEstimate().getState();
299         Assertions.assertArrayEquals(new double[M], stateN.toArray(), tol);
300 
301         // Physical state - = initialized
302         final RealVector x = model.getPhysicalEstimatedState();
303         final RealVector expX = MatrixUtils.createRealVector(M);
304         final double[] orbitState0 = new double[6];
305         orbitType.mapOrbitToArray(orbit0, positionAngleType, orbitState0, null);
306         expX.setSubVector(0, MatrixUtils.createRealVector(orbitState0));
307         expX.setEntry(6, srpCoefDriver.getReferenceValue());
308         expX.setEntry(7, satRangeBiasDriver.getReferenceValue());
309         final double[] dX = x.subtract(expX).toArray();
310         Assertions.assertArrayEquals(new double[M], dX, tol);
311 
312         // Normalized covariance - filled with 1
313         final double[][] Pn = model.getEstimate().getCovariance().getData();
314         final double[][] expPn = new double[M][M];
315         for (int i = 0; i < M; i++) {
316             for (int j = 0; j < M; j++) {
317                 expPn[i][j] = 1.;
318             }
319             Assertions.assertArrayEquals(expPn[i], Pn[i], tol, "Failed on line " + i);
320         }
321 
322         // Physical covariance = initialized
323         final RealMatrix P   = model.getPhysicalEstimatedCovarianceMatrix();
324         final RealMatrix expP = covMatrixProvider.getInitialCovarianceMatrix(new SpacecraftState(orbit0));
325         final double[][] dP = P.subtract(expP).getData();
326         for (int i = 0; i < M; i++) {
327             Assertions.assertArrayEquals(new double[M], dP[i], tol, "Failed on line " + i);
328         }
329 
330         // Check that other "physical" matrices are null
331         Assertions.assertNull(model.getEstimate().getInnovationCovariance());
332         Assertions.assertNull(model.getPhysicalInnovationCovarianceMatrix());
333         Assertions.assertNull(model.getEstimate().getKalmanGain());
334         Assertions.assertNull(model.getPhysicalKalmanGain());
335         Assertions.assertNull(model.getEstimate().getMeasurementJacobian());
336         Assertions.assertNull(model.getPhysicalMeasurementJacobian());
337         Assertions.assertNull(model.getEstimate().getStateTransitionMatrix());
338         Assertions.assertNull(model.getPhysicalStateTransitionMatrix());
339     }
340 
341     /** Add a measurement to the Kalman filter.
342      * Check model physical outputs afterwards.
343      */
344     private void checkModelAfterMeasurementAdded(final int expMeasurementNumber,
345                                                 final ObservedMeasurement<?> meas,
346                                                 final RealMatrix expPpred,
347                                                 final Orbit expOrbitPred,
348                                                 final RealMatrix expPhi,
349                                                 final RealMatrix expH) {
350 
351         // Predicted value of SRP coef and sat range bias
352         // (= value before adding measurement to the filter)
353         final double srpCoefPred = srpCoefDriver.getValue();
354         final double satRangeBiasPred = satRangeBiasDriver.getValue();
355         
356         // Expected predicted measurement
357         final double[] expMeasPred =
358                         meas.estimate(0, 0,
359                                       new SpacecraftState[] {new SpacecraftState(expOrbitPred)}).getEstimatedValue();
360 
361         // Process PV measurement in Kalman and get model
362         kalman.processMeasurements(Collections.singletonList(meas));
363         KalmanEstimation model = modelLogger.estimation;
364 
365         // Measurement size
366         final int N = meas.getDimension();
367 
368         // Time
369         Assertions.assertEquals(0., model.getCurrentDate().durationFrom(expOrbitPred.getDate()), 0.);
370 
371         // Measurement number
372         Assertions.assertEquals(expMeasurementNumber, model.getCurrentMeasurementNumber());
373 
374         // State transition matrix
375         final RealMatrix phi    = model.getPhysicalStateTransitionMatrix();
376         final double[][] dPhi   = phi.subtract(expPhi).getData();
377         for (int i = 0; i < M; i++) {
378             Assertions.assertArrayEquals(new double[M], dPhi[i], tol*100, "Failed on line " + i);
379         }
380 
381         // Measurement matrix
382         final RealMatrix H = model.getPhysicalMeasurementJacobian();
383         final double[][] dH = H.subtract(expH).getData();
384         for (int i = 0; i < N; i++) {
385             Assertions.assertArrayEquals(new double[M], dH[i], tol, "Failed on line " + i);
386         }
387 
388         // Measurement covariance matrix
389         final double[] measSigmas = meas.getTheoreticalStandardDeviation();
390         final RealMatrix R = MatrixUtils.createRealMatrix(N, N);
391         for (int i = 0; i < N; i++) {
392             R.setEntry(i, i, measSigmas[i] * measSigmas[i]);
393         }
394 
395         // Innovation matrix
396         final RealMatrix expS = expH.multiply(expPpred.multiplyTransposed(expH)).add(R);
397         final RealMatrix S = model.getPhysicalInnovationCovarianceMatrix();
398         final double[][] dS = S.subtract(expS).getData();
399         for (int i = 0; i < N; i++) {
400             Assertions.assertArrayEquals(new double[N], dS[i], tol*1e4, "Failed on line \" + i");
401         }
402 
403         // Kalman gain
404         final RealMatrix expK = expPpred.multiplyTransposed(expH).multiply(new LUDecomposition(expS).getSolver().getInverse());
405         final RealMatrix K = model.getPhysicalKalmanGain();
406         final double[][] dK = K.subtract(expK).getData();
407         for (int i = 0; i < M; i++) {
408             Assertions.assertArrayEquals(new double[N], dK[i], tol*1e5, "Failed on line " + i);
409         }
410 
411         // Predicted orbit
412         final Orbit orbitPred = model.getPredictedSpacecraftStates()[0].getOrbit();
413         final PVCoordinates pvOrbitPred = orbitPred.getPVCoordinates();
414         final PVCoordinates expPVOrbitPred = expOrbitPred.getPVCoordinates();
415         final double dpOrbitPred = Vector3D.distance(expPVOrbitPred.getPosition(), pvOrbitPred.getPosition());
416         final double dvOrbitPred = Vector3D.distance(expPVOrbitPred.getVelocity(), pvOrbitPred.getVelocity());
417         Assertions.assertEquals(0., dpOrbitPred, tol);
418         Assertions.assertEquals(0., dvOrbitPred, tol);
419 
420         // Predicted measurement
421         final double[] measPred = model.getPredictedMeasurement().getEstimatedValue();
422         Assertions.assertArrayEquals(expMeasPred, measPred, tol);
423 
424         // Predicted state
425         final double[] orbitPredState = new double[6];
426         orbitPred.getType().mapOrbitToArray(orbitPred, PositionAngleType.TRUE, orbitPredState, null);
427         final RealVector expXpred = MatrixUtils.createRealVector(M);
428         for (int i = 0; i < 6; i++) {
429             expXpred.setEntry(i, orbitPredState[i]);
430         }
431         expXpred.setEntry(6, srpCoefPred);
432         expXpred.setEntry(7, satRangeBiasPred);
433 
434         // Innovation vector
435         final RealVector observedMeas  = MatrixUtils.createRealVector(model.getPredictedMeasurement().getObservedValue());
436         final RealVector predictedMeas = MatrixUtils.createRealVector(model.getPredictedMeasurement().getEstimatedValue());
437         final RealVector innovation = observedMeas.subtract(predictedMeas);
438 
439         // Corrected state
440         final RealVector expectedXcor = expXpred.add(expK.operate(innovation));
441         final RealVector Xcor = model.getPhysicalEstimatedState();
442         final double[] dXcor = Xcor.subtract(expectedXcor).toArray();
443         Assertions.assertArrayEquals(new double[M], dXcor, tol);
444 
445         // Corrected covariance
446         final RealMatrix expectedPcor =
447                         (MatrixUtils.createRealIdentityMatrix(M).
448                         subtract(expK.multiply(expH))).multiply(expPpred);
449         final RealMatrix Pcor = model.getPhysicalEstimatedCovarianceMatrix();
450         final double[][] dPcor = Pcor.subtract(expectedPcor).getData();
451         for (int i = 0; i < M; i++) {
452             Assertions.assertArrayEquals(new double[M], dPcor[i], tol*1e5, "Failed on line " + i);
453         }
454     }
455 
456     /** Return expected Kalman current state using propagator builder and measurement driver list.
457      * @param stateSize the size of the state vector
458      * @param builder propagator builder
459      * @param estimatedMeasurementsParameters estimated measurements parameters
460      * @return expected Kalman state
461      */
462     /*private RealVector getExpectedKalmanState(final int stateSize,
463                                              final NumericalPropagatorBuilder builder,
464                                              ParameterDriversList estimatedMeasurementsParameters) {
465 
466         final double[] kalmanState = new double[stateSize];
467 
468         int i = 0;
469         // Orbital parameters
470         for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
471             if (driver.isSelected()) {
472                 kalmanState[i++] = driver.getValue();
473             }
474         }
475 
476         // Propagation parameters
477         for (ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
478             if (driver.isSelected()) {
479                 kalmanState[i++] = driver.getValue();
480             }
481         }
482 
483         // Measurement parameters
484         for (ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
485             if (driver.isSelected()) {
486                 kalmanState[i++] = driver.getValue();
487             }
488         }
489 
490         return MatrixUtils.createRealVector(kalmanState);
491     }*/
492 
493     /** Get an array of the scales of the estimated parameters.
494      * @param builder propagator builder
495      * @param estimatedMeasurementsParameters estimated measurements parameters
496      * @return array containing the scales of the estimated parameter
497      */
498     private double[] getParametersScale(final NumericalPropagatorBuilder builder,
499                                        ParameterDriversList estimatedMeasurementsParameters) {
500         final List<Double> scaleList = new ArrayList<>();
501 
502         // Orbital parameters
503         for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
504             if (driver.isSelected()) {
505             	for (Span<Double> span = driver.getValueSpanMap().getFirstSpan(); span != null; span = span.next()) {
506                     scaleList.add(driver.getScale());
507             	}
508             }
509         }
510 
511         // Propagation parameters
512         for (ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
513             if (driver.isSelected()) {
514             	for (Span<Double> span = driver.getValueSpanMap().getFirstSpan(); span != null; span = span.next()) {
515                     scaleList.add(driver.getScale());
516             	}
517             }
518         }
519 
520         // Measurement parameters
521         for (ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
522             if (driver.isSelected()) {
523             	for (Span<Double> span = driver.getValueSpanMap().getFirstSpan(); span != null; span = span.next()) {
524                     scaleList.add(driver.getScale());
525             	}
526             }
527         }
528 
529         final double[] scales = new double[scaleList.size()];
530         for (int i = 0; i < scaleList.size(); i++) {
531             scales[i] = scaleList.get(i);
532         }
533         return scales;
534     }
535 
536     /** Create a covariance matrix provider with initial and process noise matrix constant and identical.
537      * Each element Pij of the initial covariance matrix equals:
538      * Pij = scales[i]*scales[j]
539      * @param scales scales of the estimated parameters
540      * @return covariance matrix provider
541      */
542     private CovarianceMatrixProvider setInitialCovarianceMatrix(final double[] scales) {
543 
544         final int n = scales.length;
545         final RealMatrix cov = MatrixUtils.createRealMatrix(n, n);
546         for (int i = 0; i < n; i++) {
547             for (int j = 0; j < n; j++) {
548                 cov.setEntry(i, j, scales[i] * scales[j]);
549             }
550         }
551         return new ConstantProcessNoise(cov);
552     }
553 
554 
555     /** Observer allowing to get Kalman model after a measurement was processed in the Kalman filter. */
556     public class ModelLogger implements KalmanObserver {
557         KalmanEstimation estimation;
558 
559         @Override
560         public void evaluationPerformed(KalmanEstimation estimation) {
561             this.estimation = estimation;
562         }
563     }
564 
565     /** 
566      * Test that mass depletions during the processing of a measurement are propagated back to the PropagatorBuilder
567      * so that they are taken into account during processing of subsequent measurements.
568      */
569     @Test
570     public void MassDepletionTest()  {
571         // Add a maneuver with nonzero mass expenditure between the first and second measurement dates.
572         AbsoluteDate initialDate = this.propagatorBuilder.getInitialOrbitDate();
573         double initialMass = this.propagatorBuilder.getMass();
574         AbsoluteDate maneuverDate = initialDate.shiftedBy(5.0);
575 
576         ConstantThrustManeuver maneuver = new ConstantThrustManeuver(maneuverDate, 10.0,  10.0, 100.0, new Vector3D(1.0, 0.0, 0.0));
577         this.propagatorBuilder.addForceModel(maneuver);
578     
579         kalman.processMeasurements(Collections.singletonList(pv));
580         kalman.processMeasurements(Collections.singletonList(range));
581         
582         double postManeuverPredictedStateMass = modelLogger.estimation.getPredictedSpacecraftStates()[0].getMass();
583         double postManeuverCorrectedStateMass = modelLogger.estimation.getCorrectedSpacecraftStates()[0].getMass();
584         double postManeuverPropagatorMass = propagatorBuilder.getMass();
585 
586         // The mass of the predicted state should be less than the initial mass, because mass should have been expended during the maneuver.
587         Assertions.assertTrue(postManeuverPredictedStateMass < initialMass);
588 
589         // The mass for the propagator builder should have been updated to be equal to the mass from the 
590         // newly predicted state, which includes the maneuver.
591         Assertions.assertEquals(postManeuverPredictedStateMass, postManeuverPropagatorMass); 
592 
593         // The mass of the corrected state should be equal to the mass from the propagator builder,
594         Assertions.assertEquals(postManeuverCorrectedStateMass, postManeuverPropagatorMass);    
595     }
596 }