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