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;
18  
19  import org.hipparchus.analysis.UnivariateVectorFunction;
20  import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
21  import org.hipparchus.analysis.differentiation.Gradient;
22  import org.hipparchus.analysis.differentiation.GradientField;
23  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction;
24  import org.hipparchus.geometry.euclidean.threed.Rotation;
25  import org.hipparchus.geometry.euclidean.threed.Vector3D;
26  import org.hipparchus.linear.RealMatrix;
27  import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
28  import org.hipparchus.random.RandomGenerator;
29  import org.hipparchus.random.Well19937a;
30  import org.hipparchus.util.FastMath;
31  import org.hipparchus.util.Precision;
32  import org.junit.jupiter.api.Assertions;
33  import org.junit.jupiter.api.Test;
34  import org.orekit.Utils;
35  import org.orekit.bodies.BodyShape;
36  import org.orekit.bodies.GeodeticPoint;
37  import org.orekit.bodies.OneAxisEllipsoid;
38  import org.orekit.errors.OrekitException;
39  import org.orekit.errors.OrekitMessages;
40  import org.orekit.estimation.Context;
41  import org.orekit.estimation.EstimationTestUtils;
42  import org.orekit.estimation.leastsquares.BatchLSEstimator;
43  import org.orekit.frames.FieldTransform;
44  import org.orekit.frames.Frame;
45  import org.orekit.frames.FramesFactory;
46  import org.orekit.frames.StaticTransform;
47  import org.orekit.frames.TopocentricFrame;
48  import org.orekit.frames.Transform;
49  import org.orekit.orbits.OrbitType;
50  import org.orekit.orbits.PositionAngleType;
51  import org.orekit.propagation.Propagator;
52  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
53  import org.orekit.time.AbsoluteDate;
54  import org.orekit.utils.Constants;
55  import org.orekit.utils.FieldPVCoordinates;
56  import org.orekit.utils.IERSConventions;
57  import org.orekit.utils.PVCoordinates;
58  import org.orekit.utils.ParameterDriver;
59  
60  import java.util.Arrays;
61  import java.util.HashMap;
62  import java.util.List;
63  import java.util.Locale;
64  import java.util.Map;
65  
66  class GroundStationTest {
67  
68      @Test
69      void testEstimateClockOffset() {
70  
71          Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
72  
73          final NumericalPropagatorBuilder propagatorBuilder =
74                          context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
75                                                1.0e-6, 60.0, 0.001);
76  
77          // create perfect range measurements
78          final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
79                                                                             propagatorBuilder);
80          final List<ObservedMeasurement<?>> measurements =
81                          EstimationTestUtils.createMeasurements(propagator,
82                                                                 new TwoWayRangeMeasurementCreator(context),
83                                                                 1.0, 3.0, 300.0);
84  
85          // change one station clock
86          final TopocentricFrame base  = context.stations.get(0).getBaseFrame();
87          final BodyShape parent       = base.getParentShape();
88          final double deltaClock      = 0.00084532;
89          final String changedSuffix   = "-changed";
90          final GroundStation changed  = new GroundStation(new TopocentricFrame(parent, base.getPoint(),
91                                                                                base.getName() + changedSuffix), context.ut1.getEOPHistory(),
92                                                           context.stations.get(0).getDisplacements());
93  
94          // create orbit estimator
95          final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
96                                                                  propagatorBuilder);
97          for (final ObservedMeasurement<?> measurement : measurements) {
98              final Range range = (Range) measurement;
99              final String name = range.getStation().getBaseFrame().getName() + changedSuffix;
100                 if (changed.getBaseFrame().getName().equals(name)) {
101                     estimator.addMeasurement(new Range(changed, range.isTwoWay(),
102                                                        range.getDate().shiftedBy(deltaClock),
103                                                        range.getObservedValue()[0],
104                                                        range.getTheoreticalStandardDeviation()[0],
105                                                        range.getBaseWeight()[0],
106                                                        range.getSatellites().get(0)));
107                 } else {
108                     estimator.addMeasurement(range);
109                 }
110         }
111         estimator.setParametersConvergenceThreshold(1.0e-3);
112         estimator.setMaxIterations(100);
113         estimator.setMaxEvaluations(200);
114 
115         // we want to estimate station clock offset
116         changed.getClockOffsetDriver().setSelected(true);
117         changed.getEastOffsetDriver().setSelected(false);
118         changed.getNorthOffsetDriver().setSelected(false);
119         changed.getZenithOffsetDriver().setSelected(false);
120 
121         EstimationTestUtils.checkFit(context, estimator, 2, 3,
122                                      0.0, 6.8e-7,
123                                      0.0, 2.0e-6,
124                                      0.0, 1.8e-7,
125                                      0.0, 8e-11);
126         Assertions.assertEquals(deltaClock, changed.getClockOffsetDriver().getValue(), 9.6e-11);
127 
128         RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
129         RealMatrix physicalCovariances   = estimator.getPhysicalCovariances(1.0e-10);
130         Assertions.assertEquals(7,        normalizedCovariances.getRowDimension());
131         Assertions.assertEquals(7,        normalizedCovariances.getColumnDimension());
132         Assertions.assertEquals(7,        physicalCovariances.getRowDimension());
133         Assertions.assertEquals(7,        physicalCovariances.getColumnDimension());
134         Assertions.assertEquals(4.185e-9, physicalCovariances.getEntry(6, 6), 3.0e-13);
135 
136     }
137 
138     @Test
139     void testEstimateStationPosition() {
140 
141         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
142 
143         final NumericalPropagatorBuilder propagatorBuilder =
144                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
145                                               1.0e-6, 60.0, 0.001);
146 
147         // create perfect range measurements
148         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
149                                                                            propagatorBuilder);
150         final List<ObservedMeasurement<?>> measurements =
151                         EstimationTestUtils.createMeasurements(propagator,
152                                                                new TwoWayRangeMeasurementCreator(context),
153                                                                1.0, 3.0, 300.0);
154 
155         // move one station
156         final RandomGenerator random = new Well19937a(0x4adbecfc743bda60L);
157         final TopocentricFrame base = context.stations.get(0).getBaseFrame();
158         final BodyShape parent = base.getParentShape();
159         final Vector3D baseOrigin = parent.transform(base.getPoint());
160         final Vector3D deltaTopo = new Vector3D(2 * random.nextDouble() - 1,
161                                                 2 * random.nextDouble() - 1,
162                                                 2 * random.nextDouble() - 1);
163         final StaticTransform topoToParent = base.getStaticTransformTo(parent.getBodyFrame(), (AbsoluteDate) null);
164         final Vector3D deltaParent   = topoToParent.transformVector(deltaTopo);
165         final String movedSuffix     = "-moved";
166         final GroundStation moved = new GroundStation(new TopocentricFrame(parent,
167                                                                            parent.transform(baseOrigin.subtract(deltaParent),
168                                                                                             parent.getBodyFrame(),
169                                                                                             null),
170                                                                            base.getName() + movedSuffix), context.ut1.getEOPHistory(),
171                                                       context.stations.get(0).getDisplacements());
172 
173         // create orbit estimator
174         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
175                                                                 propagatorBuilder);
176         for (final ObservedMeasurement<?> measurement : measurements) {
177             final Range range = (Range) measurement;
178             final String name = range.getStation().getBaseFrame().getName() + movedSuffix;
179                 if (moved.getBaseFrame().getName().equals(name)) {
180                     estimator.addMeasurement(new Range(moved, range.isTwoWay(), range.getDate(),
181                                                        range.getObservedValue()[0],
182                                                        range.getTheoreticalStandardDeviation()[0],
183                                                        range.getBaseWeight()[0],
184                                                        range.getSatellites().get(0)));
185                 } else {
186                     estimator.addMeasurement(range);
187                 }
188         }
189         estimator.setParametersConvergenceThreshold(1.0e-3);
190         estimator.setMaxIterations(100);
191         estimator.setMaxEvaluations(200);
192 
193         // we want to estimate station offsets
194         moved.getClockOffsetDriver().setSelected(false);
195         moved.getEastOffsetDriver().setSelected(true);
196         moved.getNorthOffsetDriver().setSelected(true);
197         moved.getZenithOffsetDriver().setSelected(true);
198 
199         EstimationTestUtils.checkFit(context, estimator, 2, 7,
200                                      0.0, 6.3e-7,
201                                      0.0, 1.8e-6,
202                                      0.0, 9e-7,
203                                      0.0, 5e-10);
204         Assertions.assertEquals(deltaTopo.getX(), moved.getEastOffsetDriver().getValue(),   4.5e-7);
205         Assertions.assertEquals(deltaTopo.getY(), moved.getNorthOffsetDriver().getValue(),  6.2e-7);
206         Assertions.assertEquals(deltaTopo.getZ(), moved.getZenithOffsetDriver().getValue(), 2.6e-7);
207 
208         GeodeticPoint result = moved.getOffsetGeodeticPoint((AbsoluteDate) null);
209 
210         GeodeticPoint reference = context.stations.get(0).getBaseFrame().getPoint();
211         Assertions.assertEquals(reference.getLatitude(),  result.getLatitude(),  3.3e-14);
212         Assertions.assertEquals(reference.getLongitude(), result.getLongitude(), 1e-13);
213         Assertions.assertEquals(reference.getAltitude(),  result.getAltitude(),  2.6e-7);
214 
215         RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
216         RealMatrix physicalCovariances   = estimator.getPhysicalCovariances(1.0e-10);
217         Assertions.assertEquals(9,       normalizedCovariances.getRowDimension());
218         Assertions.assertEquals(9,       normalizedCovariances.getColumnDimension());
219         Assertions.assertEquals(9,       physicalCovariances.getRowDimension());
220         Assertions.assertEquals(9,       physicalCovariances.getColumnDimension());
221         Assertions.assertEquals(0.55431, physicalCovariances.getEntry(6, 6), 1.0e-5);
222         Assertions.assertEquals(0.22694, physicalCovariances.getEntry(7, 7), 1.0e-5);
223         Assertions.assertEquals(0.13106, physicalCovariances.getEntry(8, 8), 1.0e-5);
224     }
225 
226     @Test
227     void testEstimateEOP() {
228 
229         Context linearEOPContext = EstimationTestUtils.eccentricContext("linear-EOP:regular-data/de431-ephemerides:potential:tides");
230 
231         final AbsoluteDate refDate = new AbsoluteDate(2000, 2, 24, linearEOPContext.utc);
232         final double dut10 = 0.3079738;
233         final double lod   = 0.0011000;
234         final double xp0   = 68450.0e-6;
235         final double xpDot =   -50.0e-6;
236         final double yp0   =    60.0e-6;
237         final double ypDot =     2.0e-6;
238         for (double dt = -2 * Constants.JULIAN_DAY; dt < 2 * Constants.JULIAN_DAY; dt += 300.0) {
239             AbsoluteDate date = refDate.shiftedBy(dt);
240             Assertions.assertEquals(dut10 - dt * lod / Constants.JULIAN_DAY,
241                                 linearEOPContext.ut1.getEOPHistory().getUT1MinusUTC(date),
242                                 1.0e-15);
243             Assertions.assertEquals(lod,
244                                 linearEOPContext.ut1.getEOPHistory().getLOD(date),
245                                 1.0e-15);
246             Assertions.assertEquals((xp0 + xpDot * dt / Constants.JULIAN_DAY) * Constants.ARC_SECONDS_TO_RADIANS,
247                                 linearEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getXp(),
248                                 1.0e-15);
249             Assertions.assertEquals((yp0 + ypDot * dt / Constants.JULIAN_DAY) * Constants.ARC_SECONDS_TO_RADIANS,
250                                 linearEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getYp(),
251                                 1.0e-15);
252         }
253         final NumericalPropagatorBuilder linearPropagatorBuilder =
254                         linearEOPContext.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
255                                               1.0e-6, 60.0, 0.001);
256 
257         // create perfect range measurements
258         final Propagator propagator = EstimationTestUtils.createPropagator(linearEOPContext.initialOrbit,
259                                                                            linearPropagatorBuilder);
260         final List<ObservedMeasurement<?>> linearMeasurements =
261                         EstimationTestUtils.createMeasurements(propagator,
262                                                                new TwoWayRangeMeasurementCreator(linearEOPContext),
263                                                                1.0, 5.0, 60.0);
264 
265         Utils.clearFactories();
266         Context zeroEOPContext = EstimationTestUtils.eccentricContext("zero-EOP:regular-data/de431-ephemerides:potential:potential:tides");
267         for (double dt = -2 * Constants.JULIAN_DAY; dt < 2 * Constants.JULIAN_DAY; dt += 300.0) {
268             AbsoluteDate date = refDate.shiftedBy(dt);
269             Assertions.assertEquals(0.0,
270                                 zeroEOPContext.ut1.getEOPHistory().getUT1MinusUTC(date),
271                                 1.0e-15);
272             Assertions.assertEquals(0.0,
273                                 zeroEOPContext.ut1.getEOPHistory().getLOD(date),
274                                 1.0e-15);
275             Assertions.assertEquals(0.0,
276                                 zeroEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getXp(),
277                                 1.0e-15);
278             Assertions.assertEquals(0.0,
279                                 zeroEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getYp(),
280                                 1.0e-15);
281         }
282 
283         // create orbit estimator
284         final NumericalPropagatorBuilder zeroPropagatorBuilder =
285                         linearEOPContext.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
286                                               1.0e-6, 60.0, 0.001);
287         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
288                                                                 zeroPropagatorBuilder);
289         for (final ObservedMeasurement<?> linearMeasurement : linearMeasurements) {
290             Range linearRange = (Range) linearMeasurement;
291             for (final GroundStation station : zeroEOPContext.stations) {
292                 if (station.getBaseFrame().getName().equals(linearRange.getStation().getBaseFrame().getName())) {
293                     Range zeroRange = new Range(station, linearRange.isTwoWay(),
294                                                 linearRange.getDate(),
295                                                 linearRange.getObservedValue()[0],
296                                                 linearRange.getTheoreticalStandardDeviation()[0],
297                                                 linearRange.getBaseWeight()[0],
298                                                 linearRange.getSatellites().get(0));
299                     estimator.addMeasurement(zeroRange);
300                 }
301             }
302         }
303         estimator.setParametersConvergenceThreshold(1.0e-3);
304         estimator.setMaxIterations(100);
305         estimator.setMaxEvaluations(200);
306 
307         // we want to estimate pole and prime meridian
308         GroundStation station = zeroEOPContext.stations.get(0);
309         station.getPrimeMeridianOffsetDriver().setReferenceDate(refDate);
310         station.getPrimeMeridianOffsetDriver().setSelected(true);
311         station.getPrimeMeridianDriftDriver().setSelected(true);
312         station.getPolarOffsetXDriver().setReferenceDate(refDate);
313         station.getPolarOffsetXDriver().setSelected(true);
314         station.getPolarDriftXDriver().setSelected(true);
315         station.getPolarOffsetYDriver().setReferenceDate(refDate);
316         station.getPolarOffsetYDriver().setSelected(true);
317         station.getPolarDriftYDriver().setSelected(true);
318 
319         // just for the fun and to speed up test, we will use orbit determination, *without* estimating orbit
320         for (final ParameterDriver driver : zeroPropagatorBuilder.getOrbitalParametersDrivers().getDrivers()) {
321             driver.setSelected(false);
322         }
323 
324         estimator.estimate();
325 
326         final double computedDut1  = station.getPrimeMeridianOffsetDriver().getValue() / EstimatedEarthFrameProvider.EARTH_ANGULAR_VELOCITY;
327         final double computedLOD   = station.getPrimeMeridianDriftDriver().getValue() * (-Constants.JULIAN_DAY / EstimatedEarthFrameProvider.EARTH_ANGULAR_VELOCITY);
328         final double computedXp    = station.getPolarOffsetXDriver().getValue() / Constants.ARC_SECONDS_TO_RADIANS;
329         final double computedXpDot = station.getPolarDriftXDriver().getValue()  / Constants.ARC_SECONDS_TO_RADIANS * Constants.JULIAN_DAY;
330         final double computedYp    = station.getPolarOffsetYDriver().getValue() / Constants.ARC_SECONDS_TO_RADIANS;
331         final double computedYpDot = station.getPolarDriftYDriver().getValue()  / Constants.ARC_SECONDS_TO_RADIANS * Constants.JULIAN_DAY;
332         Assertions.assertEquals(0.0, FastMath.abs(dut10 - computedDut1),  4.3e-10);
333         Assertions.assertEquals(0.0, FastMath.abs(lod - computedLOD),     4.9e-10);
334         Assertions.assertEquals(0.0, FastMath.abs(xp0 - computedXp),      5.9e-9);
335         Assertions.assertEquals(0.0, FastMath.abs(xpDot - computedXpDot), 7.5e-9);
336         Assertions.assertEquals(0.0, FastMath.abs(yp0 - computedYp),      1.1e-9);
337         Assertions.assertEquals(0.0, FastMath.abs(ypDot - computedYpDot), 1.1e-10);
338 
339         // thresholds to use if orbit is estimated
340         // (i.e. when commenting out the loop above that sets orbital parameters drivers to "not selected")
341 //         Assertions.assertEquals(dut10, computedDut1,  6.6e-3);
342 //         Assertions.assertEquals(lod,   computedLOD,   1.1e-9);
343 //         Assertions.assertEquals(xp0,   computedXp,    3.3e-8);
344 //         Assertions.assertEquals(xpDot, computedXpDot, 2.2e-8);
345 //         Assertions.assertEquals(yp0,   computedYp,    3.3e-8);
346 //         Assertions.assertEquals(ypDot, computedYpDot, 3.8e-8);
347 
348     }
349 
350     @Test
351     void testClockOffsetCartesianDerivativesOctantPxPyPz() {
352         double relativeTolerancePositionValue      =  2.3e-15;
353         double relativeTolerancePositionDerivative =  2.5e-10;
354         double relativeToleranceVelocityValue      =  3.0e-15;
355         double relativeToleranceVelocityDerivative =  1.7e-10;
356         doTestCartesianDerivatives(FastMath.toRadians(35), FastMath.toRadians(20), 1200.0, 100.0,
357                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
358                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
359                                    ".*-clock");
360     }
361 
362     @Test
363     void testClockOffsetAngularDerivativesOctantPxPyPz() {
364         double toleranceRotationValue              =  1.9e-15;
365         double toleranceRotationDerivative         =  4.9e-15;
366         double toleranceRotationRateValue          =  1.1e-19;
367         double toleranceRotationRateDerivative     =  6.3e-19;
368         doTestAngularDerivatives(FastMath.toRadians(35), FastMath.toRadians(20), 1200.0, 100.0,
369                                  toleranceRotationValue,     toleranceRotationDerivative,
370                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
371                                  ".*-clock");
372     }
373 
374     @Test
375     void testClockOffsetCartesianDerivativesOctantPxPyMz() {
376         double relativeTolerancePositionValue      =  1.4e-15;
377         double relativeTolerancePositionDerivative =  1.7e-10;
378         double relativeToleranceVelocityValue      =  2.5e-15;
379         double relativeToleranceVelocityDerivative =  1.8e-10;
380         doTestCartesianDerivatives(FastMath.toRadians(-35), FastMath.toRadians(20), 1200.0, 100.0,
381                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
382                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
383                                    ".*-clock");
384     }
385 
386     @Test
387     void testClockOffsetAngularDerivativesOctantPxPyMz() {
388         double toleranceRotationValue              =  1.7e-15;
389         double toleranceRotationDerivative         =  5.0e-15;
390         double toleranceRotationRateValue          =  1.1e-19;
391         double toleranceRotationRateDerivative     =  6.3e-19;
392         doTestAngularDerivatives(FastMath.toRadians(-35), FastMath.toRadians(20), 1200.0, 100.0,
393                                  toleranceRotationValue,     toleranceRotationDerivative,
394                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
395                                  ".*-clock");
396     }
397 
398     @Test
399     void testClockOffsetCartesianDerivativesOctantPxMyPz() {
400         double relativeTolerancePositionValue      =  1.7e-15;
401         double relativeTolerancePositionDerivative =  2.6e-10;
402         double relativeToleranceVelocityValue      =  2.8e-15;
403         double relativeToleranceVelocityDerivative =  1.8e-10;
404         doTestCartesianDerivatives(FastMath.toRadians(35), FastMath.toRadians(-20), 1200.0, 100.0,
405                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
406                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
407                                    ".*-clock");
408     }
409 
410     @Test
411     void testClockOffsetAngularDerivativesOctantPxMyPz() {
412         double toleranceRotationValue              =  1.6e-15;
413         double toleranceRotationDerivative         =  4.9e-15;
414         double toleranceRotationRateValue          =  1.1e-19;
415         double toleranceRotationRateDerivative     =  6.3e-19;
416         doTestAngularDerivatives(FastMath.toRadians(35), FastMath.toRadians(-20), 1200.0, 100.0,
417                                  toleranceRotationValue,     toleranceRotationDerivative,
418                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
419                                  ".*-clock");
420     }
421 
422     @Test
423     void testClockOffsetCartesianDerivativesOctantPxMyMz() {
424         double relativeTolerancePositionValue      =  1.5e-15;
425         double relativeTolerancePositionDerivative =  1.6e-10;
426         double relativeToleranceVelocityValue      =  2.3e-15;
427         double relativeToleranceVelocityDerivative =  1.7e-10;
428         doTestCartesianDerivatives(FastMath.toRadians(-35), FastMath.toRadians(-20), 1200.0, 100.0,
429                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
430                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
431                                    ".*-clock");
432     }
433 
434     @Test
435     void testClockOffsetAngularDerivativesOctantPxMyMz() {
436         double toleranceRotationValue              =  1.5e-15;
437         double toleranceRotationDerivative         =  5.0e-15;
438         double toleranceRotationRateValue          =  1.1e-19;
439         double toleranceRotationRateDerivative     =  6.3e-19;
440         doTestAngularDerivatives(FastMath.toRadians(-35), FastMath.toRadians(-20), 1200.0, 100.0,
441                                  toleranceRotationValue,     toleranceRotationDerivative,
442                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
443                                  ".*-clock");
444     }
445 
446     @Test
447     void testClockOffsetCartesianDerivativesOctantMxPyPz() {
448         double relativeTolerancePositionValue      =  1.9e-15;
449         double relativeTolerancePositionDerivative =  2.6e-10;
450         double relativeToleranceVelocityValue      =  2.6e-15;
451         double relativeToleranceVelocityDerivative =  2.0e-10;
452         doTestCartesianDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 100.0,
453                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
454                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
455                                    ".*-clock");
456     }
457 
458     @Test
459     void testClockOffsetAngularDerivativesOctantMxPyPz() {
460         double toleranceRotationValue              =  1.4e-15;
461         double toleranceRotationDerivative         =  5.2e-15;
462         double toleranceRotationRateValue          =  1.1e-19;
463         double toleranceRotationRateDerivative     =  6.3e-19;
464         doTestAngularDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 100.0,
465                                  toleranceRotationValue,     toleranceRotationDerivative,
466                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
467                                  ".*-clock");
468     }
469 
470     @Test
471     void testClockOffsetCartesianDerivativesOctantMxPyMz() {
472         double relativeTolerancePositionValue      =  1.9e-15;
473         double relativeTolerancePositionDerivative =  2.6e-10;
474         double relativeToleranceVelocityValue      =  2.6e-15;
475         double relativeToleranceVelocityDerivative =  2.0e-10;
476         doTestCartesianDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 100.0,
477                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
478                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
479                                    ".*-clock");
480     }
481 
482     @Test
483     void testClockOffsetAngularDerivativesOctantMxPyMz() {
484         double toleranceRotationValue              =  1.4e-15;
485         double toleranceRotationDerivative         =  5.2e-15;
486         double toleranceRotationRateValue          =  1.1e-19;
487         double toleranceRotationRateDerivative     =  6.3e-19;
488         doTestAngularDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 100.0,
489                                  toleranceRotationValue,     toleranceRotationDerivative,
490                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
491                                  ".*-clock");
492     }
493 
494     @Test
495     void testClockOffsetCartesianDerivativesOctantMxMyPz() {
496         double relativeTolerancePositionValue      =  1.5e-15;
497         double relativeTolerancePositionDerivative =  1.7e-10;
498         double relativeToleranceVelocityValue      =  2.9e-15;
499         double relativeToleranceVelocityDerivative =  1.9e-10;
500         doTestCartesianDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 100.0,
501                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
502                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
503                                    ".*-clock");
504     }
505 
506     @Test
507     void testClockOffsetAngularDerivativesOctantMxMyPz() {
508         double toleranceRotationValue              =  1.6e-15;
509         double toleranceRotationDerivative         =  4.9e-15;
510         double toleranceRotationRateValue          =  1.1e-19;
511         double toleranceRotationRateDerivative     =  6.3e-19;
512         doTestAngularDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 100.0,
513                                  toleranceRotationValue,     toleranceRotationDerivative,
514                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
515                                  ".*-clock");
516     }
517 
518     @Test
519     void testClockOffsetCartesianDerivativesOctantMxMyMz() {
520         double relativeTolerancePositionValue      =  1.5e-15;
521         double relativeTolerancePositionDerivative =  1.7e-10;
522         double relativeToleranceVelocityValue      =  2.9e-15;
523         double relativeToleranceVelocityDerivative =  1.9e-10;
524         doTestCartesianDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 100.0,
525                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
526                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
527                                    ".*-clock");
528     }
529 
530     @Test
531     void testClockOffsetAngularDerivativesOctantMxMyMz() {
532         double toleranceRotationValue              =  1.6e-15;
533         double toleranceRotationDerivative         =  4.9e-15;
534         double toleranceRotationRateValue          =  1.1e-19;
535         double toleranceRotationRateDerivative     =  6.3e-19;
536         doTestAngularDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 100.0,
537                                  toleranceRotationValue,     toleranceRotationDerivative,
538                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
539                                  ".*-clock");
540     }
541 
542     @Test
543     void testClockOffsetCartesianDerivativesNearPole() {
544         double relativeTolerancePositionValue      =  1.2e-15;
545         double relativeTolerancePositionDerivative =  1.6e-04;
546         double relativeToleranceVelocityValue      =  1.0e-13;
547         double relativeToleranceVelocityDerivative =  4.3e-09;
548         doTestCartesianDerivatives(FastMath.toRadians(89.99995), FastMath.toRadians(90), 1200.0, 100.0,
549                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
550                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
551                                    ".*-clock");
552     }
553 
554     @Test
555     void testClockOffsetAngularDerivativesNearPole() {
556         double toleranceRotationValue              =  1.5e-15;
557         double toleranceRotationDerivative         =  5.7e-15;
558         double toleranceRotationRateValue          =  1.1e-19;
559         double toleranceRotationRateDerivative     =  6.3e-19;
560         doTestAngularDerivatives(FastMath.toRadians(89.99995), FastMath.toRadians(90), 1200.0, 100.0,
561                                  toleranceRotationValue,     toleranceRotationDerivative,
562                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
563                                  ".*-clock");
564     }
565 
566     @Test
567     void testStationOffsetCartesianDerivativesOctantPxPyPz() {
568         double relativeTolerancePositionValue      =  2.3e-15;
569         double relativeTolerancePositionDerivative =  1.1e-10;
570         double relativeToleranceVelocityValue      =  3.3e-15;
571         double relativeToleranceVelocityDerivative =  1.2e-10;
572         doTestCartesianDerivatives(FastMath.toRadians(35), FastMath.toRadians(20), 1200.0, 100.0,
573                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
574                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
575                                    ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
576     }
577 
578     @Test
579     void testStationOffsetAngularDerivativesOctantPxPyPz() {
580         double toleranceRotationValue              =  1.9e-15;
581         double toleranceRotationDerivative         =  3.1e-18;
582         double toleranceRotationRateValue          =  1.5e-19;
583         double toleranceRotationRateDerivative     =  Precision.SAFE_MIN;
584         doTestAngularDerivatives(FastMath.toRadians(35), FastMath.toRadians(20), 1200.0, 100.0,
585                                  toleranceRotationValue,     toleranceRotationDerivative,
586                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
587                                  ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
588     }
589 
590     @Test
591     void testStationOffsetCartesianDerivativesOctantPxPyMz() {
592         double relativeTolerancePositionValue      =  1.4e-15;
593         double relativeTolerancePositionDerivative =  7.3e-11;
594         double relativeToleranceVelocityValue      =  2.8e-15;
595         double relativeToleranceVelocityDerivative =  1.3e-10;
596         doTestCartesianDerivatives(FastMath.toRadians(-35), FastMath.toRadians(20), 1200.0, 100.0,
597                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
598                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
599                                    ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
600     }
601 
602     @Test
603     void testStationOffsetAngularDerivativesOctantPxPyMz() {
604         double toleranceRotationValue            =  1.7e-15;
605         double toleranceRotationDerivative       =  3.6e-18;
606         double toleranceRotationRateValue        =  1.5e-19;
607         double toleranceRotationRateDerivative   =  Precision.SAFE_MIN;
608         doTestAngularDerivatives(FastMath.toRadians(-35), FastMath.toRadians(20), 1200.0, 100.0,
609                                  toleranceRotationValue,     toleranceRotationDerivative,
610                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
611                                  ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
612     }
613 
614     @Test
615     void testStationOffsetCartesianDerivativesOctantPxMyPz() {
616         double relativeTolerancePositionValue      =  1.7e-15;
617         double relativeTolerancePositionDerivative =  9.0e-11;
618         double relativeToleranceVelocityValue      =  2.9e-15;
619         double relativeToleranceVelocityDerivative =  8.8e-11;
620         doTestCartesianDerivatives(FastMath.toRadians(35), FastMath.toRadians(-20), 1200.0, 100.0,
621                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
622                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
623                                    ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
624     }
625 
626     @Test
627     void testStationOffsetAngularDerivativesOctantPxMyPz() {
628         double toleranceRotationValue              =  1.6e-15;
629         double toleranceRotationDerivative         =  3.6e-18;
630         double toleranceRotationRateValue          =  1.5e-19;
631         double toleranceRotationRateDerivative     =  Precision.SAFE_MIN;
632         doTestAngularDerivatives(FastMath.toRadians(35), FastMath.toRadians(-20), 1200.0, 100.0,
633                                  toleranceRotationValue,     toleranceRotationDerivative,
634                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
635                                  ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
636     }
637 
638     @Test
639     void testStationOffsetCartesianDerivativesOctantPxMyMz() {
640         double relativeTolerancePositionValue      =  1.5e-15;
641         double relativeTolerancePositionDerivative =  4.2e-11;
642         double relativeToleranceVelocityValue      =  2.7e-15;
643         double relativeToleranceVelocityDerivative =  6.8e-11;
644         doTestCartesianDerivatives(FastMath.toRadians(-35), FastMath.toRadians(-20), 1200.0, 100.0,
645                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
646                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
647                                    ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
648     }
649 
650     @Test
651     void testStationOffsetAngularDerivativesOctantPxMyMz() {
652         double toleranceRotationValue            =  1.6e-15;
653         double toleranceRotationDerivative       =  2.6e-18;
654         double toleranceRotationRateValue        =  1.5e-19;
655         double toleranceRotationRateDerivative   =  Precision.SAFE_MIN;
656         doTestAngularDerivatives(FastMath.toRadians(-35), FastMath.toRadians(-20), 1200.0, 100.0,
657                                  toleranceRotationValue,     toleranceRotationDerivative,
658                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
659                                  ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
660     }
661 
662     @Test
663     void testStationOffsetCartesianDerivativesOctantMxPyPz() {
664         double relativeTolerancePositionValue      =  1.6e-15;
665         double relativeTolerancePositionDerivative =  9.9e-11;
666         double relativeToleranceVelocityValue      =  2.6e-15;
667         double relativeToleranceVelocityDerivative =  9.5e-11;
668         doTestCartesianDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 100.0,
669                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
670                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
671                                    ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
672     }
673 
674     @Test
675     void testStationOffsetAngularDerivativesOctantMxPyPz() {
676         double toleranceRotationValue            =  1.5e-15;
677         double toleranceRotationDerivative       =  3.1e-18;
678         double toleranceRotationRateValue        =  1.5e-19;
679         double toleranceRotationRateDerivative   =  Precision.SAFE_MIN;
680         doTestAngularDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 100.0,
681                                  toleranceRotationValue,     toleranceRotationDerivative,
682                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
683                                  ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
684     }
685 
686     @Test
687     void testStationOffsetCartesianDerivativesOctantMxPyMz() {
688         double relativeTolerancePositionValue      =  1.6e-15;
689         double relativeTolerancePositionDerivative =  9.9e-11;
690         double relativeToleranceVelocityValue      =  2.6e-15;
691         double relativeToleranceVelocityDerivative =  9.5e-11;
692         doTestCartesianDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 100.0,
693                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
694                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
695                                    ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
696     }
697 
698     @Test
699     void testStationOffsetAngularDerivativesOctantMxPyMz() {
700         double toleranceRotationValue            =  1.5e-15;
701         double toleranceRotationDerivative       =  3.1e-18;
702         double toleranceRotationRateValue        =  1.5e-19;
703         double toleranceRotationRateDerivative   =  Precision.SAFE_MIN;
704         doTestAngularDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 100.0,
705                                  toleranceRotationValue,     toleranceRotationDerivative,
706                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
707                                  ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
708     }
709 
710     @Test
711     void testStationOffsetCartesianDerivativesOctantMxMyPz() {
712         double relativeTolerancePositionValue      =  1.5e-15;
713         double relativeTolerancePositionDerivative =  6.2e-11;
714         double relativeToleranceVelocityValue      =  3.2e-15;
715         double relativeToleranceVelocityDerivative =  1.1e-10;
716         doTestCartesianDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 100.0,
717                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
718                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
719                                    ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
720     }
721 
722     @Test
723     void testStationOffsetAngularDerivativesOctantMxMyPz() {
724         double toleranceRotationValue            =  1.6e-15;
725         double toleranceRotationDerivative       =  3.4e-18;
726         double toleranceRotationRateValue        =  1.5e-19;
727         double toleranceRotationRateDerivative   =  Precision.SAFE_MIN;
728         doTestAngularDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 100.0,
729                                  toleranceRotationValue,     toleranceRotationDerivative,
730                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
731                                  ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
732     }
733 
734     @Test
735     void testStationOffsetCartesianDerivativesOctantMxMyMz() {
736         double relativeTolerancePositionValue      =  1.5e-15;
737         double relativeTolerancePositionDerivative =  6.2e-11;
738         double relativeToleranceVelocityValue      =  3.2e-15;
739         double relativeToleranceVelocityDerivative =  1.1e-10;
740         doTestCartesianDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 100.0,
741                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
742                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
743                                    ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
744     }
745 
746     @Test
747     void testStationOffsetAngularDerivativesOctantMxMyMz() {
748         double toleranceRotationValue            =  1.6e-15;
749         double toleranceRotationDerivative       =  3.4e-18;
750         double toleranceRotationRateValue        =  1.5e-19;
751         double toleranceRotationRateDerivative   =  Precision.SAFE_MIN;
752         doTestAngularDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 100.0,
753                                  toleranceRotationValue,     toleranceRotationDerivative,
754                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
755                                  ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
756     }
757 
758     @Test
759     void testStationOffsetCartesianDerivativesNearPole() {
760         double relativeTolerancePositionValue      =  2.4e-15;
761         double relativeTolerancePositionDerivative =  2e-9;
762         double relativeToleranceVelocityValue      =  7.5e-14;
763         double relativeToleranceVelocityDerivative =  3.9e-10;
764         doTestCartesianDerivatives(FastMath.toRadians(89.99995), FastMath.toRadians(90), 1200.0, 100.0,
765                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
766                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
767                                    ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
768     }
769 
770     @Test
771     void testStationOffsetAngularDerivativesNearPole() {
772         double toleranceRotationValue              =  3.9e-15;
773         double toleranceRotationDerivative         =  8.0e-02; // near pole, the East and North directions are singular
774         double toleranceRotationRateValue          =  1.5e-19;
775         double toleranceRotationRateDerivative     =  Precision.SAFE_MIN;
776         doTestAngularDerivatives(FastMath.toRadians(89.99995), FastMath.toRadians(90), 1200.0, 100.0,
777                                  toleranceRotationValue,     toleranceRotationDerivative,
778                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
779                                  ".*-offset-East", ".*-offset-North", ".*-offset-Zenith");
780     }
781 
782     @Test
783     void testPolarMotionCartesianDerivativesOctantPxPyPz() {
784         double relativeTolerancePositionValue      =  2.3e-15;
785         double relativeTolerancePositionDerivative =  7.3e-09;
786         double relativeToleranceVelocityValue      =  3.3e-15;
787         double relativeToleranceVelocityDerivative =  4.8e-09;
788         doTestCartesianDerivatives(FastMath.toRadians(35), FastMath.toRadians(20), 1200.0, 1.0,
789                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
790                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
791                                    "polar-offset-X", "polar-drift-X",
792                                    "polar-offset-Y", "polar-drift-Y");
793     }
794 
795     @Test
796     void testPolarMotionAngularDerivativesOctantPxPyPz() {
797         double toleranceRotationValue              =  1.9e-15;
798         double toleranceRotationDerivative         =  6.6e-09;
799         double toleranceRotationRateValue          =  1.5e-19;
800         double toleranceRotationRateDerivative     =  2.2e-13;
801         doTestAngularDerivatives(FastMath.toRadians(35), FastMath.toRadians(20), 1200.0, 0.2,
802                                  toleranceRotationValue,     toleranceRotationDerivative,
803                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
804                                  "polar-offset-X", "polar-drift-X",
805                                  "polar-offset-Y", "polar-drift-Y");
806     }
807 
808     @Test
809     void testPolarMotionCartesianDerivativesOctantPxPyMz() {
810         double relativeTolerancePositionValue      =  1.4e-15;
811         double relativeTolerancePositionDerivative =  4.1e-09;
812         double relativeToleranceVelocityValue      =  2.8e-15;
813         double relativeToleranceVelocityDerivative =  4.7e-09;
814         doTestCartesianDerivatives(FastMath.toRadians(-35), FastMath.toRadians(20), 1200.0, 1.0,
815                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
816                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
817                                    "polar-offset-X", "polar-drift-X",
818                                    "polar-offset-Y", "polar-drift-Y");
819     }
820 
821     @Test
822     void testPolarMotionAngularDerivativesOctantPxPyMz() {
823         double toleranceRotationValue              =  1.7e-15;
824         double toleranceRotationDerivative         =  6.9e-09;
825         double toleranceRotationRateValue          =  1.5e-19;
826         double toleranceRotationRateDerivative     =  2.2e-13;
827         doTestAngularDerivatives(FastMath.toRadians(-35), FastMath.toRadians(20), 1200.0, 0.2,
828                                  toleranceRotationValue,     toleranceRotationDerivative,
829                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
830                                  "polar-offset-X", "polar-drift-X",
831                                  "polar-offset-Y", "polar-drift-Y");
832     }
833 
834     @Test
835     void testPolarMotionCartesianDerivativesOctantPxMyPz() {
836         double relativeTolerancePositionValue      =  1.7e-15;
837         double relativeTolerancePositionDerivative =  7.8e-09;
838         double relativeToleranceVelocityValue      =  2.9e-15;
839         double relativeToleranceVelocityDerivative =  4.4e-09;
840         doTestCartesianDerivatives(FastMath.toRadians(35), FastMath.toRadians(-20), 1200.0, 1.0,
841                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
842                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
843                                    "polar-offset-X", "polar-drift-X",
844                                    "polar-offset-Y", "polar-drift-Y");
845     }
846 
847     @Test
848     void testPolarMotionAngularDerivativesOctantPxMyPz() {
849         double toleranceRotationValue              =  1.6e-15;
850         double toleranceRotationDerivative         =  7.5e-09;
851         double toleranceRotationRateValue          =  1.5e-19;
852         double toleranceRotationRateDerivative     =  2.2e-13;
853         doTestAngularDerivatives(FastMath.toRadians(35), FastMath.toRadians(-20), 1200.0, 0.2,
854                                  toleranceRotationValue,     toleranceRotationDerivative,
855                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
856                                  "polar-offset-X", "polar-drift-X",
857                                  "polar-offset-Y", "polar-drift-Y");
858     }
859 
860     @Test
861     void testPolarMotionCartesianDerivativesOctantPxMyMz() {
862         double relativeTolerancePositionValue      =  1.5e-15;
863         double relativeTolerancePositionDerivative =  4.1e-09;
864         double relativeToleranceVelocityValue      =  2.7e-15;
865         double relativeToleranceVelocityDerivative =  4.4e-09;
866         doTestCartesianDerivatives(FastMath.toRadians(-35), FastMath.toRadians(-20), 1200.0, 1.0,
867                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
868                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
869                                    "polar-offset-X", "polar-drift-X",
870                                    "polar-offset-Y", "polar-drift-Y");
871     }
872 
873     @Test
874     void testPolarMotionAngularDerivativesOctantPxMyMz() {
875         double toleranceRotationValue            =  1.6e-15;
876         double toleranceRotationDerivative       =  7.3e-09;
877         double toleranceRotationRateValue        =  1.5e-19;
878         double toleranceRotationRateDerivative   =  2.2e-13;
879         doTestAngularDerivatives(FastMath.toRadians(-35), FastMath.toRadians(-20), 1200.0, 0.2,
880                                  toleranceRotationValue,     toleranceRotationDerivative,
881                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
882                                  "polar-offset-X", "polar-drift-X",
883                                  "polar-offset-Y", "polar-drift-Y");
884     }
885 
886     @Test
887     void testPolarMotionCartesianDerivativesOctantMxPyPz() {
888         double relativeTolerancePositionValue      =  1.6e-15;
889         double relativeTolerancePositionDerivative =  8.4e-09;
890         double relativeToleranceVelocityValue      =  2.6e-15;
891         double relativeToleranceVelocityDerivative =  5.0e-09;
892         doTestCartesianDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 1.0,
893                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
894                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
895                                    "polar-offset-X", "polar-drift-X",
896                                    "polar-offset-Y", "polar-drift-Y");
897     }
898 
899     @Test
900     void testPolarMotionAngularDerivativesOctantMxPyPz() {
901         double toleranceRotationValue              =  1.4e-15;
902         double toleranceRotationDerivative         =  6.3e-09;
903         double toleranceRotationRateValue          =  1.5e-19;
904         double toleranceRotationRateDerivative     =  2.2e-13;
905         doTestAngularDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 0.2,
906                                  toleranceRotationValue,     toleranceRotationDerivative,
907                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
908                                  "polar-offset-X", "polar-drift-X",
909                                  "polar-offset-Y", "polar-drift-Y");
910     }
911 
912     @Test
913     void testPolarMotionCartesianDerivativesOctantMxPyMz() {
914         double relativeTolerancePositionValue      =  1.6e-15;
915         double relativeTolerancePositionDerivative =  8.4e-09;
916         double relativeToleranceVelocityValue      =  2.6e-15;
917         double relativeToleranceVelocityDerivative =  5.0e-09;
918         doTestCartesianDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 1.0,
919                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
920                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
921                                    "polar-offset-X", "polar-drift-X",
922                                    "polar-offset-Y", "polar-drift-Y");
923     }
924 
925     @Test
926     void testPolarMotionAngularDerivativesOctantMxPyMz() {
927         double toleranceRotationValue              =  1.4e-15;
928         double toleranceRotationDerivative         =  6.3e-09;
929         double toleranceRotationRateValue          =  1.5e-19;
930         double toleranceRotationRateDerivative     =  2.2e-13;
931         doTestAngularDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 0.2,
932                                  toleranceRotationValue,     toleranceRotationDerivative,
933                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
934                                  "polar-offset-X", "polar-drift-X",
935                                  "polar-offset-Y", "polar-drift-Y");
936     }
937 
938     @Test
939     void testPolarMotionCartesianDerivativesOctantMxMyPz() {
940         double relativeTolerancePositionValue      =  1.5e-15;
941         double relativeTolerancePositionDerivative =  5.0e-09;
942         double relativeToleranceVelocityValue      =  3.2e-15;
943         double relativeToleranceVelocityDerivative =  5.3e-09;
944         doTestCartesianDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 1.0,
945                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
946                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
947                                    "polar-offset-X", "polar-drift-X",
948                                    "polar-offset-Y", "polar-drift-Y");
949     }
950 
951     @Test
952     void testPolarMotionAngularDerivativesOctantMxMyPz() {
953         double toleranceRotationValue              =  1.6e-15;
954         double toleranceRotationDerivative         =  7.7e-09;
955         double toleranceRotationRateValue          =  1.5e-19;
956         double toleranceRotationRateDerivative     =  2.2e-13;
957         doTestAngularDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 0.2,
958                                  toleranceRotationValue,     toleranceRotationDerivative,
959                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
960                                  "polar-offset-X", "polar-drift-X",
961                                  "polar-offset-Y", "polar-drift-Y");
962     }
963 
964     @Test
965     void testPolarMotionCartesianDerivativesOctantMxMyMz() {
966         double relativeTolerancePositionValue      =  1.5e-15;
967         double relativeTolerancePositionDerivative =  5.0e-09;
968         double relativeToleranceVelocityValue      =  3.2e-15;
969         double relativeToleranceVelocityDerivative =  5.3e-09;
970         doTestCartesianDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 1.0,
971                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
972                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
973                                    "polar-offset-X", "polar-drift-X",
974                                    "polar-offset-Y", "polar-drift-Y");
975     }
976 
977     @Test
978     void testPolarMotionAngularDerivativesOctantMxMyMz() {
979         double toleranceRotationValue              =  1.6e-15;
980         double toleranceRotationDerivative         =  7.7e-09;
981         double toleranceRotationRateValue          =  1.5e-19;
982         double toleranceRotationRateDerivative     =  2.2e-13;
983         doTestAngularDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 0.2,
984                                  toleranceRotationValue,     toleranceRotationDerivative,
985                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
986                                  "polar-offset-X", "polar-drift-X",
987                                  "polar-offset-Y", "polar-drift-Y");
988     }
989 
990     @Test
991     void testPolarMotionCartesianDerivativesNearPole() {
992         double relativeTolerancePositionValue      =  1.2e-15;
993         double relativeTolerancePositionDerivative =  5.7e-09;
994         double relativeToleranceVelocityValue      =  9.7e-13;
995         double relativeToleranceVelocityDerivative =  1.2e-09;
996         doTestCartesianDerivatives(FastMath.toRadians(89.99995), FastMath.toRadians(90), 1200.0, 1.0,
997                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
998                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
999                                    "polar-offset-X", "polar-drift-X",
1000                                    "polar-offset-Y", "polar-drift-Y");
1001     }
1002 
1003     @Test
1004     void testPolarMotionAngularDerivativesNearPole() {
1005         double toleranceRotationValue              =  1.5e-15;
1006         double toleranceRotationDerivative         =  6.5e-09;
1007         double toleranceRotationRateValue          =  1.5e-19;
1008         double toleranceRotationRateDerivative     =  2.2e-13;
1009         doTestAngularDerivatives(FastMath.toRadians(89.99995), FastMath.toRadians(90), 1200.0, 0.2,
1010                                  toleranceRotationValue,     toleranceRotationDerivative,
1011                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
1012                                  "polar-offset-X", "polar-drift-X",
1013                                  "polar-offset-Y", "polar-drift-Y");
1014     }
1015 
1016     @Test
1017     void testPrimeMeridianCartesianDerivativesOctantPxPyPz() {
1018         double relativeTolerancePositionValue      =  2.3e-15;
1019         double relativeTolerancePositionDerivative =  5.7e-09;
1020         double relativeToleranceVelocityValue      =  3.3e-15;
1021         double relativeToleranceVelocityDerivative =  2.7e-09;
1022         doTestCartesianDerivatives(FastMath.toRadians(35), FastMath.toRadians(20), 1200.0, 1.0,
1023                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
1024                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
1025                                    "prime-meridian-offset", "prime-meridian-drift");
1026     }
1027 
1028     @Test
1029     void testPrimeMeridianAngularDerivativesOctantPxPyPz() {
1030         double toleranceRotationValue              =  1.9e-15;
1031         double toleranceRotationDerivative         =  4.9e-09;
1032         double toleranceRotationRateValue          =  1.5e-19;
1033         double toleranceRotationRateDerivative     =  2.0e-13;
1034         doTestAngularDerivatives(FastMath.toRadians(35), FastMath.toRadians(20), 1200.0, 0.2,
1035                                  toleranceRotationValue,     toleranceRotationDerivative,
1036                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
1037                                  "prime-meridian-offset", "prime-meridian-drift");
1038     }
1039 
1040     @Test
1041     void testPrimeMeridianCartesianDerivativesOctantPxPyMz() {
1042         double relativeTolerancePositionValue      =  1.4e-15;
1043         double relativeTolerancePositionDerivative =  3.1e-09;
1044         double relativeToleranceVelocityValue      =  2.8e-15;
1045         double relativeToleranceVelocityDerivative =  3.3e-09;
1046         doTestCartesianDerivatives(FastMath.toRadians(-35), FastMath.toRadians(20), 1200.0, 1.0,
1047                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
1048                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
1049                                    "prime-meridian-offset", "prime-meridian-drift");
1050     }
1051 
1052     @Test
1053     void testPrimeMeridianAngularDerivativesOctantPxPyMz() {
1054         double toleranceRotationValue            =  1.7e-15;
1055         double toleranceRotationDerivative       =  5.7e-09;
1056         double toleranceRotationRateValue        =  1.5e-19;
1057         double toleranceRotationRateDerivative   =  2.0e-13;
1058         doTestAngularDerivatives(FastMath.toRadians(-35), FastMath.toRadians(20), 1200.0, 0.2,
1059                                  toleranceRotationValue,     toleranceRotationDerivative,
1060                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
1061                                  "prime-meridian-offset", "prime-meridian-drift");
1062     }
1063 
1064     @Test
1065     void testPrimeMeridianCartesianDerivativesOctantPxMyPz() {
1066         double relativeTolerancePositionValue      =  1.7e-15;
1067         double relativeTolerancePositionDerivative =  5.3e-09;
1068         double relativeToleranceVelocityValue      =  2.9e-15;
1069         double relativeToleranceVelocityDerivative =  3.4e-09;
1070         doTestCartesianDerivatives(FastMath.toRadians(35), FastMath.toRadians(-20), 1200.0, 1.0,
1071                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
1072                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
1073                                    "prime-meridian-offset", "prime-meridian-drift");
1074     }
1075 
1076     @Test
1077     void testPrimeMeridianAngularDerivativesOctantPxMyPz() {
1078         double toleranceRotationValue              =  1.6e-15;
1079         double toleranceRotationDerivative         =  6.4e-09;
1080         double toleranceRotationRateValue          =  1.5e-19;
1081         double toleranceRotationRateDerivative     =  2.0e-13;
1082         doTestAngularDerivatives(FastMath.toRadians(35), FastMath.toRadians(-20), 1200.0, 0.2,
1083                                  toleranceRotationValue,     toleranceRotationDerivative,
1084                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
1085                                  "prime-meridian-offset", "prime-meridian-drift");
1086     }
1087 
1088     @Test
1089     void testPrimeMeridianCartesianDerivativesOctantPxMyMz() {
1090         double relativeTolerancePositionValue      =  1.5e-15;
1091         double relativeTolerancePositionDerivative =  3.1e-09;
1092         double relativeToleranceVelocityValue      =  2.7e-15;
1093         double relativeToleranceVelocityDerivative =  3.0e-09;
1094         doTestCartesianDerivatives(FastMath.toRadians(-35), FastMath.toRadians(-20), 1200.0, 1.0,
1095                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
1096                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
1097                                    "prime-meridian-offset", "prime-meridian-drift");
1098     }
1099 
1100     @Test
1101     void testPrimeMeridianAngularDerivativesOctantPxMyMz() {
1102         double toleranceRotationValue              =  1.5e-15;
1103         double toleranceRotationDerivative         =  6.3e-09;
1104         double toleranceRotationRateValue          =  1.5e-19;
1105         double toleranceRotationRateDerivative     =  2.0e-13;
1106         doTestAngularDerivatives(FastMath.toRadians(-35), FastMath.toRadians(-20), 1200.0, 0.2,
1107                                  toleranceRotationValue,     toleranceRotationDerivative,
1108                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
1109                                  "prime-meridian-offset", "prime-meridian-drift");
1110     }
1111 
1112     @Test
1113     void testPrimeMeridianCartesianDerivativesOctantMxPyPz() {
1114         double relativeTolerancePositionValue      =  1.6e-15;
1115         double relativeTolerancePositionDerivative =  5.5e-09;
1116         double relativeToleranceVelocityValue      =  2.6e-15;
1117         double relativeToleranceVelocityDerivative =  2.8e-09;
1118         doTestCartesianDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 1.0,
1119                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
1120                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
1121                                    "prime-meridian-offset", "prime-meridian-drift");
1122     }
1123 
1124     @Test
1125     void testPrimeMeridianAngularDerivativesOctantMxPyPz() {
1126         double toleranceRotationValue            =  1.5e-15;
1127         double toleranceRotationDerivative       =  6.8e-09;
1128         double toleranceRotationRateValue        =  1.5e-19;
1129         double toleranceRotationRateDerivative   =  2.0e-13;
1130         doTestAngularDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 0.2,
1131                                  toleranceRotationValue,     toleranceRotationDerivative,
1132                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
1133                                  "prime-meridian-offset", "prime-meridian-drift");
1134     }
1135 
1136     @Test
1137     void testPrimeMeridianCartesianDerivativesOctantMxPyMz() {
1138         double relativeTolerancePositionValue      =  1.6e-15;
1139         double relativeTolerancePositionDerivative =  5.5e-09;
1140         double relativeToleranceVelocityValue      =  2.6e-15;
1141         double relativeToleranceVelocityDerivative =  2.8e-09;
1142         doTestCartesianDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 1.0,
1143                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
1144                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
1145                                    "prime-meridian-offset", "prime-meridian-drift");
1146     }
1147 
1148     @Test
1149     void testPrimeMeridianAngularDerivativesOctantMxPyMz() {
1150         double toleranceRotationValue            =  1.5e-15;
1151         double toleranceRotationDerivative       =  6.8e-09;
1152         double toleranceRotationRateValue        =  1.5e-19;
1153         double toleranceRotationRateDerivative   =  2.0e-13;
1154         doTestAngularDerivatives(FastMath.toRadians(150), FastMath.toRadians(20), 1200.0, 0.2,
1155                                  toleranceRotationValue,     toleranceRotationDerivative,
1156                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
1157                                  "prime-meridian-offset", "prime-meridian-drift");
1158     }
1159 
1160     @Test
1161     void testPrimeMeridianCartesianDerivativesOctantMxMyPz() {
1162         double relativeTolerancePositionValue      =  1.5e-15;
1163         double relativeTolerancePositionDerivative =  3.1e-09;
1164         double relativeToleranceVelocityValue      =  3.2e-15;
1165         double relativeToleranceVelocityDerivative =  3.2e-09;
1166         doTestCartesianDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 1.0,
1167                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
1168                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
1169                                    "prime-meridian-offset", "prime-meridian-drift");
1170     }
1171 
1172     @Test
1173     void testPrimeMeridianAngularDerivativesOctantMxMyPz() {
1174         double toleranceRotationValue              =  1.6e-15;
1175         double toleranceRotationDerivative         =  6.2e-09;
1176         double toleranceRotationRateValue          =  1.5e-19;
1177         double toleranceRotationRateDerivative     =  2.0e-13;
1178         doTestAngularDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 0.2,
1179                                  toleranceRotationValue,     toleranceRotationDerivative,
1180                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
1181                                  "prime-meridian-offset", "prime-meridian-drift");
1182     }
1183 
1184     @Test
1185     void testPrimeMeridianCartesianDerivativesOctantMxMyMz() {
1186         double relativeTolerancePositionValue      =  1.5e-15;
1187         double relativeTolerancePositionDerivative =  3.1e-09;
1188         double relativeToleranceVelocityValue      =  3.2e-15;
1189         double relativeToleranceVelocityDerivative =  3.2e-09;
1190         doTestCartesianDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 1.0,
1191                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
1192                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
1193                                    "prime-meridian-offset", "prime-meridian-drift");
1194     }
1195 
1196     @Test
1197     void testPrimeMeridianAngularDerivativesOctantMxMyMz() {
1198         double toleranceRotationValue              =  1.6e-15;
1199         double toleranceRotationDerivative         =  6.2e-09;
1200         double toleranceRotationRateValue          =  1.5e-19;
1201         double toleranceRotationRateDerivative     =  2.0e-13;
1202         doTestAngularDerivatives(FastMath.toRadians(-150), FastMath.toRadians(-20), 1200.0, 0.2,
1203                                  toleranceRotationValue,     toleranceRotationDerivative,
1204                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
1205                                  "prime-meridian-offset", "prime-meridian-drift");
1206     }
1207 
1208     @Test
1209     void testPrimeMeridianCartesianDerivativesNearPole() {
1210         double relativeTolerancePositionValue      =  1.2e-15;
1211         double relativeTolerancePositionDerivative =  1.6e-03;
1212         double relativeToleranceVelocityValue      =  5.8e-14;
1213         double relativeToleranceVelocityDerivative =  2.6e-08;
1214         doTestCartesianDerivatives(FastMath.toRadians(89.99995), FastMath.toRadians(90), 1200.0, 1.0,
1215                                    relativeTolerancePositionValue, relativeTolerancePositionDerivative,
1216                                    relativeToleranceVelocityValue, relativeToleranceVelocityDerivative,
1217                                    "prime-meridian-offset", "prime-meridian-drift");
1218     }
1219 
1220     @Test
1221     void testPrimeMeridianAngularDerivativesNearPole() {
1222         double toleranceRotationValue              =  1.5e-15;
1223         double toleranceRotationDerivative         =  7.3e-09;
1224         double toleranceRotationRateValue          =  1.5e-19;
1225         double toleranceRotationRateDerivative     =  2.0e-13;
1226         doTestAngularDerivatives(FastMath.toRadians(89.99995), FastMath.toRadians(90), 1200.0, 0.2,
1227                                  toleranceRotationValue,     toleranceRotationDerivative,
1228                                  toleranceRotationRateValue, toleranceRotationRateDerivative,
1229                                  "prime-meridian-offset", "prime-meridian-drift");
1230     }
1231 
1232     @Test
1233     void testNoReferenceDateGradient() {
1234         Utils.setDataRoot("regular-data");
1235         final Frame eme2000 = FramesFactory.getEME2000();
1236         final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
1237         final OneAxisEllipsoid earth =
1238                         new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1239                                              Constants.WGS84_EARTH_FLATTENING,
1240                                              FramesFactory.getITRF(IERSConventions.IERS_2010, true));
1241         final GroundStation station = new GroundStation(new TopocentricFrame(earth,
1242                                                                              new GeodeticPoint(0.1, 0.2, 100),
1243                                                                              "dummy"));
1244         try {
1245             station.getOffsetToInertial(eme2000, date, false);
1246             Assertions.fail("an exception should have been thrown");
1247         } catch (OrekitException oe) {
1248             Assertions.assertEquals(OrekitMessages.NO_REFERENCE_DATE_FOR_PARAMETER, oe.getSpecifier());
1249             Assertions.assertEquals("prime-meridian-offset", oe.getParts()[0]);
1250         }
1251 
1252         try {
1253             int freeParameters = 9;
1254             Map<String, Integer> indices = new HashMap<>();
1255             for (final ParameterDriver driver : Arrays.asList(station.getPrimeMeridianOffsetDriver(),
1256                                                               station.getPrimeMeridianDriftDriver(),
1257                                                               station.getPolarOffsetXDriver(),
1258                                                               station.getPolarDriftXDriver(),
1259                                                               station.getPolarOffsetYDriver(),
1260                                                               station.getPolarDriftYDriver(),
1261                                                               station.getClockOffsetDriver(),
1262                                                               station.getEastOffsetDriver(),
1263                                                               station.getNorthOffsetDriver(),
1264                                                               station.getZenithOffsetDriver())) {
1265                 indices.put(driver.getNameSpan(date), indices.size());
1266             }
1267             station.getOffsetToInertial(eme2000, date, freeParameters, indices);
1268             Assertions.fail("an exception should have been thrown");
1269         } catch (OrekitException oe) {
1270             Assertions.assertEquals(OrekitMessages.NO_REFERENCE_DATE_FOR_PARAMETER, oe.getSpecifier());
1271             Assertions.assertEquals("prime-meridian-offset", oe.getParts()[0]);
1272         }
1273 
1274     }
1275 
1276     private void doTestCartesianDerivatives(double latitude, double longitude, double altitude, double stepFactor,
1277                                             double relativeTolerancePositionValue, double relativeTolerancePositionDerivative,
1278                                             double relativeToleranceVelocityValue, double relativeToleranceVelocityDerivative,
1279                                             String... parameterPattern) {
1280         Utils.setDataRoot("regular-data");
1281         final Frame eme2000 = FramesFactory.getEME2000();
1282         final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
1283         final AbsoluteDate date0 = date.shiftedBy(50000);
1284         final OneAxisEllipsoid earth =
1285                         new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1286                                              Constants.WGS84_EARTH_FLATTENING,
1287                                              FramesFactory.getITRF(IERSConventions.IERS_2010, true));
1288         final GroundStation station = new GroundStation(new TopocentricFrame(earth,
1289                                                                              new GeodeticPoint(latitude, longitude, altitude),
1290                                                                              "dummy"));
1291         final GradientField gradientField = GradientField.getField(parameterPattern.length);
1292         ParameterDriver[] selectedDrivers = new ParameterDriver[parameterPattern.length];
1293         UnivariateDifferentiableVectorFunction[] dFCartesian = new UnivariateDifferentiableVectorFunction[parameterPattern.length];
1294         final ParameterDriver[] allDrivers = selectAllDrivers(station);
1295         for (ParameterDriver driver : allDrivers) {
1296             driver.setReferenceDate(date0);
1297         }
1298         Map<String, Integer> indices = new HashMap<>();
1299         for (int k = 0; k < dFCartesian.length; ++k) {
1300             for (final ParameterDriver allDriver : allDrivers) {
1301                 if (allDriver.getName().matches(parameterPattern[k])) {
1302                     selectedDrivers[k] = allDriver;
1303                     dFCartesian[k] = differentiatedStationPV(station, eme2000, date, selectedDrivers[k], stepFactor);
1304                     indices.put(selectedDrivers[k].getNameSpan(date0), k);
1305                 }
1306             }
1307         }
1308 
1309         RandomGenerator generator = new Well19937a(0x084d58a19c498a54L);
1310 
1311         double maxPositionValueRelativeError      = 0;
1312         double maxPositionDerivativeRelativeError = 0;
1313         double maxVelocityValueRelativeError      = 0;
1314         double maxVelocityDerivativeRelativeError = 0;
1315         for (int i = 0; i < 1000; ++i) {
1316 
1317             // randomly change one parameter
1318             ParameterDriver changed = allDrivers[generator.nextInt(allDrivers.length)];
1319             changed.setNormalizedValue(2 * generator.nextDouble() - 1);
1320 
1321             // transform to check
1322             FieldTransform<Gradient> t = station.getOffsetToInertial(eme2000, date, parameterPattern.length, indices);
1323             FieldPVCoordinates<Gradient> pv = t.transformPVCoordinates(FieldPVCoordinates.getZero(gradientField));
1324             for (int k = 0; k < dFCartesian.length; ++k) {
1325 
1326                 // reference values and derivatives computed using finite differences
1327                 Gradient[] refCartesian = dFCartesian[k].value(Gradient.variable(1, 0, selectedDrivers[k].getValue()));
1328 
1329                 // position
1330                 final Vector3D refP = new Vector3D(refCartesian[0].getValue(),
1331                                                    refCartesian[1].getValue(),
1332                                                    refCartesian[2].getValue());
1333                 final Vector3D resP = new Vector3D(pv.getPosition().getX().getValue(),
1334                                                    pv.getPosition().getY().getValue(),
1335                                                    pv.getPosition().getZ().getValue());
1336                 maxPositionValueRelativeError =
1337                                 FastMath.max(maxPositionValueRelativeError, Vector3D.distance(refP, resP) / refP.getNorm());
1338                 final Vector3D refPD = new Vector3D(refCartesian[0].getPartialDerivative(0),
1339                                                     refCartesian[1].getPartialDerivative(0),
1340                                                     refCartesian[2].getPartialDerivative(0));
1341                 final Vector3D resPD = new Vector3D(pv.getPosition().getX().getPartialDerivative(k),
1342                                                     pv.getPosition().getY().getPartialDerivative(k),
1343                                                     pv.getPosition().getZ().getPartialDerivative(k));
1344                 maxPositionDerivativeRelativeError =
1345                                 FastMath.max(maxPositionDerivativeRelativeError, Vector3D.distance(refPD, resPD) / refPD.getNorm());
1346 
1347                 // velocity
1348                 final Vector3D refV = new Vector3D(refCartesian[3].getValue(),
1349                                                    refCartesian[4].getValue(),
1350                                                    refCartesian[5].getValue());
1351                 final Vector3D resV = new Vector3D(pv.getVelocity().getX().getValue(),
1352                                                    pv.getVelocity().getY().getValue(),
1353                                                    pv.getVelocity().getZ().getValue());
1354                 maxVelocityValueRelativeError =
1355                                 FastMath.max(maxVelocityValueRelativeError, Vector3D.distance(refV, resV) / refV.getNorm());
1356                 final Vector3D refVD = new Vector3D(refCartesian[3].getPartialDerivative(0),
1357                                                     refCartesian[4].getPartialDerivative(0),
1358                                                     refCartesian[5].getPartialDerivative(0));
1359                 final Vector3D resVD = new Vector3D(pv.getVelocity().getX().getPartialDerivative(k),
1360                                                     pv.getVelocity().getY().getPartialDerivative(k),
1361                                                     pv.getVelocity().getZ().getPartialDerivative(k));
1362                 maxVelocityDerivativeRelativeError =
1363                                 FastMath.max(maxVelocityDerivativeRelativeError, Vector3D.distance(refVD, resVD) / refVD.getNorm());
1364 
1365             }
1366 
1367         }
1368 
1369         if (maxPositionValueRelativeError      > relativeTolerancePositionValue      ||
1370             maxPositionDerivativeRelativeError > relativeTolerancePositionDerivative ||
1371             maxVelocityValueRelativeError      > relativeToleranceVelocityValue      ||
1372             maxVelocityDerivativeRelativeError > relativeToleranceVelocityDerivative) {
1373             print("relativeTolerancePositionValue",      maxPositionValueRelativeError);
1374             print("relativeTolerancePositionDerivative", maxPositionDerivativeRelativeError);
1375             print("relativeToleranceVelocityValue",      maxVelocityValueRelativeError);
1376             print("relativeToleranceVelocityDerivative", maxVelocityDerivativeRelativeError);
1377         }
1378         Assertions.assertEquals(0.0, maxPositionValueRelativeError,      relativeTolerancePositionValue);
1379         Assertions.assertEquals(0.0, maxPositionDerivativeRelativeError, relativeTolerancePositionDerivative);
1380         Assertions.assertEquals(0.0, maxVelocityValueRelativeError,      relativeToleranceVelocityValue);
1381         Assertions.assertEquals(0.0, maxVelocityDerivativeRelativeError, relativeToleranceVelocityDerivative);
1382 
1383     }
1384 
1385     private void doTestAngularDerivatives(double latitude, double longitude, double altitude, double stepFactor,
1386                                           double toleranceRotationValue,     double toleranceRotationDerivative,
1387                                           double toleranceRotationRateValue, double toleranceRotationRateDerivative,
1388                                           String... parameterPattern) {
1389         Utils.setDataRoot("regular-data");
1390         final Frame eme2000 = FramesFactory.getEME2000();
1391         final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
1392         final AbsoluteDate date0 = date.shiftedBy(50000);
1393         final OneAxisEllipsoid earth =
1394                         new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1395                                              Constants.WGS84_EARTH_FLATTENING,
1396                                              FramesFactory.getITRF(IERSConventions.IERS_2010, true));
1397         final GroundStation station = new GroundStation(new TopocentricFrame(earth,
1398                                                                              new GeodeticPoint(latitude, longitude, altitude),
1399                                                                              "dummy"));
1400         ParameterDriver[] selectedDrivers = new ParameterDriver[parameterPattern.length];
1401         UnivariateDifferentiableVectorFunction[] dFAngular   = new UnivariateDifferentiableVectorFunction[parameterPattern.length];
1402         final ParameterDriver[] allDrivers = selectAllDrivers(station);
1403         for (ParameterDriver driver : allDrivers) {
1404             driver.setReferenceDate(date0);
1405         }
1406         Map<String, Integer> indices = new HashMap<>();
1407         for (int k = 0; k < dFAngular.length; ++k) {
1408             for (final ParameterDriver allDriver : allDrivers) {
1409                 if (allDriver.getName().matches(parameterPattern[k])) {
1410                     selectedDrivers[k] = allDriver;
1411                     dFAngular[k] = differentiatedTransformAngular(station, eme2000, date, selectedDrivers[k], stepFactor);
1412                     indices.put(selectedDrivers[k].getNameSpan(date0), k);
1413                 }
1414             }
1415         }
1416         RandomGenerator generator = new Well19937a(0xa01a1d8fe5d80af7L);
1417 
1418         double maxRotationValueError          = 0;
1419         double maxRotationDerivativeError     = 0;
1420         double maxRotationRateValueError      = 0;
1421         double maxRotationRateDerivativeError = 0;
1422         for (int i = 0; i < 1000; ++i) {
1423 
1424             // randomly change one parameter
1425             ParameterDriver changed = allDrivers[generator.nextInt(allDrivers.length)];
1426             changed.setNormalizedValue(2 * generator.nextDouble() - 1);
1427 
1428             // transform to check
1429             FieldTransform<Gradient> t = station.getOffsetToInertial(eme2000, date, parameterPattern.length, indices);
1430             for (int k = 0; k < dFAngular.length; ++k) {
1431 
1432                 // reference values and derivatives computed using finite differences
1433                 Gradient[] refAngular = dFAngular[k].value(Gradient.variable(1, 0, selectedDrivers[k].getValue()));
1434 
1435                 // rotation
1436                 final Rotation refQ = new Rotation(refAngular[0].getValue(),
1437                                                    refAngular[1].getValue(),
1438                                                    refAngular[2].getValue(),
1439                                                    refAngular[3].getValue(),
1440                                                    true);
1441                 final Rotation resQ = t.getRotation().toRotation();
1442                 maxRotationValueError      = FastMath.max(maxRotationValueError, Rotation.distance(refQ, resQ));
1443                 double sign = FastMath.copySign(1.0,
1444                                                 refAngular[0].getValue() * t.getRotation().getQ0().getValue() +
1445                                                 refAngular[1].getValue() * t.getRotation().getQ1().getValue() +
1446                                                 refAngular[2].getValue() * t.getRotation().getQ2().getValue() +
1447                                                 refAngular[3].getValue() * t.getRotation().getQ3().getValue());
1448                 maxRotationDerivativeError = FastMath.max(maxRotationDerivativeError,
1449                                                           FastMath.abs(sign * refAngular[0].getPartialDerivative(0) -
1450                                                                        t.getRotation().getQ0().getPartialDerivative(k)));
1451                 maxRotationDerivativeError = FastMath.max(maxRotationDerivativeError,
1452                                                           FastMath.abs(sign * refAngular[1].getPartialDerivative(0) -
1453                                                                        t.getRotation().getQ1().getPartialDerivative(k)));
1454                 maxRotationDerivativeError = FastMath.max(maxRotationDerivativeError,
1455                                                           FastMath.abs(sign * refAngular[2].getPartialDerivative(0) -
1456                                                                        t.getRotation().getQ2().getPartialDerivative(k)));
1457                 maxRotationDerivativeError = FastMath.max(maxRotationDerivativeError,
1458                                                           FastMath.abs(sign * refAngular[3].getPartialDerivative(0) -
1459                                                                        t.getRotation().getQ3().getPartialDerivative(k)));
1460 
1461                 // rotation rate
1462                 final Vector3D refRate  = new Vector3D(refAngular[4].getValue(), refAngular[5].getValue(), refAngular[6].getValue());
1463                 final Vector3D resRate  = t.getRotationRate().toVector3D();
1464                 final Vector3D refRateD = new Vector3D(refAngular[4].getPartialDerivative(0),
1465                                                        refAngular[5].getPartialDerivative(0),
1466                                                        refAngular[6].getPartialDerivative(0));
1467                 final Vector3D resRateD = new Vector3D(t.getRotationRate().getX().getPartialDerivative(k),
1468                                                        t.getRotationRate().getY().getPartialDerivative(k),
1469                                                        t.getRotationRate().getZ().getPartialDerivative(k));
1470                 maxRotationRateValueError      = FastMath.max(maxRotationRateValueError, Vector3D.distance(refRate, resRate));
1471                 maxRotationRateDerivativeError = FastMath.max(maxRotationRateDerivativeError, Vector3D.distance(refRateD, resRateD));
1472 
1473             }
1474 
1475         }
1476 
1477         if (maxRotationValueError          > toleranceRotationValue          ||
1478             maxRotationDerivativeError     > toleranceRotationDerivative     ||
1479             maxRotationRateValueError      > toleranceRotationRateValue      ||
1480             maxRotationRateDerivativeError > toleranceRotationRateDerivative) {
1481             print("toleranceRotationValue",          maxRotationValueError);
1482             print("toleranceRotationDerivative",     maxRotationDerivativeError);
1483             print("toleranceRotationRateValue",      maxRotationRateValueError);
1484             print("toleranceRotationRateDerivative", maxRotationRateDerivativeError);
1485         }
1486         Assertions.assertEquals(0.0, maxRotationValueError,           toleranceRotationValue);
1487         Assertions.assertEquals(0.0, maxRotationDerivativeError,      toleranceRotationDerivative);
1488         Assertions.assertEquals(0.0, maxRotationRateValueError,       toleranceRotationRateValue);
1489         Assertions.assertEquals(0.0, maxRotationRateDerivativeError,  toleranceRotationRateDerivative);
1490 
1491     }
1492     private void print(String name, double v) {
1493         if (v < Precision.SAFE_MIN) {
1494             System.out.format(Locale.US, "            double %-35s =  Precision.SAFE_MIN;%n", name);
1495         } else {
1496             double s = FastMath.pow(10.0, 1 - FastMath.floor(FastMath.log(v) / FastMath.log(10.0)));
1497             System.out.format(Locale.US, "            double %-35s = %8.1e;%n",
1498                               name, FastMath.ceil(s * v) / s);
1499         }
1500     }
1501 
1502     private UnivariateDifferentiableVectorFunction differentiatedStationPV(final GroundStation station,
1503                                                                            final Frame eme2000,
1504                                                                            final AbsoluteDate date,
1505                                                                            final ParameterDriver driver,
1506                                                                            final double stepFactor) {
1507 
1508         final FiniteDifferencesDifferentiator differentiator =
1509                         new FiniteDifferencesDifferentiator(5, stepFactor * driver.getScale());
1510 
1511         return differentiator.differentiate((UnivariateVectorFunction) x -> {
1512             final double[] result = new double[6];
1513             try {
1514                 final double previouspI = driver.getValue(date);
1515                 driver.setValue(x, new AbsoluteDate());
1516                 Transform t = station.getOffsetToInertial(eme2000, date, false);
1517                 driver.setValue(previouspI, date);
1518                 PVCoordinates stationPV = t.transformPVCoordinates(PVCoordinates.ZERO);
1519                 result[ 0] = stationPV.getPosition().getX();
1520                 result[ 1] = stationPV.getPosition().getY();
1521                 result[ 2] = stationPV.getPosition().getZ();
1522                 result[ 3] = stationPV.getVelocity().getX();
1523                 result[ 4] = stationPV.getVelocity().getY();
1524                 result[ 5] = stationPV.getVelocity().getZ();
1525             } catch (OrekitException oe) {
1526                 Assertions.fail(oe.getLocalizedMessage());
1527             }
1528             return result;
1529         });
1530     }
1531 
1532     private UnivariateDifferentiableVectorFunction differentiatedTransformAngular(final GroundStation station,
1533                                                                                   final Frame eme2000,
1534                                                                                   final AbsoluteDate date,
1535                                                                                   final ParameterDriver driver,
1536                                                                                   final double stepFactor) {
1537 
1538         final FiniteDifferencesDifferentiator differentiator =
1539                         new FiniteDifferencesDifferentiator(5, stepFactor * driver.getScale());
1540 
1541         return differentiator.differentiate(new UnivariateVectorFunction() {
1542             private double previous0 = Double.NaN;
1543             private double previous1 = Double.NaN;
1544             private double previous2 = Double.NaN;
1545             private double previous3 = Double.NaN;
1546             @Override
1547             public double[] value(double x) {
1548                 final double[] result = new double[7];
1549                 try {
1550                     final double previouspI = driver.getValue(date);
1551                     driver.setValue(x, date);
1552                     Transform t = station.getOffsetToInertial(eme2000, date, false);
1553                     driver.setValue(previouspI, date);
1554                     final double sign;
1555                     if (Double.isNaN(previous0)) {
1556                         sign = +1;
1557                     } else {
1558                         sign = FastMath.copySign(1.0,
1559                                                  previous0 * t.getRotation().getQ0() +
1560                                                  previous1 * t.getRotation().getQ1() +
1561                                                  previous2 * t.getRotation().getQ2() +
1562                                                  previous3 * t.getRotation().getQ3());
1563                     }
1564                     previous0 = sign * t.getRotation().getQ0();
1565                     previous1 = sign * t.getRotation().getQ1();
1566                     previous2 = sign * t.getRotation().getQ2();
1567                     previous3 = sign * t.getRotation().getQ3();
1568                     result[0] = previous0;
1569                     result[1] = previous1;
1570                     result[2] = previous2;
1571                     result[3] = previous3;
1572                     result[4] = t.getRotationRate().getX();
1573                     result[5] = t.getRotationRate().getY();
1574                     result[6] = t.getRotationRate().getZ();
1575                 } catch (OrekitException oe) {
1576                     Assertions.fail(oe.getLocalizedMessage());
1577                 }
1578                 return result;
1579             }
1580         });
1581     }
1582 
1583     private ParameterDriver[] selectAllDrivers(final GroundStation station) {
1584         return new ParameterDriver[] {
1585             station.getPrimeMeridianOffsetDriver(),
1586             station.getPrimeMeridianDriftDriver(),
1587             station.getPolarOffsetXDriver(),
1588             station.getPolarDriftXDriver(),
1589             station.getPolarOffsetYDriver(),
1590             station.getPolarDriftYDriver(),
1591             station.getClockOffsetDriver(),
1592             station.getEastOffsetDriver(),
1593             station.getNorthOffsetDriver(),
1594             station.getZenithOffsetDriver()
1595         };
1596     }
1597 
1598 }
1599