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.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.Assert;
29  import org.junit.Before;
30  import org.junit.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.InterSatellitesRange;
35  import org.orekit.estimation.measurements.ObservableSatellite;
36  import org.orekit.estimation.measurements.ObservedMeasurement;
37  import org.orekit.estimation.measurements.modifiers.Bias;
38  import org.orekit.orbits.KeplerianOrbit;
39  import org.orekit.orbits.Orbit;
40  import org.orekit.orbits.OrbitType;
41  import org.orekit.orbits.PositionAngle;
42  import org.orekit.propagation.Propagator;
43  import org.orekit.propagation.SpacecraftState;
44  import org.orekit.propagation.analytical.KeplerianPropagator;
45  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
46  import org.orekit.propagation.events.InterSatDirectViewDetector;
47  import org.orekit.time.AbsoluteDate;
48  import org.orekit.time.FixedStepSelector;
49  import org.orekit.time.TimeScalesFactory;
50  import org.orekit.utils.PVCoordinates;
51  
52  public class InterSatellitesRangeBuilderTest {
53  
54      private static final double SIGMA =  0.5;
55      private static final double BIAS  = -0.01;
56  
57      private MeasurementBuilder<InterSatellitesRange> getBuilder(final RandomGenerator random,
58                                                                  final ObservableSatellite receiver,
59                                                                  final ObservableSatellite remote) {
60          final RealMatrix covariance = MatrixUtils.createRealDiagonalMatrix(new double[] { SIGMA * SIGMA });
61          MeasurementBuilder<InterSatellitesRange> isrb =
62                          new InterSatellitesRangeBuilder(random == null ? null : new CorrelatedRandomVectorGenerator(covariance,
63                                                                                                                      1.0e-10,
64                                                                                                                      new GaussianRandomGenerator(random)),
65                                                          receiver, remote, true, SIGMA, 1.0);
66          isrb.addModifier(new Bias<>(new String[] { "bias" },
67                           new double[] { BIAS },
68                           new double[] { 1.0 },
69                           new double[] { Double.NEGATIVE_INFINITY },
70                           new double[] { Double.POSITIVE_INFINITY }));
71          return isrb;
72      }
73  
74      @Test
75      public void testForward() {
76          doTest(0xc82a56322345dc25l, 0.0, 1.2, 2.8 * SIGMA);
77      }
78  
79      @Test
80      public void testBackward() {
81          doTest(0x95c10149c4891232l, 0.0, -1.0, 2.6 * SIGMA);
82      }
83  
84      private Propagator buildPropagator() {
85          return EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
86      }
87  
88      private void doTest(long seed, double startPeriod, double endPeriod, double tolerance) {
89          Generator generator = new Generator();
90          generator.addPropagator(buildPropagator()); // dummy first propagator
91          generator.addPropagator(buildPropagator()); // dummy second propagator
92          ObservableSatellite receiver = generator.addPropagator(buildPropagator()); // useful third propagator
93          generator.addPropagator(buildPropagator()); // dummy fourth propagator
94          final Orbit o1 = context.initialOrbit;
95          // for the second satellite, we simply reverse velocity
96          final Orbit o2 = new KeplerianOrbit(new PVCoordinates(o1.getPVCoordinates().getPosition(),
97                                                                o1.getPVCoordinates().getVelocity().negate()),
98                                              o1.getFrame(), o1.getDate(), o1.getMu());
99          ObservableSatellite remote = generator.addPropagator(new KeplerianPropagator(o2)); // useful sixth propagator
100         final double step = 60.0;
101 
102         // beware that in order to avoid deadlocks, the secondary PV coordinates provider
103         // in InterSatDirectViewDetector must be *different* from the second propagator
104         // added to generator above! The reason is the event detector will be bound
105         // to the first propagator, so it cannot also refer to the second one at the same time
106         // this is the reason why we create a *new* KeplerianPropagator below
107         generator.addScheduler(new EventBasedScheduler<>(getBuilder(new Well19937a(seed), receiver, remote),
108                                                          new FixedStepSelector(step, TimeScalesFactory.getUTC()),
109                                                          generator.getPropagator(receiver),
110                                                          new InterSatDirectViewDetector(context.earth, new KeplerianPropagator(o2)),
111                                                          SignSemantic.FEASIBLE_MEASUREMENT_WHEN_POSITIVE));
112 
113         final double period = o1.getKeplerianPeriod();
114         AbsoluteDate t0     = o1.getDate().shiftedBy(startPeriod * period);
115         AbsoluteDate t1     = o1.getDate().shiftedBy(endPeriod   * period);
116         SortedSet<ObservedMeasurement<?>> measurements = generator.generate(t0, t1);
117 
118         // and yet another set of propagators for reference
119         Propagator propagator1 = buildPropagator();
120         Propagator propagator2 = new KeplerianPropagator(o2);
121 
122         double maxError = 0;
123         AbsoluteDate previous = null;
124         AbsoluteDate tInf = t0.compareTo(t1) < 0 ? t0 : t1;
125         AbsoluteDate tSup = t0.compareTo(t1) < 0 ? t1 : t0;
126         for (ObservedMeasurement<?> measurement : measurements) {
127             AbsoluteDate date = measurement.getDate();
128             double[] m = measurement.getObservedValue();
129             Assert.assertTrue(date.compareTo(tInf) >= 0);
130             Assert.assertTrue(date.compareTo(tSup) <= 0);
131             if (previous != null) {
132                 // measurements are always chronological, even with backward propagation,
133                 // due to the SortedSet (which is intended for combining several
134                 // measurements types with different builders and schedulers)
135                 Assert.assertTrue(date.durationFrom(previous) >= 0.999999 * step);
136             }
137             previous = date;
138             double[] e = measurement.estimate(0, 0,
139                                               new SpacecraftState[] {
140                                                   propagator1.propagate(date),
141                                                   propagator2.propagate(date)
142                                               }).getEstimatedValue();
143             for (int i = 0; i < m.length; ++i) {
144                 maxError = FastMath.max(maxError, FastMath.abs(e[i] - m[i]));
145             }
146         }
147         Assert.assertEquals(0.0, maxError, tolerance);
148      }
149 
150      @Before
151      public void setUp() {
152          context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
153 
154          propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
155                                                    1.0e-6, 300.0, 0.001, Force.POTENTIAL,
156                                                    Force.THIRD_BODY_SUN, Force.THIRD_BODY_MOON);
157      }
158 
159      Context context;
160      NumericalPropagatorBuilder propagatorBuilder;
161 
162 }