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