1   /* Copyright 2022-2025 Luc Maisonobe
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.generation;
18  
19  import java.util.SortedSet;
20  
21  import org.hipparchus.linear.MatrixUtils;
22  import org.hipparchus.linear.RealMatrix;
23  import org.hipparchus.random.CorrelatedRandomVectorGenerator;
24  import org.hipparchus.random.GaussianRandomGenerator;
25  import org.hipparchus.random.RandomGenerator;
26  import org.hipparchus.random.Well19937a;
27  import org.hipparchus.util.FastMath;
28  import org.junit.jupiter.api.Assertions;
29  import org.junit.jupiter.api.BeforeEach;
30  import org.junit.jupiter.api.Test;
31  import org.orekit.estimation.Context;
32  import org.orekit.estimation.EstimationTestUtils;
33  import org.orekit.estimation.Force;
34  import org.orekit.estimation.measurements.EstimatedMeasurementBase;
35  import org.orekit.estimation.measurements.ObservableSatellite;
36  import org.orekit.estimation.measurements.gnss.AmbiguityCache;
37  import org.orekit.estimation.measurements.gnss.OneWayGNSSPhase;
38  import org.orekit.estimation.measurements.modifiers.Bias;
39  import org.orekit.gnss.PredefinedGnssSignal;
40  import org.orekit.orbits.KeplerianOrbit;
41  import org.orekit.orbits.Orbit;
42  import org.orekit.orbits.OrbitType;
43  import org.orekit.orbits.PositionAngleType;
44  import org.orekit.propagation.Propagator;
45  import org.orekit.propagation.SpacecraftState;
46  import org.orekit.propagation.analytical.KeplerianPropagator;
47  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
48  import org.orekit.propagation.events.InterSatDirectViewDetector;
49  import org.orekit.time.AbsoluteDate;
50  import org.orekit.time.FixedStepSelector;
51  import org.orekit.time.TimeScalesFactory;
52  import org.orekit.utils.PVCoordinates;
53  
54  public class OneWayGNSSPhaseBuilderTest {
55  
56      private static final double SIGMA =  0.5;
57      private static final double BIAS  = -0.01;
58      private static final double WAVELENGTH = PredefinedGnssSignal.G01.getWavelength();
59  
60      private MeasurementBuilder<OneWayGNSSPhase> getBuilder(final RandomGenerator random,
61                                                             final ObservableSatellite receiver,
62                                                             final ObservableSatellite remote) {
63          final RealMatrix covariance = MatrixUtils.createRealDiagonalMatrix(new double[] { SIGMA * SIGMA });
64          remote.getClockOffsetDriver().setValue(1.0e-16);
65          remote.getClockDriftDriver().setValue(0);
66          remote.getClockAccelerationDriver().setValue(0);
67          MeasurementBuilder<OneWayGNSSPhase> b =
68                          new OneWayGNSSPhaseBuilder(random == null ? null : new CorrelatedRandomVectorGenerator(covariance,
69                                                                                                                 1.0e-10,
70                                                                                                                 new GaussianRandomGenerator(random)),
71                                                     receiver, remote,
72                                                     WAVELENGTH, SIGMA, 1.0,
73                                                     new AmbiguityCache());
74          b.addModifier(new Bias<>(new String[] { "bias" },
75                           new double[] { BIAS },
76                           new double[] { 1.0 },
77                           new double[] { Double.NEGATIVE_INFINITY },
78                           new double[] { Double.POSITIVE_INFINITY }));
79          return b;
80      }
81  
82      @Test
83      public void testForward() {
84          doTest(0x812bfe784826bab3L, 0.0, 1.2, 3.4 * SIGMA);
85      }
86  
87      @Test
88      public void testBackward() {
89          doTest(0xb54896c361d88441L, 0.0, -1.0, 2.8 * SIGMA);
90      }
91  
92      private Propagator buildPropagator() {
93          return EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
94      }
95  
96      private void doTest(long seed, double startPeriod, double endPeriod, double tolerance) {
97          Generator generator = new Generator();
98          generator.addPropagator(buildPropagator()); // dummy first propagator
99          generator.addPropagator(buildPropagator()); // dummy second propagator
100         ObservableSatellite receiver = generator.addPropagator(buildPropagator()); // useful third propagator
101         generator.addPropagator(buildPropagator()); // dummy fourth propagator
102         final Orbit o1 = context.initialOrbit;
103         // for the second satellite, we simply reverse velocity
104         final Orbit o2 = new KeplerianOrbit(new PVCoordinates(o1.getPosition(),
105                                                               o1.getPVCoordinates().getVelocity().negate()),
106                                             o1.getFrame(), o1.getDate(), o1.getMu());
107         ObservableSatellite remote = generator.addPropagator(new KeplerianPropagator(o2)); // useful sixth propagator
108         final double step = 60.0;
109 
110         // beware that in order to avoid deadlocks, the secondary PV coordinates provider
111         // in InterSatDirectViewDetector must be *different* from the second propagator
112         // added to generator above! The reason is the event detector will be bound
113         // to the first propagator, so it cannot also refer to the second one at the same time
114         // this is the reason why we create a *new* KeplerianPropagator below
115         generator.addScheduler(new EventBasedScheduler<>(getBuilder(new Well19937a(seed), receiver, remote),
116                                                          new FixedStepSelector(step, TimeScalesFactory.getUTC()),
117                                                          generator.getPropagator(receiver),
118                                                          new InterSatDirectViewDetector(context.earth, new KeplerianPropagator(o2)),
119                                                          SignSemantic.FEASIBLE_MEASUREMENT_WHEN_POSITIVE));
120 
121         final GatheringSubscriber gatherer = new GatheringSubscriber();
122         generator.addSubscriber(gatherer);
123         final double period = o1.getKeplerianPeriod();
124         AbsoluteDate t0     = o1.getDate().shiftedBy(startPeriod * period);
125         AbsoluteDate t1     = o1.getDate().shiftedBy(endPeriod   * period);
126         generator.generate(t0, t1);
127         SortedSet<EstimatedMeasurementBase<?>> measurements = gatherer.getGeneratedMeasurements();
128 
129         // and yet another set of propagators for reference
130         Propagator propagator1 = buildPropagator();
131         Propagator propagator2 = new KeplerianPropagator(o2);
132 
133         double maxError = 0;
134         AbsoluteDate previous = null;
135         AbsoluteDate tInf = t0.isBefore(t1) ? t0 : t1;
136         AbsoluteDate tSup = t0.isBefore(t1) ? t1 : t0;
137         for (EstimatedMeasurementBase<?> measurement : measurements) {
138             AbsoluteDate date = measurement.getDate();
139             double[] m = measurement.getObservedValue();
140             Assertions.assertTrue(date.compareTo(tInf) >= 0);
141             Assertions.assertTrue(date.compareTo(tSup) <= 0);
142             if (previous != null) {
143                 if (t0.isBefore(t1)) {
144                     // measurements are expected to be chronological
145                     Assertions.assertTrue(date.durationFrom(previous) >= 0.999999 * step);
146                 } else {
147                     // measurements are expected to be reverse chronological
148                     Assertions.assertTrue(previous.durationFrom(date) >= 0.999999 * step);
149                 }
150             }
151             previous = date;
152             double[] e = measurement.
153                 getObservedMeasurement().
154                 estimateWithoutDerivatives(new SpacecraftState[] {
155                                                propagator1.propagate(date),
156                                                propagator2.propagate(date)
157                                            }).
158                 getEstimatedValue();
159             for (int i = 0; i < m.length; ++i) {
160                 maxError = FastMath.max(maxError, FastMath.abs(e[i] - m[i]));
161             }
162         }
163         Assertions.assertEquals(0.0, maxError, tolerance);
164      }
165 
166      @BeforeEach
167      public void setUp() {
168          context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
169 
170          propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
171                                                    1.0e-6, 300.0, 0.001, Force.POTENTIAL,
172                                                    Force.THIRD_BODY_SUN, Force.THIRD_BODY_MOON);
173      }
174 
175      Context context;
176      NumericalPropagatorBuilder propagatorBuilder;
177 
178 }