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.ObservableSatellite;
34  import org.orekit.estimation.measurements.Position;
35  import org.orekit.estimation.measurements.modifiers.Bias;
36  import org.orekit.orbits.OrbitType;
37  import org.orekit.orbits.PositionAngleType;
38  import org.orekit.propagation.Propagator;
39  import org.orekit.propagation.SpacecraftState;
40  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
41  import org.orekit.time.AbsoluteDate;
42  import org.orekit.time.BurstSelector;
43  import org.orekit.time.TimeScalesFactory;
44  
45  import java.util.SortedSet;
46  
47  public class PositionBuilderTest {
48  
49      private static final double SIGMA = 10.0;
50      private static final double BIAS  =  5.0;
51  
52      private MeasurementBuilder<Position> getBuilder(final RandomGenerator random, final ObservableSatellite satellite) {
53          final RealMatrix covariance = MatrixUtils.createRealDiagonalMatrix(new double[] {
54              SIGMA * SIGMA, SIGMA * SIGMA, SIGMA * SIGMA
55          });
56          MeasurementBuilder<Position> pb =
57                          new PositionBuilder(random == null ? null : new CorrelatedRandomVectorGenerator(covariance,
58                                                                                                          1.0e-10,
59                                                                                                          new GaussianRandomGenerator(random)),
60                                              SIGMA, 1.0, satellite);
61           pb.addModifier(new Bias<>(new String[] { "pxBias", "pyBias", "pzBias" },
62                                     new double[] { BIAS, BIAS, BIAS },
63                                     new double[] { 1.0, 1.0, 1.0 },
64                                     new double[] { Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY },
65                                     new double[] { Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY }));
66           return pb;
67      }
68  
69      @Test
70      public void testForward() {
71          doTest(0x292b6e87436fe4c7L, 0.0, 1.2, 3.3 * SIGMA);
72      }
73  
74      @Test
75      public void testBackward() {
76          doTest(0x2f3285aa70b83c47L, 0.0, -1.0, 3.3 * SIGMA);
77      }
78  
79      private Propagator buildPropagator() {
80          return EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
81      }
82  
83      private void doTest(long seed, double startPeriod, double endPeriod, double tolerance) {
84          Generator    generator    = new Generator();
85          final int    maxBurstSize = 10;
86          final double highRateStep = 5.0;
87          final double burstPeriod  = 300.0;
88  
89          generator.addPropagator(buildPropagator()); // dummy first propagator
90          generator.addPropagator(buildPropagator()); // dummy second propagator
91          ObservableSatellite satellite = generator.addPropagator(buildPropagator()); // useful third propagator
92          generator.addPropagator(buildPropagator()); // dummy fourth propagator
93          generator.addScheduler(new ContinuousScheduler<>(getBuilder(new Well19937a(seed), satellite),
94                                                           new BurstSelector(maxBurstSize, highRateStep, burstPeriod,
95                                                                             TimeScalesFactory.getUTC())));
96          final GatheringSubscriber gatherer = new GatheringSubscriber();
97          generator.addSubscriber(gatherer);
98          final double period = context.initialOrbit.getKeplerianPeriod();
99          AbsoluteDate t0     = context.initialOrbit.getDate().shiftedBy(startPeriod * period);
100         AbsoluteDate t1     = context.initialOrbit.getDate().shiftedBy(endPeriod   * period);
101         generator.generate(t0, t1);
102         SortedSet<EstimatedMeasurementBase<?>> measurements = gatherer.getGeneratedMeasurements();
103         Propagator propagator = buildPropagator();
104         double maxError = 0;
105         AbsoluteDate previous = null;
106         AbsoluteDate tInf = t0.isBefore(t1) ? t0 : t1;
107         AbsoluteDate tSup = t0.isBefore(t1) ? t1 : t0;
108         int count = 0;
109         for (EstimatedMeasurementBase<?> measurement : measurements) {
110             AbsoluteDate date = measurement.getDate();
111             double[] m = measurement.getObservedValue();
112             Assertions.assertTrue(date.compareTo(tInf) >= 0);
113             Assertions.assertTrue(date.compareTo(tSup) <= 0);
114             if (previous != null) {
115                 // measurements are always chronological, even with backward propagation,
116                 // due to the SortedSet (which is intended for combining several
117                 // measurements types with different builders and schedulers)
118                 final double expected = (count % maxBurstSize == 0) ?
119                                         burstPeriod - (maxBurstSize - 1) * highRateStep :
120                                         highRateStep;
121                 if (t0.isBefore(t1)) {
122                     // measurements are expected to be chronological
123                     Assertions.assertEquals(expected, date.durationFrom(previous), 1.0e-10 * expected);
124                 } else {
125                     // measurements are expected to be reverse chronological
126                     Assertions.assertEquals(expected, previous.durationFrom(date), 1.0e-10 * expected);
127                 }
128             }
129             previous = date;
130             ++count;
131             SpacecraftState state = propagator.propagate(date);
132             double[] e = measurement.
133                 getObservedMeasurement().
134                 estimateWithoutDerivatives(new SpacecraftState[] { state }).
135                 getEstimatedValue();
136             for (int i = 0; i < m.length; ++i) {
137                 maxError = FastMath.max(maxError, FastMath.abs(e[i] - m[i]));
138             }
139         }
140         Assertions.assertEquals(0.0, maxError, tolerance);
141      }
142 
143      @BeforeEach
144      public void setUp() {
145          context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
146 
147          propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
148                                                    1.0e-6, 300.0, 0.001, Force.POTENTIAL,
149                                                    Force.THIRD_BODY_SUN, Force.THIRD_BODY_MOON);
150      }
151 
152      Context context;
153      NumericalPropagatorBuilder propagatorBuilder;
154 
155 }