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.measurements.modifiers;
18  
19  import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
20  import org.hipparchus.random.RandomGenerator;
21  import org.hipparchus.random.Well19937a;
22  import org.junit.jupiter.api.Assertions;
23  import org.junit.jupiter.api.Test;
24  import org.orekit.errors.OrekitException;
25  import org.orekit.errors.OrekitMessages;
26  import org.orekit.estimation.Context;
27  import org.orekit.estimation.EstimationTestUtils;
28  import org.orekit.estimation.leastsquares.BatchLSEstimator;
29  import org.orekit.estimation.measurements.ObservedMeasurement;
30  import org.orekit.estimation.measurements.Range;
31  import org.orekit.estimation.measurements.TwoWayRangeMeasurementCreator;
32  import org.orekit.frames.TopocentricFrame;
33  import org.orekit.orbits.OrbitType;
34  import org.orekit.orbits.PositionAngleType;
35  import org.orekit.propagation.Propagator;
36  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
37  import org.orekit.utils.ParameterDriver;
38  
39  import java.util.List;
40  
41  public class BiasTest {
42  
43      @SuppressWarnings("unchecked")
44      @Test
45      public void testEstimateBias() {
46  
47          Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
48  
49          final NumericalPropagatorBuilder propagatorBuilder =
50                          context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
51                                                1.0e-6, 60.0, 0.001);
52  
53          // create perfect range measurements
54          final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
55                                                                             propagatorBuilder);
56          final List<ObservedMeasurement<?>> measurements =
57                          EstimationTestUtils.createMeasurements(propagator,
58                                                                 new TwoWayRangeMeasurementCreator(context),
59                                                                 1.0, 3.0, 300.0);
60  
61          // create range biases: one bias for each station
62          final RandomGenerator random = new Well19937a(0x0c4b69da5d64b35aL);
63          final Bias<?>[] stationsRangeBiases = new Bias<?>[context.stations.size()];
64          final double[] realStationsBiases  = new double[context.stations.size()];
65          for (int i = 0; i < context.stations.size(); ++i) {
66              final TopocentricFrame base = context.stations.get(i).getBaseFrame();
67              stationsRangeBiases[i] = new Bias<Range>(new String[] {
68                                                           base.getName() + " range bias"
69                                                       },
70                                                       new double[] {
71                                                           0.0
72                                                       },
73                                                       new double[] {
74                                                           1.0
75                                                       },
76                                                       new double[] {
77                                                           Double.NEGATIVE_INFINITY
78                                                       },
79                                                       new double[] {
80                                                           Double.POSITIVE_INFINITY
81                                                       });
82              realStationsBiases[i]  = 2 * random.nextDouble() - 1;
83              Assertions.assertEquals(base.getName() + " range bias", stationsRangeBiases[i].getEffectName());
84          }
85  
86          // create orbit estimator
87          final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
88                                                                  propagatorBuilder);
89  
90          // add the measurements, with both spacecraft and stations biases
91          for (final ObservedMeasurement<?> measurement : measurements) {
92              final Range range = (Range) measurement;
93              for (int i = 0; i < context.stations.size(); ++i) {
94                  if (range.getStation() == context.stations.get(i)) {
95                      double biasedRange = range.getObservedValue()[0] + realStationsBiases[i];
96                      final Range m = new Range(range.getStation(),
97                                                range.isTwoWay(),
98                                                range.getDate(),
99                                                biasedRange,
100                                               range.getTheoreticalStandardDeviation()[0],
101                                               range.getBaseWeight()[0],
102                                               range.getSatellites().get(0));
103                     m.addModifier((Bias<Range>) stationsRangeBiases[i]);
104                     estimator.addMeasurement(m);
105                 }
106             }
107         }
108 
109         estimator.setParametersConvergenceThreshold(1.0e-3);
110         estimator.setMaxIterations(10);
111         estimator.setMaxEvaluations(20);
112 
113         // we want to estimate the biases
114         for (Bias<?> bias : stationsRangeBiases) {
115             for (final ParameterDriver driver : bias.getParametersDrivers()) {
116                 driver.setSelected(true);
117             }
118         }
119 
120         EstimationTestUtils.checkFit(context, estimator, 2, 3,
121                                      0.0,  7.2e-7,
122                                      0.0,  2.1e-6,
123                                      0.0,  6e-7,
124                                      0.0,  2.5e-10);
125         for (int i = 0; i < stationsRangeBiases.length; ++i) {
126             Assertions.assertEquals(realStationsBiases[i],
127                                 stationsRangeBiases[i].getParametersDrivers().get(0).getValue(),
128                                 3.3e-6);
129         }
130 
131     }
132 
133     @Test
134     public void testTooSmallScale() {
135         try {
136             new Bias<Range>(new String[] { "OK", "not-OK" },
137                             new double[] { 1000.0,    1000.0 },
138                             new double[] {    1.0,    0.0 },
139                             new double[] { -10000.0, -10000.0 },
140                             new double[] { +10000.0, +10000.0 });
141             Assertions.fail("an exception should have been thrown");
142         } catch (OrekitException oe) {
143             Assertions.assertEquals(OrekitMessages.TOO_SMALL_SCALE_FOR_PARAMETER, oe.getSpecifier());
144             Assertions.assertEquals("not-OK", oe.getParts()[0]);
145         }
146     }
147 }
148 
149