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.leastsquares;
18  
19  import org.hipparchus.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.linear.RealMatrix;
21  import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation;
22  import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
23  import org.junit.jupiter.api.Assertions;
24  import org.junit.jupiter.api.Test;
25  import org.orekit.attitudes.LofOffset;
26  import org.orekit.errors.OrekitException;
27  import org.orekit.errors.OrekitMessages;
28  import org.orekit.estimation.KeplerianContext;
29  import org.orekit.estimation.KeplerianEstimationTestUtils;
30  import org.orekit.estimation.measurements.EstimationsProvider;
31  import org.orekit.estimation.measurements.GroundStation;
32  import org.orekit.estimation.measurements.ObservedMeasurement;
33  import org.orekit.estimation.measurements.PVMeasurementCreator;
34  import org.orekit.estimation.measurements.Range;
35  import org.orekit.estimation.measurements.TwoWayRangeMeasurementCreator;
36  import org.orekit.estimation.measurements.RangeRateMeasurementCreator;
37  import org.orekit.estimation.measurements.modifiers.PhaseCentersRangeModifier;
38  import org.orekit.frames.LOFType;
39  import org.orekit.gnss.antenna.FrequencyPattern;
40  import org.orekit.orbits.Orbit;
41  import org.orekit.orbits.PositionAngleType;
42  import org.orekit.propagation.Propagator;
43  import org.orekit.propagation.conversion.KeplerianPropagatorBuilder;
44  import org.orekit.utils.ParameterDriversList;
45  
46  import java.util.List;
47  
48  class KeplerianBatchLSEstimatorTest {
49  
50      /**
51       * Perfect PV measurements with a perfect start
52       */
53      @Test
54      void testPV() {
55  
56          KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
57  
58          final KeplerianPropagatorBuilder propagatorBuilder = context.createBuilder(PositionAngleType.MEAN, true, 1.0);
59  
60          // create perfect PV measurements
61          final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
62                                                                                           propagatorBuilder);
63          final List<ObservedMeasurement<?>> measurements =
64                          KeplerianEstimationTestUtils.createMeasurements(propagator,
65                                                                               new PVMeasurementCreator(),
66                                                                               0.0, 1.0, 300.0);
67  
68          // create orbit estimator
69          final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
70                                                                  propagatorBuilder);
71          for (final ObservedMeasurement<?> measurement : measurements) {
72              estimator.addMeasurement(measurement);
73          }
74          estimator.setParametersConvergenceThreshold(1.0e-2);
75          estimator.setMaxIterations(10);
76          estimator.setMaxEvaluations(20);
77  
78          KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 1,
79                                                     0.0, 1.0e-15,
80                                                     0.0, 1.0e-15,
81                                                     0.0, 1.0e-15,
82                                                     0.0, 1.0e-15);
83  
84          RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
85          RealMatrix physicalCovariances   = estimator.getPhysicalCovariances(1.0e-10);
86          Assertions.assertEquals(6,       normalizedCovariances.getRowDimension());
87          Assertions.assertEquals(6,       normalizedCovariances.getColumnDimension());
88          Assertions.assertEquals(6,       physicalCovariances.getRowDimension());
89          Assertions.assertEquals(6,       physicalCovariances.getColumnDimension());
90  
91      }
92  
93      /** Test PV measurements generation and backward propagation in least-square orbit determination. */
94      @Test
95      void testKeplerPVBackward() {
96  
97          KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
98  
99          final KeplerianPropagatorBuilder propagatorBuilder =
100                         context.createBuilder(PositionAngleType.MEAN, true, 1.0);
101 
102         // create perfect PV measurements
103         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
104                                                                                          propagatorBuilder);
105         final List<ObservedMeasurement<?>> measurements =
106                         KeplerianEstimationTestUtils.createMeasurements(propagator,
107                                                                new PVMeasurementCreator(),
108                                                                0.0, -1.0, 300.0);
109 
110         // create orbit estimator
111         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
112                                                                 propagatorBuilder);
113         for (final ObservedMeasurement<?> measurement : measurements) {
114             estimator.addMeasurement(measurement);
115         }
116         estimator.setParametersConvergenceThreshold(1.0e-2);
117         estimator.setMaxIterations(10);
118         estimator.setMaxEvaluations(20);
119 
120         KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 1,
121                                                    0.0, 1.0e-15,
122                                                    0.0, 1.0e-15,
123                                                    0.0, 1.0e-15,
124                                                    0.0, 1.0e-15);
125 
126         RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
127         RealMatrix physicalCovariances   = estimator.getPhysicalCovariances(1.0e-10);
128         Assertions.assertEquals(6,       normalizedCovariances.getRowDimension());
129         Assertions.assertEquals(6,       normalizedCovariances.getColumnDimension());
130         Assertions.assertEquals(6,       physicalCovariances.getRowDimension());
131         Assertions.assertEquals(6,       physicalCovariances.getColumnDimension());
132         Assertions.assertEquals(0.0,     physicalCovariances.getEntry(0, 0), 1.7e-15);
133 
134     }
135 
136     /**
137      * Perfect range measurements with a perfect start
138      */
139     @Test
140     void testKeplerRange() {
141 
142         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
143 
144         final KeplerianPropagatorBuilder propagatorBuilder =
145                         context.createBuilder(PositionAngleType.MEAN, true, 1.0);
146 
147         // create perfect range measurements
148         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
149                                                                            propagatorBuilder);
150         final List<ObservedMeasurement<?>> measurements =
151                         KeplerianEstimationTestUtils.createMeasurements(propagator,
152                                                                new TwoWayRangeMeasurementCreator(context),
153                                                                1.0, 3.0, 300.0);
154 
155         // create orbit estimator
156         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
157                                                                 propagatorBuilder);
158         for (final ObservedMeasurement<?> range : measurements) {
159             estimator.addMeasurement(range);
160         }
161         estimator.setParametersConvergenceThreshold(1.0e-2);
162         estimator.setMaxIterations(10);
163         estimator.setMaxEvaluations(20);
164 
165         KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 2,
166                                                    0.0, 8.5e-7,
167                                                    0.0, 2.6e-6,
168                                                    0.0, 7e-8,
169                                                    0.0, 4.0e-11);
170 
171     }
172 
173     /**
174      * Perfect range measurements with a perfect start and an on-board antenna range offset
175      */
176     @Test
177     void testKeplerRangeWithOnBoardAntennaOffset() {
178 
179         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
180 
181         final KeplerianPropagatorBuilder propagatorBuilder =
182                         context.createBuilder(PositionAngleType.MEAN, true, 1.0);
183         propagatorBuilder.setAttitudeProvider(new LofOffset(propagatorBuilder.getFrame(), LOFType.LVLH));
184         final Vector3D antennaPhaseCenter = new Vector3D(-1.2, 2.3, -0.7);
185 
186         // create perfect range measurements with antenna offset
187         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
188                                                                            propagatorBuilder);
189         final List<ObservedMeasurement<?>> measurements =
190                         KeplerianEstimationTestUtils.createMeasurements(propagator,
191                                                                new TwoWayRangeMeasurementCreator(context,
192                                                                                                  Vector3D.ZERO, null,
193                                                                                                  antennaPhaseCenter, null,
194                                                                                                  0),
195                                                                1.0, 3.0, 300.0);
196 
197         // create orbit estimator
198         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
199                                                                 propagatorBuilder);
200         final PhaseCentersRangeModifier obaModifier = new PhaseCentersRangeModifier(FrequencyPattern.ZERO_CORRECTION,
201                                                                                     new FrequencyPattern(antennaPhaseCenter,
202                                                                                                          null));
203         for (final ObservedMeasurement<?> range : measurements) {
204             ((Range) range).addModifier(obaModifier);
205             estimator.addMeasurement(range);
206         }
207         estimator.setParametersConvergenceThreshold(1.0e-2);
208         estimator.setMaxIterations(10);
209         estimator.setMaxEvaluations(20);
210 
211         KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 13,
212                                                    0.0, 5.9e-5,
213                                                    0.0, 1.5e-4,
214                                                    0.0, 4.3e-9,
215                                                    0.0, 4.2e-12);
216 
217     }
218 
219     /**
220      * Perfect range rate measurements with a perfect start
221      */
222     @Test
223     void testKeplerRangeRate() {
224 
225         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
226 
227         final KeplerianPropagatorBuilder propagatorBuilder =
228                         context.createBuilder(PositionAngleType.MEAN, true, 1.0);
229 
230         // create perfect range rate measurements
231         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
232                                                                                          propagatorBuilder);
233         final double groundClockDrift =  4.8e-9;
234         for (final GroundStation station : context.stations) {
235             station.getClockDriftDriver().setValue(groundClockDrift);
236         }
237         final double satClkDrift = 3.2e-10;
238         final List<ObservedMeasurement<?>> measurements =
239                         KeplerianEstimationTestUtils.createMeasurements(propagator,
240                                                                new RangeRateMeasurementCreator(context, false, satClkDrift),
241                                                                1.0, 3.0, 300.0);
242 
243         // create orbit estimator
244         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
245                                                                 propagatorBuilder);
246         for (final ObservedMeasurement<?> rangerate : measurements) {
247             estimator.addMeasurement(rangerate);
248         }
249         estimator.setParametersConvergenceThreshold(1.0e-3);
250         estimator.setMaxIterations(10);
251         estimator.setMaxEvaluations(20);
252 
253         KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 3,
254                                                    0.0, 5.6e-7,
255                                                    0.0, 7.7e-7,
256                                                    0.0, 8.7e-5,
257                                                    0.0, 3.5e-8);
258     }
259 
260     @Test
261     void testWrappedException() {
262 
263         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
264 
265         final KeplerianPropagatorBuilder propagatorBuilder =
266                         context.createBuilder(PositionAngleType.MEAN, true, 1.0);
267 
268         // create perfect range measurements
269         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
270                                                                                          propagatorBuilder);
271         final List<ObservedMeasurement<?>> measurements =
272                         KeplerianEstimationTestUtils.createMeasurements(propagator,
273                                                                              new TwoWayRangeMeasurementCreator(context),
274                                                                              1.0, 3.0, 300.0);
275 
276         // create orbit estimator
277         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
278                                                                 propagatorBuilder);
279         for (final ObservedMeasurement<?> range : measurements) {
280             estimator.addMeasurement(range);
281         }
282         estimator.setParametersConvergenceThreshold(1.0e-2);
283         estimator.setMaxIterations(10);
284         estimator.setMaxEvaluations(20);
285         estimator.setObserver(new BatchLSObserver() {
286             /** {@inheritDoc} */
287             @Override
288             public void evaluationPerformed(int iterationsCount, int evaluationscount,
289                                            Orbit[] orbits,
290                                            ParameterDriversList estimatedOrbitalParameters,
291                                            ParameterDriversList estimatedPropagatorParameters,
292                                            ParameterDriversList estimatedMeasurementsParameters,
293                                            EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
294                 throw new DummyException();
295             }
296         });
297 
298         try {
299             KeplerianEstimationTestUtils.checkFit(context, estimator, 3, 4,
300                                                        0.0, 1.5e-6,
301                                                        0.0, 3.2e-6,
302                                                        0.0, 3.8e-7,
303                                                        0.0, 1.5e-10);
304             Assertions.fail("an exception should have been thrown");
305         } catch (DummyException de) {
306             // expected
307         }
308 
309     }
310 
311     private static class DummyException extends OrekitException {
312         private static final long serialVersionUID = 1L;
313         public DummyException() {
314             super(OrekitMessages.INTERNAL_ERROR);
315         }
316     }
317 
318 }