1   /* Copyright 2002-2022 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.Assert;
27  import org.junit.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.RangeMeasurementCreator;
39  import org.orekit.estimation.measurements.RangeRateMeasurementCreator;
40  import org.orekit.estimation.measurements.modifiers.OnBoardAntennaRangeModifier;
41  import org.orekit.frames.LOFType;
42  import org.orekit.orbits.Orbit;
43  import org.orekit.orbits.PositionAngle;
44  import org.orekit.propagation.Propagator;
45  import org.orekit.propagation.conversion.EcksteinHechlerPropagatorBuilder;
46  import org.orekit.utils.ParameterDriversList;
47  
48  public class EcksteinHechlerBatchLSEstimatorTest {
49  
50      /**
51       * Perfect PV measurements with a perfect start
52       */
53      @Test
54      public void testPV() {
55  
56          EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
57  
58          final EcksteinHechlerPropagatorBuilder propagatorBuilder = context.createBuilder(PositionAngle.MEAN, true, 1.0);
59  
60          // create perfect PV measurements
61          final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
62                                                                                            propagatorBuilder);
63          final List<ObservedMeasurement<?>> measurements =
64                          EcksteinHechlerEstimationTestUtils.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          EcksteinHechlerEstimationTestUtils.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          Assert.assertEquals(6,       normalizedCovariances.getRowDimension());
87          Assert.assertEquals(6,       normalizedCovariances.getColumnDimension());
88          Assert.assertEquals(6,       physicalCovariances.getRowDimension());
89          Assert.assertEquals(6,       physicalCovariances.getColumnDimension());
90  
91      }
92  
93      /** Test PV measurements generation and backward propagation in least-square orbit determination. */
94      @Test
95      public void testKeplerPVBackward() {
96  
97          EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
98  
99          final EcksteinHechlerPropagatorBuilder propagatorBuilder =
100                         context.createBuilder(PositionAngle.MEAN, true, 1.0);
101 
102         // create perfect PV measurements
103         final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
104                                                                                          propagatorBuilder);
105         final List<ObservedMeasurement<?>> measurements =
106                         EcksteinHechlerEstimationTestUtils.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         EcksteinHechlerEstimationTestUtils.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         Assert.assertEquals(6,       normalizedCovariances.getRowDimension());
129         Assert.assertEquals(6,       normalizedCovariances.getColumnDimension());
130         Assert.assertEquals(6,       physicalCovariances.getRowDimension());
131         Assert.assertEquals(6,       physicalCovariances.getColumnDimension());
132 
133     }
134 
135     /**
136      * Perfect range measurements with a perfect start
137      */
138     @Test
139     public void testKeplerRange() {
140 
141         EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
142 
143         final EcksteinHechlerPropagatorBuilder propagatorBuilder =
144                         context.createBuilder(PositionAngle.MEAN, true, 1.0);
145 
146         // create perfect range measurements
147         final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
148                                                                            propagatorBuilder);
149         final List<ObservedMeasurement<?>> measurements =
150                         EcksteinHechlerEstimationTestUtils.createMeasurements(propagator,
151                                                                new RangeMeasurementCreator(context),
152                                                                1.0, 3.0, 300.0);
153 
154         // create orbit estimator
155         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
156                                                                 propagatorBuilder);
157         for (final ObservedMeasurement<?> range : measurements) {
158             estimator.addMeasurement(range);
159         }
160         estimator.setParametersConvergenceThreshold(1.0e-2);
161         estimator.setMaxIterations(10);
162         estimator.setMaxEvaluations(20);
163 
164         EcksteinHechlerEstimationTestUtils.checkFit(context, estimator, 1, 10,
165                                                     0.0, 1.1e-6,
166                                                     0.0, 3.0e-6,
167                                                     0.0, 5.6e-8,
168                                                     0.0, 5.4e-11);
169 
170     }
171 
172     /**
173      * Perfect range measurements with a perfect start and an on-board antenna range offset 
174      */
175     @Test
176     public void testKeplerRangeWithOnBoardAntennaOffset() {
177 
178         EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
179 
180         final EcksteinHechlerPropagatorBuilder propagatorBuilder =
181                         context.createBuilder(PositionAngle.MEAN, true, 1.0);
182         propagatorBuilder.setAttitudeProvider(new LofOffset(propagatorBuilder.getFrame(), LOFType.LVLH));
183         final Vector3D antennaPhaseCenter = new Vector3D(-1.2, 2.3, -0.7);
184 
185         // create perfect range measurements with antenna offset
186         final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
187                                                                            propagatorBuilder);
188         final List<ObservedMeasurement<?>> measurements =
189                         EcksteinHechlerEstimationTestUtils.createMeasurements(propagator,
190                                                                new RangeMeasurementCreator(context, antennaPhaseCenter),
191                                                                1.0, 3.0, 300.0);
192 
193         // create orbit estimator
194         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
195                                                                 propagatorBuilder);
196         final OnBoardAntennaRangeModifier obaModifier = new OnBoardAntennaRangeModifier(antennaPhaseCenter);
197         for (final ObservedMeasurement<?> range : measurements) {
198             ((Range) range).addModifier(obaModifier);
199             estimator.addMeasurement(range);
200         }
201         estimator.setParametersConvergenceThreshold(1.0e-2);
202         estimator.setMaxIterations(10);
203         estimator.setMaxEvaluations(20);
204 
205         EcksteinHechlerEstimationTestUtils.checkFit(context, estimator, 1, 11,
206                                                    0.0, 4.3e-5,
207                                                    0.0, 1.2e-4,
208                                                    0.0, 3.0e-8,
209                                                    0.0, 2.5e-11);
210 
211     }
212 
213     /**
214      * Perfect range rate measurements with a perfect start
215      */
216     @Test
217     public void testKeplerRangeRate() {
218 
219         EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
220 
221         final EcksteinHechlerPropagatorBuilder propagatorBuilder =
222                         context.createBuilder(PositionAngle.MEAN, true, 1.0);
223 
224         // create perfect range rate measurements
225         final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
226                                                                                          propagatorBuilder);
227         final double groundClockDrift =  4.8e-9;
228         for (final GroundStation station : context.stations) {
229             station.getClockDriftDriver().setValue(groundClockDrift);
230         }
231         final double satClkDrift = 3.2e-10;
232         final List<ObservedMeasurement<?>> measurements1 =
233                         EcksteinHechlerEstimationTestUtils.createMeasurements(propagator,
234                                                                new RangeRateMeasurementCreator(context, false, satClkDrift),
235                                                                1.0, 3.0, 300.0);
236 
237         final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
238         measurements.addAll(measurements1);
239 
240         // create orbit estimator
241         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
242                                                                 propagatorBuilder);
243         for (final ObservedMeasurement<?> rangerate : measurements) {
244             estimator.addMeasurement(rangerate);
245         }
246         estimator.setParametersConvergenceThreshold(1.0e-3);
247         estimator.setMaxIterations(10);
248         estimator.setMaxEvaluations(20);
249 
250         EcksteinHechlerEstimationTestUtils.checkFit(context, estimator, 1, 10,
251                                                    0.0, 7.4e-8,
252                                                    0.0, 1.3e-7,
253                                                    0.0, 3.2e-5,
254                                                    0.0, 3.2e-8);
255     }
256 
257     @Test
258     public void testWrappedException() {
259 
260         EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
261 
262         final EcksteinHechlerPropagatorBuilder propagatorBuilder =
263                         context.createBuilder(PositionAngle.MEAN, true, 1.0);
264 
265         // create perfect range measurements
266         final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
267                                                                                          propagatorBuilder);
268         final List<ObservedMeasurement<?>> measurements =
269                         EcksteinHechlerEstimationTestUtils.createMeasurements(propagator,
270                                                                              new RangeMeasurementCreator(context),
271                                                                              1.0, 3.0, 300.0);
272 
273         // create orbit estimator
274         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
275                                                                 propagatorBuilder);
276         for (final ObservedMeasurement<?> range : measurements) {
277             estimator.addMeasurement(range);
278         }
279         estimator.setParametersConvergenceThreshold(1.0e-2);
280         estimator.setMaxIterations(10);
281         estimator.setMaxEvaluations(20);
282         estimator.setObserver(new BatchLSObserver() {
283             /** {@inheritDoc} */
284             @Override
285             public void evaluationPerformed(int iterationsCount, int evaluationscount,
286                                            Orbit[] orbits,
287                                            ParameterDriversList estimatedOrbitalParameters,
288                                            ParameterDriversList estimatedPropagatorParameters,
289                                            ParameterDriversList estimatedMeasurementsParameters,
290                                            EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
291                 throw new DummyException();
292             }
293         });
294 
295         try {
296             EcksteinHechlerEstimationTestUtils.checkFit(context, estimator, 3, 4,
297                                                        0.0, 1.5e-6,
298                                                        0.0, 3.2e-6,
299                                                        0.0, 3.8e-7,
300                                                        0.0, 1.5e-10);
301             Assert.fail("an exception should have been thrown");
302         } catch (DummyException de) {
303             // expected
304         }
305 
306     }
307 
308     private static class DummyException extends OrekitException {
309         private static final long serialVersionUID = 1L;
310         public DummyException() {
311             super(OrekitMessages.INTERNAL_ERROR);
312         }
313     }
314 
315 }