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.KeplerianContext;
32  import org.orekit.estimation.KeplerianEstimationTestUtils;
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.KeplerianPropagatorBuilder;
46  import org.orekit.utils.ParameterDriversList;
47  
48  public class KeplerianBatchLSEstimatorTest {
49  
50      /**
51       * Perfect PV measurements with a perfect start
52       */
53      @Test
54      public void testPV() {
55  
56          KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
57  
58          final KeplerianPropagatorBuilder propagatorBuilder = context.createBuilder(PositionAngle.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          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          KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
98  
99          final KeplerianPropagatorBuilder propagatorBuilder =
100                         context.createBuilder(PositionAngle.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         Assert.assertEquals(6,       normalizedCovariances.getRowDimension());
129         Assert.assertEquals(6,       normalizedCovariances.getColumnDimension());
130         Assert.assertEquals(6,       physicalCovariances.getRowDimension());
131         Assert.assertEquals(6,       physicalCovariances.getColumnDimension());
132         Assert.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     public void testKeplerRange() {
141 
142         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
143 
144         final KeplerianPropagatorBuilder propagatorBuilder =
145                         context.createBuilder(PositionAngle.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 RangeMeasurementCreator(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, 5,
166                                                    0.0, 8.4e-7,
167                                                    0.0, 2.0e-6,
168                                                    0.0, 7.7e-9,
169                                                    0.0, 4.9e-12);
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         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
180 
181         final KeplerianPropagatorBuilder propagatorBuilder =
182                         context.createBuilder(PositionAngle.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 RangeMeasurementCreator(context, antennaPhaseCenter),
192                                                                1.0, 3.0, 300.0);
193 
194         // create orbit estimator
195         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
196                                                                 propagatorBuilder);
197         final OnBoardAntennaRangeModifier obaModifier = new OnBoardAntennaRangeModifier(antennaPhaseCenter);
198         for (final ObservedMeasurement<?> range : measurements) {
199             ((Range) range).addModifier(obaModifier);
200             estimator.addMeasurement(range);
201         }
202         estimator.setParametersConvergenceThreshold(1.0e-2);
203         estimator.setMaxIterations(10);
204         estimator.setMaxEvaluations(20);
205 
206         KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 11,
207                                                    0.0, 5.9e-5,
208                                                    0.0, 1.5e-4,
209                                                    0.0, 7.2e-9,
210                                                    0.0, 7.8e-12);
211 
212     }
213 
214     /**
215      * Perfect range rate measurements with a perfect start
216      */
217     @Test
218     public void testKeplerRangeRate() {
219 
220         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
221 
222         final KeplerianPropagatorBuilder propagatorBuilder =
223                         context.createBuilder(PositionAngle.MEAN, true, 1.0);
224 
225         // create perfect range rate measurements
226         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
227                                                                                          propagatorBuilder);
228         final double groundClockDrift =  4.8e-9;
229         for (final GroundStation station : context.stations) {
230             station.getClockDriftDriver().setValue(groundClockDrift);
231         }
232         final double satClkDrift = 3.2e-10;
233         final List<ObservedMeasurement<?>> measurements1 =
234                         KeplerianEstimationTestUtils.createMeasurements(propagator,
235                                                                new RangeRateMeasurementCreator(context, false, satClkDrift),
236                                                                1.0, 3.0, 300.0);
237 
238         final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
239         measurements.addAll(measurements1);
240 
241         // create orbit estimator
242         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
243                                                                 propagatorBuilder);
244         for (final ObservedMeasurement<?> rangerate : measurements) {
245             estimator.addMeasurement(rangerate);
246         }
247         estimator.setParametersConvergenceThreshold(1.0e-3);
248         estimator.setMaxIterations(10);
249         estimator.setMaxEvaluations(20);
250 
251         KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 3,
252                                                    0.0, 5.6e-7,
253                                                    0.0, 7.7e-7,
254                                                    0.0, 8.7e-5,
255                                                    0.0, 3.5e-8);
256     }
257 
258     @Test
259     public void testWrappedException() {
260 
261         KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
262 
263         final KeplerianPropagatorBuilder propagatorBuilder =
264                         context.createBuilder(PositionAngle.MEAN, true, 1.0);
265 
266         // create perfect range measurements
267         final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
268                                                                                          propagatorBuilder);
269         final List<ObservedMeasurement<?>> measurements =
270                         KeplerianEstimationTestUtils.createMeasurements(propagator,
271                                                                              new RangeMeasurementCreator(context),
272                                                                              1.0, 3.0, 300.0);
273 
274         // create orbit estimator
275         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
276                                                                 propagatorBuilder);
277         for (final ObservedMeasurement<?> range : measurements) {
278             estimator.addMeasurement(range);
279         }
280         estimator.setParametersConvergenceThreshold(1.0e-2);
281         estimator.setMaxIterations(10);
282         estimator.setMaxEvaluations(20);
283         estimator.setObserver(new BatchLSObserver() {
284             /** {@inheritDoc} */
285             @Override
286             public void evaluationPerformed(int iterationsCount, int evaluationscount,
287                                            Orbit[] orbits,
288                                            ParameterDriversList estimatedOrbitalParameters,
289                                            ParameterDriversList estimatedPropagatorParameters,
290                                            ParameterDriversList estimatedMeasurementsParameters,
291                                            EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
292                 throw new DummyException();
293             }
294         });
295 
296         try {
297             KeplerianEstimationTestUtils.checkFit(context, estimator, 3, 4,
298                                                        0.0, 1.5e-6,
299                                                        0.0, 3.2e-6,
300                                                        0.0, 3.8e-7,
301                                                        0.0, 1.5e-10);
302             Assert.fail("an exception should have been thrown");
303         } catch (DummyException de) {
304             // expected
305         }
306 
307     }
308 
309     private static class DummyException extends OrekitException {
310         private static final long serialVersionUID = 1L;
311         public DummyException() {
312             super(OrekitMessages.INTERNAL_ERROR);
313         }
314     }
315 
316 }