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