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