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.frames;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.geometry.euclidean.threed.Rotation;
21  import org.hipparchus.geometry.euclidean.threed.RotationConvention;
22  import org.hipparchus.geometry.euclidean.threed.Vector3D;
23  import org.junit.jupiter.api.Assertions;
24  import org.junit.jupiter.api.BeforeEach;
25  import org.junit.jupiter.api.Test;
26  import org.orekit.Utils;
27  import org.orekit.data.DataContext;
28  import org.orekit.data.PolynomialNutation;
29  import org.orekit.time.AbsoluteDate;
30  import org.orekit.time.DateComponents;
31  import org.orekit.time.FieldAbsoluteDate;
32  import org.orekit.time.TimeComponents;
33  import org.orekit.time.TimeScalesFactory;
34  import org.orekit.utils.Constants;
35  import org.orekit.utils.IERSConventions;
36  import org.orekit.utils.PVCoordinates;
37  
38  
39  public class MODProviderTest {
40  
41      @Test
42      public void testEuler1976() {
43  
44          TransformProvider eulerBasedProvider = new TransformProvider() {
45              private final PolynomialNutation zetaA =
46                      new PolynomialNutation(0.0,
47                              2306.2181 * Constants.ARC_SECONDS_TO_RADIANS,
48                              0.30188   * Constants.ARC_SECONDS_TO_RADIANS,
49                              0.017998  * Constants.ARC_SECONDS_TO_RADIANS);
50              private final PolynomialNutation thetaA =
51                      new PolynomialNutation(0.0,
52                              2004.3109 * Constants.ARC_SECONDS_TO_RADIANS,
53                              -0.42665  * Constants.ARC_SECONDS_TO_RADIANS,
54                              -0.041833 * Constants.ARC_SECONDS_TO_RADIANS);
55              private final PolynomialNutation zA =
56                      new PolynomialNutation(0.0,
57                              2306.2181 * Constants.ARC_SECONDS_TO_RADIANS,
58                              1.09468   * Constants.ARC_SECONDS_TO_RADIANS,
59                              0.018203  * Constants.ARC_SECONDS_TO_RADIANS);
60  
61              public Transform getTransform(AbsoluteDate date) {
62                  final double tc = IERSConventions.IERS_1996.evaluateTC(date);
63                  final Rotation r1 = new Rotation(Vector3D.PLUS_K,  zA.value(tc), RotationConvention.VECTOR_OPERATOR);
64                  final Rotation r2 = new Rotation(Vector3D.PLUS_J, -thetaA.value(tc), RotationConvention.VECTOR_OPERATOR);
65                  final Rotation r3 = new Rotation(Vector3D.PLUS_K,  zetaA.value(tc), RotationConvention.VECTOR_OPERATOR);
66                  return new Transform(date, r1.compose(r2.compose(r3,
67                                                                   RotationConvention.VECTOR_OPERATOR),
68                                                        RotationConvention.VECTOR_OPERATOR));
69              }
70  
71              public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
72                  throw new UnsupportedOperationException("never called in this test");
73              }
74  
75          };
76  
77          MODProvider modProvider = new MODProvider(IERSConventions.IERS_1996,
78                  DataContext.getDefault().getTimeScales());
79  
80          for (double dt = -5 * Constants.JULIAN_YEAR; dt < 5 * Constants.JULIAN_YEAR; dt += 10 * Constants.JULIAN_DAY) {
81              AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(dt);
82              Transform t = new Transform(date,
83                                          modProvider.getTransform(date).getInverse(),
84                                          eulerBasedProvider.getTransform(date));
85              Assertions.assertEquals(0, t.getRotation().getAngle(), 1.01e-11);
86          }
87  
88      }
89  
90      @Test
91      public void testEuler2000() {
92  
93          // this alternate representation of the transform
94          // is from equation 33 in IERS conventions 2003
95          TransformProvider eulerBasedProvider = new TransformProvider() {
96              private final PolynomialNutation zetaA =
97                      new PolynomialNutation(   2.5976176 * Constants.ARC_SECONDS_TO_RADIANS,
98                                             2306.0809506 * Constants.ARC_SECONDS_TO_RADIANS,
99                                                0.3019015 * Constants.ARC_SECONDS_TO_RADIANS,
100                                               0.0179663 * Constants.ARC_SECONDS_TO_RADIANS,
101                                              -0.0000327 * Constants.ARC_SECONDS_TO_RADIANS,
102                                              -0.0000002 * Constants.ARC_SECONDS_TO_RADIANS);
103             private final PolynomialNutation thetaA =
104                     new PolynomialNutation(0.0,
105                                            2004.1917476 * Constants.ARC_SECONDS_TO_RADIANS,
106                                              -0.4269353 * Constants.ARC_SECONDS_TO_RADIANS,
107                                              -0.0418251 * Constants.ARC_SECONDS_TO_RADIANS,
108                                              -0.0000601 * Constants.ARC_SECONDS_TO_RADIANS,
109                                              -0.0000001 * Constants.ARC_SECONDS_TO_RADIANS);
110             private final PolynomialNutation zA =
111                     new PolynomialNutation(  -2.5976176 * Constants.ARC_SECONDS_TO_RADIANS,
112                                            2306.0803226 * Constants.ARC_SECONDS_TO_RADIANS,
113                                               1.0947790 * Constants.ARC_SECONDS_TO_RADIANS,
114                                               0.0182273 * Constants.ARC_SECONDS_TO_RADIANS,
115                                               0.0000470 * Constants.ARC_SECONDS_TO_RADIANS,
116                                              -0.0000003 * Constants.ARC_SECONDS_TO_RADIANS);
117 
118             public Transform getTransform(AbsoluteDate date) {
119                 final double tc = IERSConventions.IERS_2003.evaluateTC(date);
120                 final Rotation r1 = new Rotation(Vector3D.PLUS_K,  zA.value(tc), RotationConvention.VECTOR_OPERATOR);
121                 final Rotation r2 = new Rotation(Vector3D.PLUS_J, -thetaA.value(tc), RotationConvention.VECTOR_OPERATOR);
122                 final Rotation r3 = new Rotation(Vector3D.PLUS_K,  zetaA.value(tc), RotationConvention.VECTOR_OPERATOR);
123                 return new Transform(date, r1.compose(r2.compose(r3,
124                                                                  RotationConvention.VECTOR_OPERATOR),
125                                                       RotationConvention.VECTOR_OPERATOR));
126             }
127 
128             public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
129                 throw new UnsupportedOperationException("never called in this test");
130             }
131 
132         };
133 
134         MODProvider modProvider = new MODProvider(IERSConventions.IERS_2003,
135                 DataContext.getDefault().getTimeScales());
136 
137         for (double dt = -Constants.JULIAN_CENTURY; dt < Constants.JULIAN_CENTURY; dt += 50 * Constants.JULIAN_DAY) {
138             AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(dt);
139             Transform t = new Transform(date,
140                                         modProvider.getTransform(date).getInverse(),
141                                         eulerBasedProvider.getTransform(date));
142             Assertions.assertEquals(0, t.getRotation().getAngle(), 6.6e-13);
143         }
144 
145     }
146 
147     @Test
148     public void testAASReferenceLEO() {
149 
150         // this reference test has been extracted from the following paper:
151         // Implementation Issues Surrounding the New IAU Reference Systems for Astrodynamics
152         // David A. Vallado, John H. Seago, P. Kenneth Seidelmann
153         // http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf
154         AbsoluteDate t0 = new AbsoluteDate(new DateComponents(2004, 04, 06),
155                                            new TimeComponents(07, 51, 28.386009),
156                                            TimeScalesFactory.getUTC());
157 
158         Transform tt = FramesFactory.getGCRF().getTransformTo(FramesFactory.getMOD(IERSConventions.IERS_1996), t0);
159         //GCRF iau76 w corr
160         PVCoordinates pvGCRFiau76 =
161             new PVCoordinates(new Vector3D(5102508.9579, 6123011.4007, 6378136.9282),
162                               new Vector3D(-4743.220157, 790.536497, 5533.755727));
163         double norm = pvGCRFiau76.getPosition().getNorm();
164 
165         //MOD iau76 w corr
166         PVCoordinates pvMODiau76Wcorr =
167             new PVCoordinates(new Vector3D(5094028.3745, 6127870.8164, 6380248.5164),
168                               new Vector3D(-4746.263052, 786.014045, 5531.790562));
169 
170         checkPV(pvMODiau76Wcorr, tt.transformPVCoordinates(pvGCRFiau76), 9e-12 * norm, 7.3e-7);
171 
172         Transform tf = FramesFactory.getEME2000().getTransformTo(FramesFactory.getMOD(false), t0);
173         //J2000 iau76
174         PVCoordinates pvJ2000iau76 =
175             new PVCoordinates(new Vector3D(5102509.6000, 6123011.5200, 6378136.3000),
176                               new Vector3D(-4743.219600, 790.536600, 5533.756190));
177         //MOD iau76
178         PVCoordinates pvMODiau76 =
179             new PVCoordinates(new Vector3D(5094029.0167, 6127870.9363, 6380247.8885),
180                               new Vector3D(-4746.262495, 786.014149, 5531.791025));
181         checkPV(pvMODiau76, tf.transformPVCoordinates(pvJ2000iau76), 9e-12 * norm, 3.1e-7);
182 
183     }
184 
185     @Test
186     public void testAASReferenceGEO() {
187 
188         // this reference test has been extracted from the following paper:
189         // Implementation Issues Surrounding the New IAU Reference Systems for Astrodynamics
190         // David A. Vallado, John H. Seago, P. Kenneth Seidelmann
191         // http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf
192         AbsoluteDate t0 = new AbsoluteDate(new DateComponents(2004, 06, 01),
193                                            TimeComponents.H00,
194                                            TimeScalesFactory.getUTC());
195 
196         Transform tt = FramesFactory.getGCRF().getTransformTo(FramesFactory.getMOD(IERSConventions.IERS_1996), t0);
197         //GCRF iau76 w corr
198         PVCoordinates pvGCRFiau76 =
199             new PVCoordinates(new Vector3D(-40588150.3649, -11462167.0282, 27143.2028),
200                               new Vector3D(834.787457, -2958.305691, -1.172994));
201         double norm = pvGCRFiau76.getPosition().getNorm();
202 
203         //MOD iau76 w corr
204         PVCoordinates pvMODiau76Wcorr =
205             new PVCoordinates(new Vector3D(-40576822.6395, -11502231.5015, 9733.7842),
206                               new Vector3D(837.708020, -2957.480117, -0.814253));
207         checkPV(pvMODiau76Wcorr, tt.transformPVCoordinates(pvGCRFiau76), 9e-12 * norm, 6.9e-7);
208 
209         Transform tf = FramesFactory.getEME2000().getTransformTo(FramesFactory.getMOD(false), t0);
210         //J2000 iau76
211         PVCoordinates pvJ2000iau76 =
212             new PVCoordinates(new Vector3D(-40588150.3620, -11462167.0280, 27147.6490),
213                               new Vector3D(834.787457, -2958.305691, -1.173016));
214         //MOD iau76
215         PVCoordinates pvMODiau76 =
216             new PVCoordinates(new Vector3D(-40576822.6385, -11502231.5013, 9738.2304),
217                               new Vector3D(837.708020, -2957.480118, -0.814275));
218         checkPV(pvMODiau76, tf.transformPVCoordinates(pvJ2000iau76), 9e-12 * norm, 6.9e-7);
219 
220     }
221 
222     @Test
223     public void testSofaPmat76() {
224 
225         // the reference value has been computed using the March 2012 version of the SOFA library
226         // http://www.iausofa.org/2012_0301_C.html, with the following code
227         //
228         //        double utc1, utc2, tai1, tai2, tt1, tt2, rmatp[3][3];
229         //
230         //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
231         //        utc1  = DJM0 + 53049.0;
232         //        utc2  = 0.0;
233         //        iauUtctai(utc1, utc2, &tai1, &tai2);
234         //        iauTaitt(tai1, tai2, &tt1, &tt2);
235         //
236         //        iauPmat76(tt1, tt2, rmatp);
237         //
238         //        printf("iauPmat76(%.20g, %.20g,rmatp)\n"
239         //               "  --> %.20g %.20g %.20g\n"
240         //               "      %.20g %.20g %.20g\n"
241         //               "      %.20g %.20g %.20g\n",
242         //               tt1, tt2,
243         //               rmatp[0][0], rmatp[0][1], rmatp[0][2],
244         //               rmatp[1][0], rmatp[1][1], rmatp[1][2],
245         //               rmatp[2][0], rmatp[2][1], rmatp[2][2]);
246         //
247         // the output of this test reads:
248         //      iauPmat76(2453049.5, 0.00074287037037037029902,rmatp)
249         //        --> 0.9999994956729996165 -0.00092111268696996330928 -0.00040025637336518803469
250         //            0.00092111268696944094067 0.99999957577560194544 -1.843419633938077413e-07
251         //            0.00040025637336639019806 -1.8433935312187383064e-07 0.99999991989739756004
252 
253         // the SOFA routine iauPmat76 uses the Euler angles zetaA, thetaA and zA, whereas
254         // Orekit uses the "canonical 4-rotation method" with angles epsilon0, psiA, omegaA
255         // and chiA. As seen in the thresholds of the test testEuler1976 above, the two methods
256         // drift with respect to each other at a rate about 2e-12 radians/year (i.e. about 0.42
257         // microarcseconds per year) on the time range around J2000.0). We see the same models
258         // difference in this test as in the testEuler1976 above.
259         // Note that the difference is much smaller for IERS 2003 (see another test below)
260         // because the Euler models for these conventions match much better the 4-rotation angles
261 
262         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
263         Frame mod  = FramesFactory.getFrame(Predefined.MOD_CONVENTIONS_1996);
264         Frame gcrf = FramesFactory.getFrame(Predefined.GCRF);
265         checkRotation(new double[][] {
266             { 0.9999994956729996165,     -0.00092111268696996330928, -0.00040025637336518803469 },
267             { 0.00092111268696944094067,  0.99999957577560194544,    -1.843419633938077413e-07  },
268             { 0.00040025637336639019806, -1.8433935312187383064e-07,  0.99999991989739756004    }
269 
270         }, gcrf.getTransformTo(mod, date), 9e-12);
271 
272     }
273 
274     @Test
275     public void testSofaBp00() {
276 
277         // the reference value has been computed using the March 2012 version of the SOFA library
278         // http://www.iausofa.org/2012_0301_C.html, with the following code
279         //
280         //        double utc1, utc2, tai1, tai2, tt1, tt2, rb[3][3], rp[3][3], rbp[3][3];
281         //
282         //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
283         //        utc1  = DJM0 + 53049.0;
284         //        utc2  = 0.0;
285         //        iauUtctai(utc1, utc2, &tai1, &tai2);
286         //        iauTaitt(tai1, tai2, &tt1, &tt2);
287         //
288         //        iauBp00(tt1, tt2, rb, rp, rbp);
289         //
290         //        printf("iauBp00(%.20g, %.20g, rb, rp, rbp)\n"
291         //               "  rb  --> %.20g %.20g %.20g\n          %.20g %.20g %.20g\n          %.20g %.20g %.20g\n"
292         //               "  rp  --> %.20g %.20g %.20g\n          %.20g %.20g %.20g\n          %.20g %.20g %.20g\n"
293         //               "  rbp --> %.20g %.20g %.20g\n          %.20g %.20g %.20g\n          %.20g %.20g %.20g\n",
294         //               tt1, tt2,
295         //               rb[0][0],  rb[0][1],  rb[0][2],
296         //               rb[1][0],  rb[1][1],  rb[1][2],
297         //               rb[2][0],  rb[2][1],  rb[2][2],
298         //               rp[0][0],  rp[0][1],  rp[0][2],
299         //               rp[1][0],  rp[1][1],  rp[1][2],
300         //               rp[2][0],  rp[2][1],  rp[2][2],
301         //               rbp[0][0], rbp[0][1], rbp[0][2],
302         //               rbp[1][0], rbp[1][1], rbp[1][2],
303         //               rbp[2][0], rbp[2][1], rbp[2][2]);
304         //
305         // the output of this test reads:
306         //        iauBp00(2453049.5, 0.00074287037037037029902, rb, rp, rbp)
307         //        rb  --> 0.99999999999999422684 -7.0782797441991980175e-08 8.0562171469761337802e-08
308         //                7.0782794778573375197e-08 0.99999999999999689138 3.3060414542221364117e-08
309         //                -8.0562173809869716745e-08 -3.3060408839805516801e-08 0.99999999999999622524
310         //        rp  --> 0.99999949573309343531 -0.00092105778423522924759 -0.00040023257863225548568
311         //                0.00092105778625203805956 0.99999957582617116092 -1.7927962069881782439e-07
312         //                0.00040023257399096032498 -1.8935780260465051583e-07 0.99999991990692216337
313         //        rbp --> 0.99999949570013624278 -0.00092112855376512230675 -0.00040015204695196122638
314         //                0.00092112856903123034591 0.99999957576097886491 -1.4614501776464880046e-07
315         //                0.00040015201181019732432 -2.2244653837776004327e-07 0.99999991993915571253
316 
317         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
318         Frame mod  = FramesFactory.getFrame(Predefined.MOD_CONVENTIONS_2003);
319         Frame gcrf = FramesFactory.getFrame(Predefined.GCRF);
320         checkRotation(new double[][] {
321             { 0.99999949570013624278,    -0.00092112855376512230675, -0.00040015204695196122638 },
322             { 0.00092112856903123034591,  0.99999957576097886491,    -1.4614501776464880046e-07 },
323             { 0.00040015201181019732432, -2.2244653837776004327e-07,  0.99999991993915571253    }
324 
325         }, gcrf.getTransformTo(mod, date), 4.0e-16);
326 
327     }
328 
329     @Test
330     public void testSofaBp06() {
331 
332         // the reference value has been computed using the March 2012 version of the SOFA library
333         // http://www.iausofa.org/2012_0301_C.html, with the following code
334         //
335         //        double utc1, utc2, tai1, tai2, tt1, tt2, rb[3][3], rp[3][3], rbp[3][3];
336         //
337         //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
338         //        utc1  = DJM0 + 53049.0;
339         //        utc2  = 0.0;
340         //        iauUtctai(utc1, utc2, &tai1, &tai2);
341         //        iauTaitt(tai1, tai2, &tt1, &tt2);
342         //
343         //        iauBp06(tt1, tt2, rb, rp, rbp);
344         //
345         //        printf("iauBp06(%.20g, %.20g, rb, rp, rbp)\n"
346         //               "  rb  --> %.20g %.20g %.20g\n          %.20g %.20g %.20g\n          %.20g %.20g %.20g\n"
347         //               "  rp  --> %.20g %.20g %.20g\n          %.20g %.20g %.20g\n          %.20g %.20g %.20g\n"
348         //               "  rbp --> %.20g %.20g %.20g\n          %.20g %.20g %.20g\n          %.20g %.20g %.20g\n",
349         //               tt1, tt2,
350         //               rb[0][0],  rb[0][1],  rb[0][2],
351         //               rb[1][0],  rb[1][1],  rb[1][2],
352         //               rb[2][0],  rb[2][1],  rb[2][2],
353         //               rp[0][0],  rp[0][1],  rp[0][2],
354         //               rp[1][0],  rp[1][1],  rp[1][2],
355         //               rp[2][0],  rp[2][1],  rp[2][2],
356         //               rbp[0][0], rbp[0][1], rbp[0][2],
357         //               rbp[1][0], rbp[1][1], rbp[1][2],
358         //               rbp[2][0], rbp[2][1], rbp[2][2]);
359         //
360         // the output of this test reads:
361         //      iauBp06(2453049.5, 0.00074287037037037029902, rb, rp, rbp)
362         //        rb  --> 0.99999999999999411582 -7.0783689609715561276e-08 8.0562139776131860839e-08
363         //                7.0783686946376762683e-08 0.99999999999999689138 3.3059437354321374869e-08
364         //                -8.0562142116200574817e-08 -3.3059431692183949281e-08 0.99999999999999622524
365         //        rp  --> 0.99999949573328705821 -0.00092105756953249828464 -0.00040023258864902886648
366         //                0.00092105757159046265214 0.99999957582636889164 -1.7917675327629716543e-07
367         //                0.00040023258391302344184 -1.8946059324815925032e-07 0.99999991990691827759
368         //        rbp --> 0.99999949570032919954 -0.00092112833995494939818 -0.00040015205699952106997
369         //                0.00092112835526181411488 0.9999995757611759295 -1.4604312753574433259e-07
370         //                0.00040015202176394668336 -2.2254835219115420841e-07 0.99999991993915182675
371 
372         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
373         Frame mod  = FramesFactory.getFrame(Predefined.MOD_CONVENTIONS_2010);
374         Frame gcrf = FramesFactory.getFrame(Predefined.GCRF);
375         checkRotation(new double[][] {
376             { 0.99999949570032919954,    -0.00092112833995494939818, -0.00040015205699952106997 },
377             { 0.00092112835526181411488,  0.9999995757611759295,     -1.4604312753574433259e-07 },
378             { 0.00040015202176394668336, -2.2254835219115420841e-07,  0.99999991993915182675    }
379 
380         }, gcrf.getTransformTo(mod, date), 1.1e-12);
381 
382     }
383 
384     @Test
385     public void testMOD1976vs2006() {
386 
387         final Frame mod1976 = FramesFactory.getMOD(IERSConventions.IERS_1996);
388         final Frame mod2006 = FramesFactory.getMOD(IERSConventions.IERS_2010);
389         for (double dt = 0; dt < 10 * Constants.JULIAN_YEAR; dt += 10 * Constants.JULIAN_DAY) {
390             AbsoluteDate date = new AbsoluteDate(AbsoluteDate.J2000_EPOCH, dt);
391             double delta = mod1976.getTransformTo(mod2006, date).getRotation().getAngle();
392             // MOD2006 and MOD2000 are similar to about 33 milli-arcseconds between 2000 and 2010
393             Assertions.assertEquals(0.0, delta, 2.0e-7);
394         }
395     }
396 
397     @Test
398     public void testMOD2000vs2006() {
399 
400         final Frame mod2000 = FramesFactory.getMOD(IERSConventions.IERS_2003);
401         final Frame mod2006 = FramesFactory.getMOD(IERSConventions.IERS_2010);
402         for (double dt = 0; dt < 10 * Constants.JULIAN_YEAR; dt += 10 * Constants.JULIAN_DAY) {
403             AbsoluteDate date = new AbsoluteDate(AbsoluteDate.J2000_EPOCH, dt);
404             double delta = mod2000.getTransformTo(mod2006, date).getRotation().getAngle();
405             // MOD2006 and MOD2000 are similar to about 0.15 milli-arcseconds between 2000 and 2010
406             Assertions.assertEquals(0.0, delta, 7.2e-10);
407         }
408     }
409 
410     @BeforeEach
411     public void setUp() {
412         Utils.setDataRoot("compressed-data");
413     }
414 
415     private void checkPV(PVCoordinates reference,
416                          PVCoordinates result,
417                          double positionThreshold,
418                          double velocityThreshold) {
419 
420         Vector3D dP = result.getPosition().subtract(reference.getPosition());
421         Vector3D dV = result.getVelocity().subtract(reference.getVelocity());
422         Assertions.assertEquals(0, dP.getNorm(), positionThreshold);
423         Assertions.assertEquals(0, dV.getNorm(), velocityThreshold);
424     }
425 
426     private void checkRotation(double[][] reference, Transform t, double epsilon) {
427         double[][] mat = t.getRotation().getMatrix();
428         for (int i = 0; i < 3; ++i) {
429             for (int j = 0; j < 3; ++j) {
430                 Assertions.assertEquals(reference[i][j], mat[i][j], epsilon);
431             }
432         }
433     }
434 
435 }