1   /* Copyright 2002-2024 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.ArrayList;
47  import java.util.List;
48  
49  class KeplerianBatchLSEstimatorTest {
50  
51      /**
52       * Perfect PV measurements with a perfect start
53       */
54      @Test
55      void testPV() {
56  
57          KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
58  
59          final KeplerianPropagatorBuilder propagatorBuilder = context.createBuilder(PositionAngleType.MEAN, true, 1.0);
60  
61          // create perfect PV measurements
62          final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
63                                                                                           propagatorBuilder);
64          final List<ObservedMeasurement<?>> measurements =
65                          KeplerianEstimationTestUtils.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          KeplerianEstimationTestUtils.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      void testKeplerPVBackward() {
97  
98          KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
99  
100         final KeplerianPropagatorBuilder propagatorBuilder =
101                         context.createBuilder(PositionAngleType.MEAN, true, 1.0);
102 
103         // create perfect PV measurements
104         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
105                                                                                          propagatorBuilder);
106         final List<ObservedMeasurement<?>> measurements =
107                         KeplerianEstimationTestUtils.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         KeplerianEstimationTestUtils.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         Assertions.assertEquals(0.0,     physicalCovariances.getEntry(0, 0), 1.7e-15);
134 
135     }
136 
137     /**
138      * Perfect range measurements with a perfect start
139      */
140     @Test
141     void testKeplerRange() {
142 
143         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
144 
145         final KeplerianPropagatorBuilder propagatorBuilder =
146                         context.createBuilder(PositionAngleType.MEAN, true, 1.0);
147 
148         // create perfect range measurements
149         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
150                                                                            propagatorBuilder);
151         final List<ObservedMeasurement<?>> measurements =
152                         KeplerianEstimationTestUtils.createMeasurements(propagator,
153                                                                new TwoWayRangeMeasurementCreator(context),
154                                                                1.0, 3.0, 300.0);
155 
156         // create orbit estimator
157         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
158                                                                 propagatorBuilder);
159         for (final ObservedMeasurement<?> range : measurements) {
160             estimator.addMeasurement(range);
161         }
162         estimator.setParametersConvergenceThreshold(1.0e-2);
163         estimator.setMaxIterations(10);
164         estimator.setMaxEvaluations(20);
165 
166         KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 4,
167                                                    0.0, 8.5e-7,
168                                                    0.0, 2.0e-6,
169                                                    0.0, 1.9e-8,
170                                                    0.0, 9.0e-12);
171 
172     }
173 
174     /**
175      * Perfect range measurements with a perfect start and an on-board antenna range offset
176      */
177     @Test
178     void testKeplerRangeWithOnBoardAntennaOffset() {
179 
180         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
181 
182         final KeplerianPropagatorBuilder propagatorBuilder =
183                         context.createBuilder(PositionAngleType.MEAN, true, 1.0);
184         propagatorBuilder.setAttitudeProvider(new LofOffset(propagatorBuilder.getFrame(), LOFType.LVLH));
185         final Vector3D antennaPhaseCenter = new Vector3D(-1.2, 2.3, -0.7);
186 
187         // create perfect range measurements with antenna offset
188         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
189                                                                            propagatorBuilder);
190         final List<ObservedMeasurement<?>> measurements =
191                         KeplerianEstimationTestUtils.createMeasurements(propagator,
192                                                                new TwoWayRangeMeasurementCreator(context,
193                                                                                                  Vector3D.ZERO, null,
194                                                                                                  antennaPhaseCenter, null,
195                                                                                                  0),
196                                                                1.0, 3.0, 300.0);
197 
198         // create orbit estimator
199         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
200                                                                 propagatorBuilder);
201         final PhaseCentersRangeModifier obaModifier = new PhaseCentersRangeModifier(FrequencyPattern.ZERO_CORRECTION,
202                                                                                     new FrequencyPattern(antennaPhaseCenter,
203                                                                                                          null));
204         for (final ObservedMeasurement<?> range : measurements) {
205             ((Range) range).addModifier(obaModifier);
206             estimator.addMeasurement(range);
207         }
208         estimator.setParametersConvergenceThreshold(1.0e-2);
209         estimator.setMaxIterations(10);
210         estimator.setMaxEvaluations(20);
211 
212         KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 12,
213                                                    0.0, 5.9e-5,
214                                                    0.0, 1.5e-4,
215                                                    0.0, 4.3e-9,
216                                                    0.0, 4.2e-12);
217 
218     }
219 
220     /**
221      * Perfect range rate measurements with a perfect start
222      */
223     @Test
224     void testKeplerRangeRate() {
225 
226         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
227 
228         final KeplerianPropagatorBuilder propagatorBuilder =
229                         context.createBuilder(PositionAngleType.MEAN, true, 1.0);
230 
231         // create perfect range rate measurements
232         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
233                                                                                          propagatorBuilder);
234         final double groundClockDrift =  4.8e-9;
235         for (final GroundStation station : context.stations) {
236             station.getClockDriftDriver().setValue(groundClockDrift);
237         }
238         final double satClkDrift = 3.2e-10;
239         final List<ObservedMeasurement<?>> measurements1 =
240                         KeplerianEstimationTestUtils.createMeasurements(propagator,
241                                                                new RangeRateMeasurementCreator(context, false, satClkDrift),
242                                                                1.0, 3.0, 300.0);
243 
244         final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
245         measurements.addAll(measurements1);
246 
247         // create orbit estimator
248         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
249                                                                 propagatorBuilder);
250         for (final ObservedMeasurement<?> rangerate : measurements) {
251             estimator.addMeasurement(rangerate);
252         }
253         estimator.setParametersConvergenceThreshold(1.0e-3);
254         estimator.setMaxIterations(10);
255         estimator.setMaxEvaluations(20);
256 
257         KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 3,
258                                                    0.0, 5.6e-7,
259                                                    0.0, 7.7e-7,
260                                                    0.0, 8.7e-5,
261                                                    0.0, 3.5e-8);
262     }
263 
264     @Test
265     void testWrappedException() {
266 
267         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
268 
269         final KeplerianPropagatorBuilder propagatorBuilder =
270                         context.createBuilder(PositionAngleType.MEAN, true, 1.0);
271 
272         // create perfect range measurements
273         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
274                                                                                          propagatorBuilder);
275         final List<ObservedMeasurement<?>> measurements =
276                         KeplerianEstimationTestUtils.createMeasurements(propagator,
277                                                                              new TwoWayRangeMeasurementCreator(context),
278                                                                              1.0, 3.0, 300.0);
279 
280         // create orbit estimator
281         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
282                                                                 propagatorBuilder);
283         for (final ObservedMeasurement<?> range : measurements) {
284             estimator.addMeasurement(range);
285         }
286         estimator.setParametersConvergenceThreshold(1.0e-2);
287         estimator.setMaxIterations(10);
288         estimator.setMaxEvaluations(20);
289         estimator.setObserver(new BatchLSObserver() {
290             /** {@inheritDoc} */
291             @Override
292             public void evaluationPerformed(int iterationsCount, int evaluationscount,
293                                            Orbit[] orbits,
294                                            ParameterDriversList estimatedOrbitalParameters,
295                                            ParameterDriversList estimatedPropagatorParameters,
296                                            ParameterDriversList estimatedMeasurementsParameters,
297                                            EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
298                 throw new DummyException();
299             }
300         });
301 
302         try {
303             KeplerianEstimationTestUtils.checkFit(context, estimator, 3, 4,
304                                                        0.0, 1.5e-6,
305                                                        0.0, 3.2e-6,
306                                                        0.0, 3.8e-7,
307                                                        0.0, 1.5e-10);
308             Assertions.fail("an exception should have been thrown");
309         } catch (DummyException de) {
310             // expected
311         }
312 
313     }
314 
315     private static class DummyException extends OrekitException {
316         private static final long serialVersionUID = 1L;
317         public DummyException() {
318             super(OrekitMessages.INTERNAL_ERROR);
319         }
320     }
321 
322 }