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.orbits;
18  
19  import org.hamcrest.MatcherAssert;
20  import org.hipparchus.analysis.UnivariateFunction;
21  import org.hipparchus.analysis.differentiation.DSFactory;
22  import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
23  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.linear.MatrixUtils;
26  import org.hipparchus.linear.RealMatrixPreservingVisitor;
27  import org.hipparchus.util.FastMath;
28  import org.hipparchus.util.MathUtils;
29  import org.junit.jupiter.api.*;
30  import org.junit.jupiter.params.ParameterizedTest;
31  import org.junit.jupiter.params.provider.EnumSource;
32  import org.orekit.Utils;
33  import org.orekit.errors.OrekitIllegalArgumentException;
34  import org.orekit.errors.OrekitMessages;
35  import org.orekit.frames.Frame;
36  import org.orekit.frames.FramesFactory;
37  import org.orekit.frames.Transform;
38  import org.orekit.time.AbsoluteDate;
39  import org.orekit.time.TimeScalesFactory;
40  import org.orekit.utils.Constants;
41  import org.orekit.utils.PVCoordinates;
42  import org.orekit.utils.TimeStampedPVCoordinates;
43  
44  import java.util.function.Function;
45  
46  import static org.orekit.OrekitMatchers.relativelyCloseTo;
47  
48  
49  class CircularOrbitTest {
50  
51      // Computation date
52      private AbsoluteDate date;
53  
54      // Body mu
55      private double mu;
56  
57      @ParameterizedTest
58      @EnumSource(PositionAngleType.class)
59      void testWithCachedPositionAngleType(final PositionAngleType positionAngleType) {
60          // GIVEN
61          final Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
62          final Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
63          final PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
64          final double muEarth = 3.9860047e14;
65          final CartesianOrbit cartesianOrbit = new CartesianOrbit(pvCoordinates, FramesFactory.getEME2000(), date, muEarth);
66          final CircularOrbit circularOrbit = new CircularOrbit(cartesianOrbit);
67          // WHEN
68          final CircularOrbit orbit = circularOrbit.withCachedPositionAngleType(positionAngleType);
69          // THEN
70          Assertions.assertEquals(circularOrbit.getFrame(), orbit.getFrame());
71          Assertions.assertEquals(circularOrbit.getDate(), orbit.getDate());
72          Assertions.assertEquals(circularOrbit.getMu(), orbit.getMu());
73          final Vector3D relativePosition = circularOrbit.getPosition(orbit.getFrame()).subtract(
74                  orbit.getPosition());
75          Assertions.assertEquals(0., relativePosition.getNorm(), 1e-6);
76          Assertions.assertEquals(circularOrbit.hasNonKeplerianAcceleration(),
77                  orbit.hasNonKeplerianAcceleration());
78      }
79  
80      @Test
81      void testNonKeplerianAcceleration() {
82          // GIVEN
83          final PVCoordinates pvCoordinates = new PVCoordinates(new Vector3D(1, 2, 3),
84                  Vector3D.MINUS_K.scalarMultiply(0.1), new Vector3D(0, 1));
85          final CartesianOrbit cartesianOrbit = new CartesianOrbit(pvCoordinates, FramesFactory.getEME2000(),
86                  AbsoluteDate.ARBITRARY_EPOCH, 1.);
87          final CircularOrbit circularOrbit = (CircularOrbit) OrbitType.CIRCULAR.convertType(cartesianOrbit);
88          // WHEN
89          final Vector3D nonKeplerianAcceleration = circularOrbit.nonKeplerianAcceleration();
90          // THEN
91          final Vector3D expectedVector = computeLegacyNonKeplerianAcceleration(circularOrbit);
92          Assertions.assertArrayEquals(expectedVector.toArray(), nonKeplerianAcceleration.toArray(), 1e-10);
93      }
94  
95      private Vector3D computeLegacyNonKeplerianAcceleration(final Orbit orbit) {
96          final double[][] dCdP = new double[6][6];
97          final PositionAngleType positionAngleType = PositionAngleType.MEAN;
98          orbit.getJacobianWrtParameters(positionAngleType, dCdP);
99  
100         final double[] derivatives = new double[6];
101         orbit.getType().mapOrbitToArray(orbit, positionAngleType, new double[6], derivatives);
102         derivatives[5] -= orbit.getKeplerianMeanMotion();
103         final double nonKeplerianAx = dCdP[3][0] * derivatives[0]  + dCdP[3][1] * derivatives[1] + dCdP[3][2] * derivatives[2] +
104                 dCdP[3][3] * derivatives[3] + dCdP[3][4] * derivatives[4] + dCdP[3][5] * derivatives[5];
105         final double nonKeplerianAy = dCdP[4][0] * derivatives[0]  + dCdP[4][1] * derivatives[1] + dCdP[4][2] * derivatives[2] +
106                 dCdP[4][3] * derivatives[3] + dCdP[4][4] * derivatives[4] + dCdP[4][5] * derivatives[5];
107         final double nonKeplerianAz = dCdP[5][0] * derivatives[0]  + dCdP[5][1] * derivatives[1] + dCdP[5][2] * derivatives[2] +
108                 dCdP[5][3] * derivatives[3] + dCdP[5][4] * derivatives[4] + dCdP[5][5] * derivatives[5];
109 
110         return new Vector3D(nonKeplerianAx, nonKeplerianAy, nonKeplerianAz);
111     }
112 
113     @ParameterizedTest
114     @EnumSource(PositionAngleType.class)
115     void testInFrameKeplerian(final PositionAngleType positionAngleType) {
116         testTemplateInFrame(Vector3D.ZERO, positionAngleType);
117     }
118 
119     @Test
120     void testInFrameNonKeplerian() {
121         testTemplateInFrame(Vector3D.PLUS_K, PositionAngleType.TRUE);
122     }
123 
124     private void testTemplateInFrame(final Vector3D acceleration, final PositionAngleType positionAngleType) {
125         // GIVEN
126         final Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
127         final Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
128         final PVCoordinates pvCoordinates = new PVCoordinates(position, velocity, acceleration);
129         final double muEarth = 3.9860047e14;
130         final CartesianOrbit cartesianOrbit = new CartesianOrbit(pvCoordinates, FramesFactory.getEME2000(), date, muEarth);
131         final CircularOrbit circularOrbit = new CircularOrbit(cartesianOrbit).withCachedPositionAngleType(positionAngleType);
132         // WHEN
133         final CircularOrbit orbitWithOtherFrame = circularOrbit.inFrame(FramesFactory.getGCRF());
134         // THEN
135         Assertions.assertNotEquals(circularOrbit.getFrame(), orbitWithOtherFrame.getFrame());
136         Assertions.assertEquals(circularOrbit.getDate(), orbitWithOtherFrame.getDate());
137         Assertions.assertEquals(circularOrbit.getMu(), orbitWithOtherFrame.getMu());
138         final Vector3D relativePosition = circularOrbit.getPosition(orbitWithOtherFrame.getFrame()).subtract(
139                 orbitWithOtherFrame.getPosition());
140         Assertions.assertEquals(0., relativePosition.getNorm(), 1e-6);
141         Assertions.assertEquals(circularOrbit.hasNonKeplerianAcceleration(),
142                 orbitWithOtherFrame.hasNonKeplerianAcceleration());
143     }
144 
145     @Test
146     void testCircularToEquinoctialEll() {
147 
148         double ix = 1.200e-04;
149         double iy = -1.16e-04;
150         double i  = 2 * FastMath.asin(FastMath.sqrt((ix * ix + iy * iy) / 4));
151         double raan = FastMath.atan2(iy, ix);
152 
153         // elliptic orbit
154         CircularOrbit circ =
155                 new CircularOrbit(42166712.0, 0.5, -0.5, i, raan,
156                         5.300 - raan, PositionAngleType.MEAN,
157                         FramesFactory.getEME2000(), date, mu);
158         Vector3D pos = circ.getPosition();
159         Vector3D vit = circ.getVelocity();
160 
161         PVCoordinates pvCoordinates = new PVCoordinates( pos, vit);
162 
163         EquinoctialOrbit param = new EquinoctialOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
164         Assertions.assertEquals(param.getA(),  circ.getA(), Utils.epsilonTest * circ.getA());
165         Assertions.assertEquals(param.getEquinoctialEx(), circ.getEquinoctialEx(), Utils.epsilonE * FastMath.abs(circ.getE()));
166         Assertions.assertEquals(param.getEquinoctialEy(), circ.getEquinoctialEy(), Utils.epsilonE * FastMath.abs(circ.getE()));
167         Assertions.assertEquals(param.getHx(), circ.getHx(), Utils.epsilonAngle * FastMath.abs(circ.getI()));
168         Assertions.assertEquals(param.getHy(), circ.getHy(), Utils.epsilonAngle * FastMath.abs(circ.getI()));
169         Assertions.assertEquals(MathUtils.normalizeAngle(param.getLv(), circ.getLv()), circ.getLv(), Utils.epsilonAngle * FastMath.abs(circ.getLv()));
170 
171     }
172 
173     @Test
174     void testCircularToEquinoctialCirc() {
175 
176         double ix = 1.200e-04;
177         double iy = -1.16e-04;
178         double i  = 2 * FastMath.asin(FastMath.sqrt((ix * ix + iy * iy) / 4));
179         double raan = FastMath.atan2(iy, ix);
180 
181         // circular orbit
182         EquinoctialOrbit circCir =
183                 new EquinoctialOrbit(42166712.0, 0.1e-10, -0.1e-10, i, raan,
184                         5.300 - raan, PositionAngleType.MEAN,
185                         FramesFactory.getEME2000(), date, mu);
186         Vector3D posCir = circCir.getPosition();
187         Vector3D vitCir = circCir.getVelocity();
188 
189         PVCoordinates pvCoordinates = new PVCoordinates( posCir, vitCir);
190 
191         EquinoctialOrbit paramCir = new EquinoctialOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
192         Assertions.assertEquals(paramCir.getA(), circCir.getA(), Utils.epsilonTest * circCir.getA());
193         Assertions.assertEquals(paramCir.getEquinoctialEx(), circCir.getEquinoctialEx(), Utils.epsilonEcir * FastMath.abs(circCir.getE()));
194         Assertions.assertEquals(paramCir.getEquinoctialEy(), circCir.getEquinoctialEy(), Utils.epsilonEcir * FastMath.abs(circCir.getE()));
195         Assertions.assertEquals(paramCir.getHx(), circCir.getHx(), Utils.epsilonAngle * FastMath.abs(circCir.getI()));
196         Assertions.assertEquals(paramCir.getHy(), circCir.getHy(), Utils.epsilonAngle * FastMath.abs(circCir.getI()));
197         Assertions.assertEquals(MathUtils.normalizeAngle(paramCir.getLv(), circCir.getLv()), circCir.getLv(), Utils.epsilonAngle * FastMath.abs(circCir.getLv()));
198 
199     }
200 
201     @Test
202     void testCircularToCartesian() {
203 
204         double ix = 1.200e-04;
205         double iy = -1.16e-04;
206         double i  = 2 * FastMath.asin(FastMath.sqrt((ix * ix + iy * iy) / 4));
207         double raan = FastMath.atan2(iy, ix);
208         double cosRaan = FastMath.cos(raan);
209         double sinRaan = FastMath.sin(raan);
210         double exTilde = -7.900e-6;
211         double eyTilde = 1.100e-4;
212         double ex = exTilde * cosRaan + eyTilde * sinRaan;
213         double ey = eyTilde * cosRaan - exTilde * sinRaan;
214 
215         CircularOrbit circ=
216                 new CircularOrbit(42166712.0, ex, ey, i, raan,
217                         5.300 - raan, PositionAngleType.MEAN,
218                         FramesFactory.getEME2000(), date, mu);
219         Vector3D pos = circ.getPosition();
220         Vector3D vel = circ.getVelocity();
221 
222         // check 1/a = 2/r  - V2/mu
223         double r = pos.getNorm();
224         double v = vel.getNorm();
225         Assertions.assertEquals(2 / r - v * v / mu, 1 / circ.getA(), 1.0e-7);
226 
227         Assertions.assertEquals( 0.233745668678733e+08, pos.getX(), Utils.epsilonTest * r);
228         Assertions.assertEquals(-0.350998914352669e+08, pos.getY(), Utils.epsilonTest * r);
229         Assertions.assertEquals(-0.150053723123334e+04, pos.getZ(), Utils.epsilonTest * r);
230 
231         Assertions.assertEquals(2558.7096558809967, vel.getX(), Utils.epsilonTest * v);
232         Assertions.assertEquals(1704.1586039092576, vel.getY(), Utils.epsilonTest * v);
233         Assertions.assertEquals(   0.5013093577879, vel.getZ(), Utils.epsilonTest * v);
234 
235     }
236 
237     @Test
238     void testCircularToKeplerian() {
239 
240         double ix   =  1.20e-4;
241         double iy   = -1.16e-4;
242         double i    = 2 * FastMath.asin(FastMath.sqrt((ix * ix + iy * iy) / 4));
243         double raan = FastMath.atan2(iy, ix);
244         double cosRaan = FastMath.cos(raan);
245         double sinRaan = FastMath.sin(raan);
246         double exTilde = -7.900e-6;
247         double eyTilde = 1.100e-4;
248         double ex = exTilde * cosRaan + eyTilde * sinRaan;
249         double ey = eyTilde * cosRaan - exTilde * sinRaan;
250 
251         CircularOrbit circ=
252                 new CircularOrbit(42166712.0, ex, ey, i, raan,
253                         5.300 - raan, PositionAngleType.MEAN,
254                         FramesFactory.getEME2000(), date, mu);
255         KeplerianOrbit kep = new KeplerianOrbit(circ);
256 
257         Assertions.assertEquals(42166712.000, circ.getA(), Utils.epsilonTest * kep.getA());
258         Assertions.assertEquals(0.110283316961361e-03, kep.getE(), Utils.epsilonE * FastMath.abs(kep.getE()));
259         Assertions.assertEquals(0.166901168553917e-03, kep.getI(),
260                 Utils.epsilonAngle * FastMath.abs(kep.getI()));
261         Assertions.assertEquals(MathUtils.normalizeAngle(-3.87224326008837, kep.getPerigeeArgument()),
262                 kep.getPerigeeArgument(),
263                 Utils.epsilonTest * FastMath.abs(kep.getPerigeeArgument()));
264         Assertions.assertEquals(MathUtils.normalizeAngle(5.51473467358854, kep.getRightAscensionOfAscendingNode()),
265                 kep.getRightAscensionOfAscendingNode(),
266                 Utils.epsilonTest * FastMath.abs(kep.getRightAscensionOfAscendingNode()));
267         Assertions.assertEquals(MathUtils.normalizeAngle(3.65750858649982, kep.getMeanAnomaly()),
268                 kep.getMeanAnomaly(),
269                 Utils.epsilonTest * FastMath.abs(kep.getMeanAnomaly()));
270 
271     }
272 
273     @Test
274     void testHyperbolic1() {
275         try {
276             new CircularOrbit(-42166712.0, 1.9, 0.5, 0.01, -0.02, 5.300,
277                     PositionAngleType.MEAN,  FramesFactory.getEME2000(), date, mu);
278             Assertions.fail("an exception should have been thrown");
279         } catch (OrekitIllegalArgumentException oe) {
280             Assertions.assertEquals(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS, oe.getSpecifier());
281         }
282     }
283 
284     @Test
285     void testHyperbolic2() {
286         Orbit orbit = new KeplerianOrbit(-42166712.0, 1.9, 0.5, 0.01, -0.02, 5.300,
287                 PositionAngleType.MEAN,  FramesFactory.getEME2000(), date, mu);
288         try {
289             new CircularOrbit(orbit.getPVCoordinates(), orbit.getFrame(), orbit.getMu());
290             Assertions.fail("an exception should have been thrown");
291         } catch (OrekitIllegalArgumentException oe) {
292             Assertions.assertEquals(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS, oe.getSpecifier());
293         }
294     }
295 
296     @Test
297     void testAnomalyEll() {
298 
299         // elliptic orbit
300         Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
301         Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
302 
303         PVCoordinates pvCoordinates = new PVCoordinates( position, velocity);
304 
305         CircularOrbit  p   = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
306         KeplerianOrbit kep = new KeplerianOrbit(p);
307 
308         double e       = p.getE();
309         double eRatio  = FastMath.sqrt((1 - e) / (1 + e));
310         double raan    = kep.getRightAscensionOfAscendingNode();
311         double paPraan = kep.getPerigeeArgument() + raan;
312 
313         double lv = 1.1;
314         // formulations for elliptic case
315         double lE = 2 * FastMath.atan(eRatio * FastMath.tan((lv - paPraan) / 2)) + paPraan;
316         double lM = lE - e * FastMath.sin(lE - paPraan);
317 
318         p = new CircularOrbit(p.getA() , p.getCircularEx(), p.getCircularEy(),
319                 p.getRightAscensionOfAscendingNode(),
320                 p.getAlphaV(), lv - raan, PositionAngleType.TRUE, p.getFrame(), date, mu);
321         Assertions.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
322         Assertions.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
323         Assertions.assertEquals(p.getAlphaM() + raan, lM, Utils.epsilonAngle * FastMath.abs(lM));
324         p = new CircularOrbit(p.getA() , p.getCircularEx(), p.getCircularEy(),
325                 p.getRightAscensionOfAscendingNode(),
326                 p.getAlphaV(), 0, PositionAngleType.TRUE, p.getFrame(), date, mu);
327 
328 
329         p = new CircularOrbit(p.getA() , p.getCircularEx(), p.getCircularEy(),
330                 p.getRightAscensionOfAscendingNode(),
331                 p.getAlphaV(), lE - raan, PositionAngleType.ECCENTRIC, p.getFrame(), date, mu);
332         Assertions.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
333         Assertions.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
334         Assertions.assertEquals(p.getAlphaM() + raan, lM, Utils.epsilonAngle * FastMath.abs(lM));
335         p = new CircularOrbit(p.getA() , p.getCircularEx(), p.getCircularEy(),
336                 p.getRightAscensionOfAscendingNode(),
337                 p.getAlphaV(), 0, PositionAngleType.TRUE, p.getFrame(), date, mu);
338 
339         p = new CircularOrbit(p.getA() , p.getCircularEx(), p.getCircularEy(),
340                 p.getRightAscensionOfAscendingNode(),
341                 p.getAlphaV(), lM - raan, PositionAngleType.MEAN, p.getFrame(), date, mu);
342         Assertions.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
343         Assertions.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
344         Assertions.assertEquals(p.getAlphaM() + raan, lM, Utils.epsilonAngle * FastMath.abs(lM));
345 
346     }
347 
348     @Test
349     void testAnomalyCirc() {
350 
351         Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
352         Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
353         PVCoordinates pvCoordinates = new PVCoordinates( position, velocity);
354         CircularOrbit  p = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
355         double raan = p.getRightAscensionOfAscendingNode();
356 
357         // circular orbit
358         p = new CircularOrbit(p.getA() , 0, 0, p.getRightAscensionOfAscendingNode(),
359                 p.getAlphaV(), p.getAlphaV(), PositionAngleType.TRUE, p.getFrame(), date, mu);
360 
361         double lv = 1.1;
362         double lE = lv;
363         double lM = lE;
364 
365         p = new CircularOrbit(p.getA() , p.getCircularEx(), p.getCircularEy(),
366                 p.getRightAscensionOfAscendingNode(),
367                 p.getAlphaV(), lv - raan, PositionAngleType.TRUE, p.getFrame(), date, mu);
368         Assertions.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
369         Assertions.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
370         Assertions.assertEquals(p.getAlphaM() + raan, lM, Utils.epsilonAngle * FastMath.abs(lM));
371         p = new CircularOrbit(p.getA() , p.getCircularEx(), p.getCircularEy(),
372                 p.getRightAscensionOfAscendingNode(),
373                 p.getAlphaV(), 0, PositionAngleType.TRUE, p.getFrame(), date, mu);
374 
375         p = new CircularOrbit(p.getA() , p.getCircularEx(), p.getCircularEy(),
376                 p.getRightAscensionOfAscendingNode(),
377                 p.getAlphaV(), lE - raan, PositionAngleType.ECCENTRIC, p.getFrame(), date, mu);
378 
379         Assertions.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
380         Assertions.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
381         Assertions.assertEquals(p.getAlphaM() + raan, lM, Utils.epsilonAngle * FastMath.abs(lM));
382         p = new CircularOrbit(p.getA() , p.getCircularEx(), p.getCircularEy(),
383                 p.getRightAscensionOfAscendingNode(),
384                 p.getAlphaV(), 0, PositionAngleType.TRUE, p.getFrame(), date, mu);
385 
386         p = new CircularOrbit(p.getA() , p.getCircularEx(), p.getCircularEy(),
387                 p.getRightAscensionOfAscendingNode(),
388                 p.getAlphaV(), lM - raan, PositionAngleType.MEAN, p.getFrame(), date, mu);
389         Assertions.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
390         Assertions.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
391         Assertions.assertEquals(p.getAlphaM() + raan, lM, Utils.epsilonAngle * FastMath.abs(lM));
392 
393     }
394 
395     @Test
396     void testPositionVelocityNormsEll() {
397 
398         // elliptic and non equatorial (i retrograde) orbit
399         double hx =  1.2;
400         double hy =  2.1;
401         double i  = 2 * FastMath.atan(FastMath.sqrt(hx * hx + hy * hy));
402         double raan = FastMath.atan2(hy, hx);
403         CircularOrbit p =
404                 new CircularOrbit(42166712.0, 0.5, -0.5, i, raan,
405                         0.67 - raan, PositionAngleType.TRUE,
406                         FramesFactory.getEME2000(), date, mu);
407 
408         double ex = p.getEquinoctialEx();
409         double ey = p.getEquinoctialEy();
410         double lv = p.getLv();
411         double ksi     = 1 + ex * FastMath.cos(lv) + ey * FastMath.sin(lv);
412         double nu      = ex * FastMath.sin(lv) - ey * FastMath.cos(lv);
413         double epsilon = FastMath.sqrt(1 - ex * ex - ey * ey);
414 
415         double a  = p.getA();
416         double na = FastMath.sqrt(mu / a);
417 
418         Assertions.assertEquals(a * epsilon * epsilon / ksi,
419                 p.getPosition().getNorm(),
420                 Utils.epsilonTest * FastMath.abs(p.getPosition().getNorm()));
421         Assertions.assertEquals(na * FastMath.sqrt(ksi * ksi + nu * nu) / epsilon,
422                 p.getVelocity().getNorm(),
423                 Utils.epsilonTest * FastMath.abs(p.getVelocity().getNorm()));
424 
425     }
426 
427     @Test
428     void testNumericalIssue25() {
429         Vector3D position = new Vector3D(3782116.14107698, 416663.11924914, 5875541.62103057);
430         Vector3D velocity = new Vector3D(-6349.7848910501, 288.4061811651, 4066.9366759691);
431         CircularOrbit orbit = new CircularOrbit(new PVCoordinates(position, velocity),
432                 FramesFactory.getEME2000(),
433                 new AbsoluteDate("2004-01-01T23:00:00.000",
434                         TimeScalesFactory.getUTC()),
435                 3.986004415E14);
436         Assertions.assertEquals(0.0, orbit.getE(), 2.0e-14);
437     }
438 
439     @Test
440     void testPerfectlyEquatorial() {
441         Vector3D position = new Vector3D(-7293947.695148368, 5122184.668436634, 0.0);
442         Vector3D velocity = new Vector3D(-3890.4029433398, -5369.811285264604, 0.0);
443         CircularOrbit orbit = new CircularOrbit(new PVCoordinates(position, velocity),
444                 FramesFactory.getEME2000(),
445                 new AbsoluteDate("2004-01-01T23:00:00.000",
446                         TimeScalesFactory.getUTC()),
447                 3.986004415E14);
448         Assertions.assertEquals(0.0, orbit.getI(), 2.0e-14);
449         Assertions.assertEquals(0.0, orbit.getRightAscensionOfAscendingNode(), 2.0e-14);
450     }
451 
452     @Test
453     void testPositionVelocityNormsCirc() {
454 
455         // elliptic and non equatorial (i retrograde) orbit
456         double hx =  0.1e-8;
457         double hy =  0.1e-8;
458         double i  = 2 * FastMath.atan(FastMath.sqrt(hx * hx + hy * hy));
459         double raan = FastMath.atan2(hy, hx);
460         CircularOrbit pCirEqua =
461                 new CircularOrbit(42166712.0, 0.1e-8, 0.1e-8, i, raan,
462                         0.67 - raan, PositionAngleType.TRUE,
463                         FramesFactory.getEME2000(), date, mu);
464 
465         double ex = pCirEqua.getEquinoctialEx();
466         double ey = pCirEqua.getEquinoctialEy();
467         double lv = pCirEqua.getLv();
468         double ksi     = 1 + ex * FastMath.cos(lv) + ey * FastMath.sin(lv);
469         double nu      = ex * FastMath.sin(lv) - ey * FastMath.cos(lv);
470         double epsilon = FastMath.sqrt(1 - ex * ex - ey * ey);
471 
472         double a  = pCirEqua.getA();
473         double na = FastMath.sqrt(mu / a);
474 
475         Assertions.assertEquals(a * epsilon * epsilon / ksi,
476                 pCirEqua.getPosition().getNorm(),
477                 Utils.epsilonTest * FastMath.abs(pCirEqua.getPosition().getNorm()));
478         Assertions.assertEquals(na * FastMath.sqrt(ksi * ksi + nu * nu) / epsilon,
479                 pCirEqua.getVelocity().getNorm(),
480                 Utils.epsilonTest * FastMath.abs(pCirEqua.getVelocity().getNorm()));
481     }
482 
483     @Test
484     void testGeometryEll() {
485 
486         // elliptic and non equatorial (i retrograde) orbit
487         double hx =  1.2;
488         double hy =  2.1;
489         double i  = 2 * FastMath.atan(FastMath.sqrt(hx * hx + hy * hy));
490         double raan = FastMath.atan2(hy, hx);
491         CircularOrbit p =
492                 new CircularOrbit(42166712.0, 0.5, -0.5, i, raan,
493                         0.67 - raan, PositionAngleType.TRUE,
494                         FramesFactory.getEME2000(), date, mu);
495 
496         Vector3D position = p.getPosition();
497         Vector3D velocity = p.getVelocity();
498         Vector3D momentum = p.getPVCoordinates().getMomentum().normalize();
499 
500         double apogeeRadius  = p.getA() * (1 + p.getE());
501         double perigeeRadius = p.getA() * (1 - p.getE());
502 
503         for (double alphaV = 0; alphaV <= 2 * FastMath.PI; alphaV += 2 * FastMath.PI/100.) {
504             p = new CircularOrbit(p.getA() , p.getCircularEx(), p.getCircularEy(), p.getI(),
505                     p.getRightAscensionOfAscendingNode(),
506                     alphaV, PositionAngleType.TRUE, p.getFrame(), date, mu);
507             position = p.getPosition();
508             // test if the norm of the position is in the range [perigee radius, apogee radius]
509             // Warning: these tests are without absolute value by choice
510             Assertions.assertTrue((position.getNorm() - apogeeRadius)  <= (  apogeeRadius * Utils.epsilonTest));
511             Assertions.assertTrue((position.getNorm() - perigeeRadius) >= (- perigeeRadius * Utils.epsilonTest));
512 
513             position= position.normalize();
514             velocity = p.getVelocity();
515             velocity= velocity.normalize();
516 
517             // at this stage of computation, all the vectors (position, velocity and momemtum) are normalized here
518 
519             // test of orthogonality between position and momentum
520             Assertions.assertTrue(FastMath.abs(Vector3D.dotProduct(position, momentum)) < Utils.epsilonTest);
521             // test of orthogonality between velocity and momentum
522             Assertions.assertTrue(FastMath.abs(Vector3D.dotProduct(velocity, momentum)) < Utils.epsilonTest);
523         }
524 
525     }
526 
527     @Test
528     void testGeometryCirc() {
529 
530         //  circular and equatorial orbit
531         double hx =  0.1e-8;
532         double hy =  0.1e-8;
533         double i  = 2 * FastMath.atan(FastMath.sqrt(hx * hx + hy * hy));
534         double raan = FastMath.atan2(hy, hx);
535         CircularOrbit pCirEqua =
536                 new CircularOrbit(42166712.0, 0.1e-8, 0.1e-8, i, raan,
537                         0.67 - raan, PositionAngleType.TRUE,
538                         FramesFactory.getEME2000(), date, mu);
539 
540         Vector3D position = pCirEqua.getPosition();
541         Vector3D velocity = pCirEqua.getVelocity();
542         Vector3D momentum = pCirEqua.getPVCoordinates().getMomentum().normalize();
543 
544         double apogeeRadius  = pCirEqua.getA() * (1 + pCirEqua.getE());
545         double perigeeRadius = pCirEqua.getA() * (1 - pCirEqua.getE());
546         // test if apogee equals perigee
547         Assertions.assertEquals(perigeeRadius, apogeeRadius, 1.e+4 * Utils.epsilonTest * apogeeRadius);
548 
549         for (double alphaV = 0; alphaV <= 2 * FastMath.PI; alphaV += 2 * FastMath.PI/100.) {
550             pCirEqua = new CircularOrbit(pCirEqua.getA() , pCirEqua.getCircularEx(), pCirEqua.getCircularEy(), pCirEqua.getI(),
551                     pCirEqua.getRightAscensionOfAscendingNode(),
552                     alphaV, PositionAngleType.TRUE, pCirEqua.getFrame(), date, mu);
553             position = pCirEqua.getPosition();
554 
555             // test if the norm pf the position is in the range [perigee radius, apogee radius]
556             Assertions.assertTrue((position.getNorm() - apogeeRadius)  <= (  apogeeRadius * Utils.epsilonTest));
557             Assertions.assertTrue((position.getNorm() - perigeeRadius) >= (- perigeeRadius * Utils.epsilonTest));
558 
559             position= position.normalize();
560             velocity = pCirEqua.getVelocity();
561             velocity= velocity.normalize();
562 
563             // at this stage of computation, all the vectors (position, velocity and momemtum) are normalized here
564 
565             // test of orthogonality between position and momentum
566             Assertions.assertTrue(FastMath.abs(Vector3D.dotProduct(position, momentum)) < Utils.epsilonTest);
567             // test of orthogonality between velocity and momentum
568             Assertions.assertTrue(FastMath.abs(Vector3D.dotProduct(velocity, momentum)) < Utils.epsilonTest);
569         }
570     }
571 
572     @Test
573     void testSymmetryEll() {
574 
575         // elliptic and non equatorail orbit
576         Vector3D position = new Vector3D(4512.9, 18260., -5127.);
577         Vector3D velocity = new Vector3D(134664.6, 90066.8, 72047.6);
578         PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
579 
580         CircularOrbit p = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
581 
582         Vector3D positionOffset = p.getPosition();
583         Vector3D velocityOffset = p.getVelocity();
584 
585         positionOffset = positionOffset.subtract(position);
586         velocityOffset = velocityOffset.subtract(velocity);
587 
588         Assertions.assertEquals(0.0, positionOffset.getNorm(), position.getNorm() * Utils.epsilonTest);
589         Assertions.assertEquals(0.0, velocityOffset.getNorm(), velocity.getNorm() * Utils.epsilonTest);
590 
591     }
592 
593     @Test
594     void testSymmetryCir() {
595         // circular and equatorial orbit
596         Vector3D position = new Vector3D(33051.2, 26184.9, -1.3E-5);
597         Vector3D velocity = new Vector3D(-60376.2, 76208., 2.7E-4);
598         PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
599 
600         CircularOrbit p = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
601 
602         Vector3D positionOffset = p.getPosition().subtract(position);
603         Vector3D velocityOffset = p.getVelocity().subtract(velocity);
604 
605         Assertions.assertEquals(0.0, positionOffset.getNorm(), position.getNorm() * Utils.epsilonTest);
606         Assertions.assertEquals(0.0, velocityOffset.getNorm(), velocity.getNorm() * Utils.epsilonTest);
607 
608     }
609 
610     @Test
611     void testNonInertialFrame() throws IllegalArgumentException {
612         Assertions.assertThrows(IllegalArgumentException.class, () -> {
613             Vector3D position = new Vector3D(33051.2, 26184.9, -1.3E-5);
614             Vector3D velocity = new Vector3D(-60376.2, 76208., 2.7E-4);
615             PVCoordinates pvCoordinates = new PVCoordinates( position, velocity);
616             new CircularOrbit(pvCoordinates,
617                     new Frame(FramesFactory.getEME2000(), Transform.IDENTITY, "non-inertial", false),
618                     date, mu);
619         });
620     }
621 
622     @Test
623     void testJacobianReference() {
624 
625         AbsoluteDate dateTca = new AbsoluteDate(2000, 4, 1, 0, 0, 0.000, TimeScalesFactory.getUTC());
626         double mu =  3.986004415e+14;
627         CircularOrbit orbCir = new CircularOrbit(7000000.0, 0.01, -0.02, 1.2, 2.1,
628                 0.7, PositionAngleType.MEAN,
629                 FramesFactory.getEME2000(), dateTca, mu);
630 
631         // the following reference values have been computed using the free software
632         // version 6.2 of the MSLIB fortran library by the following program:
633         //        program cir_jacobian
634         //
635         //        use mslib
636         //        implicit none
637         //
638         //        integer, parameter :: nb = 11
639         //        integer :: i,j
640         //        type(tm_code_retour)      ::  code_retour
641         //
642         //        real(pm_reel), parameter :: mu= 3.986004415e+14_pm_reel
643         //        real(pm_reel),dimension(3)::vit_car,pos_car
644         //        type(tm_orb_cir)::cir
645         //        real(pm_reel), dimension(6,6)::jacob
646         //        real(pm_reel)::norme
647         //
648         //
649         //        cir%a=7000000_pm_reel
650         //        cir%ex=0.01_pm_reel
651         //        cir%ey=-0.02_pm_reel
652         //        cir%i=1.2_pm_reel
653         //        cir%gom=2.1_pm_reel
654         //        cir%pso_M=0.7_pm_reel
655         //
656         //        call mv_cir_car(mu,cir,pos_car,vit_car,code_retour)
657         //        write(*,*)code_retour%valeur
658         //        write(*,1000)pos_car,vit_car
659         //
660         //
661         //        call mu_norme(pos_car,norme,code_retour)
662         //        write(*,*)norme
663         //
664         //        call mv_car_cir (mu, pos_car, vit_car, cir, code_retour, jacob)
665         //        write(*,*)code_retour%valeur
666         //
667         //        write(*,*)"circular = ", cir%a, cir%ex, cir%ey, cir%i, cir%gom, cir%pso_M
668         //
669         //        do i = 1,6
670         //           write(*,*) " ",(jacob(i,j),j=1,6)
671         //        end do
672         //
673         //        1000 format (6(f24.15,1x))
674         //        end program cir_jacobian
675         Vector3D pRef = new Vector3D(-4106905.105389204807580, 3603162.539798960555345, 4439730.167038885876536);
676         Vector3D vRef = new Vector3D(740.132407342422994, -5308.773280141396754, 5250.338353483879473);
677         double[][] jRef = {
678                 { -1.1535467596325562,        1.0120556393573172,        1.2470306024626943,        181.96913090864561,       -1305.2162699469984,        1290.8494448855752      },
679                 { -5.07367368325471104E-008, -1.27870567070456834E-008,  1.31544531338558113E-007, -3.09332106417043592E-005, -9.60781276304445404E-005,  1.91506964883791605E-004 },
680                 { -6.59428471712402018E-008,  1.24561703203882533E-007, -1.41907027322388158E-008,  7.63442601186485441E-005, -1.77446722746170009E-004,  5.99464401287846734E-005 },
681                 {  7.55079920652274275E-008,  4.41606835295069131E-008,  3.40079310688458225E-008,  7.89724635377817962E-005,  4.61868720707717372E-005,  3.55682891687782599E-005 },
682                 { -9.20788748896973282E-008, -5.38521280004949642E-008, -4.14712660805579618E-008,  7.78626692360739821E-005,  4.55378113077967091E-005,  3.50684505810897702E-005 },
683                 {  1.85082436324531617E-008,  1.20506219457886855E-007, -8.31277842285972640E-008,  1.27364008345789645E-004, -1.54770720974742483E-004, -1.78589436862677754E-004 }
684         };
685 
686         PVCoordinates pv = orbCir.getPVCoordinates();
687         Assertions.assertEquals(0, pv.getPosition().subtract(pRef).getNorm(), 3.0e-16 * pRef.getNorm());
688         Assertions.assertEquals(0, pv.getVelocity().subtract(vRef).getNorm(), 2.0e-16 * vRef.getNorm());
689 
690         double[][] jacobian = new double[6][6];
691         orbCir.getJacobianWrtCartesian(PositionAngleType.MEAN, jacobian);
692 
693         for (int i = 0; i < jacobian.length; i++) {
694             double[] row    = jacobian[i];
695             double[] rowRef = jRef[i];
696             for (int j = 0; j < row.length; j++) {
697                 Assertions.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 5.0e-15);
698             }
699         }
700 
701     }
702 
703     @Test
704     void testJacobianFinitedifferences() {
705 
706         AbsoluteDate dateTca = new AbsoluteDate(2000, 4, 1, 0, 0, 0.000, TimeScalesFactory.getUTC());
707         double mu =  3.986004415e+14;
708         CircularOrbit orbCir = new CircularOrbit(7000000.0, 0.01, -0.02, 1.2, 2.1,
709                 0.7, PositionAngleType.MEAN,
710                 FramesFactory.getEME2000(), dateTca, mu);
711 
712         for (PositionAngleType type : PositionAngleType.values()) {
713             double hP = 2.0;
714             double[][] finiteDiffJacobian = finiteDifferencesJacobian(type, orbCir, hP);
715             double[][] jacobian = new double[6][6];
716             orbCir.getJacobianWrtCartesian(type, jacobian);
717 
718             for (int i = 0; i < jacobian.length; i++) {
719                 double[] row    = jacobian[i];
720                 double[] rowRef = finiteDiffJacobian[i];
721                 for (int j = 0; j < row.length; j++) {
722                     Assertions.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 8.0e-9);
723                 }
724             }
725 
726             double[][] invJacobian = new double[6][6];
727             orbCir.getJacobianWrtParameters(type, invJacobian);
728             MatrixUtils.createRealMatrix(jacobian).
729                     multiply(MatrixUtils.createRealMatrix(invJacobian)).
730                     walkInRowOrder(new RealMatrixPreservingVisitor() {
731                         public void start(int rows, int columns,
732                                           int startRow, int endRow, int startColumn, int endColumn) {
733                         }
734 
735                         public void visit(int row, int column, double value) {
736                             Assertions.assertEquals(row == column ? 1.0 : 0.0, value, 4.0e-9);
737                         }
738 
739                         public double end() {
740                             return Double.NaN;
741                         }
742                     });
743 
744         }
745 
746     }
747 
748     private double[][] finiteDifferencesJacobian(PositionAngleType type, CircularOrbit orbit, double hP)
749     {
750         double[][] jacobian = new double[6][6];
751         for (int i = 0; i < 6; ++i) {
752             fillColumn(type, i, orbit, hP, jacobian);
753         }
754         return jacobian;
755     }
756 
757     private void fillColumn(PositionAngleType type, int i, CircularOrbit orbit, double hP, double[][] jacobian) {
758 
759         // at constant energy (i.e. constant semi major axis), we have dV = -mu dP / (V * r^2)
760         // we use this to compute a velocity step size from the position step size
761         Vector3D p = orbit.getPosition();
762         Vector3D v = orbit.getVelocity();
763         double hV = orbit.getMu() * hP / (v.getNorm() * p.getNormSq());
764 
765         double h;
766         Vector3D dP = Vector3D.ZERO;
767         Vector3D dV = Vector3D.ZERO;
768         switch (i) {
769             case 0:
770                 h = hP;
771                 dP = new Vector3D(hP, 0, 0);
772                 break;
773             case 1:
774                 h = hP;
775                 dP = new Vector3D(0, hP, 0);
776                 break;
777             case 2:
778                 h = hP;
779                 dP = new Vector3D(0, 0, hP);
780                 break;
781             case 3:
782                 h = hV;
783                 dV = new Vector3D(hV, 0, 0);
784                 break;
785             case 4:
786                 h = hV;
787                 dV = new Vector3D(0, hV, 0);
788                 break;
789             default:
790                 h = hV;
791                 dV = new Vector3D(0, 0, hV);
792                 break;
793         }
794 
795         CircularOrbit oM4h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, -4, dP), new Vector3D(1, v, -4, dV)),
796                 orbit.getFrame(), orbit.getDate(), orbit.getMu());
797         CircularOrbit oM3h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, -3, dP), new Vector3D(1, v, -3, dV)),
798                 orbit.getFrame(), orbit.getDate(), orbit.getMu());
799         CircularOrbit oM2h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, -2, dP), new Vector3D(1, v, -2, dV)),
800                 orbit.getFrame(), orbit.getDate(), orbit.getMu());
801         CircularOrbit oM1h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, -1, dP), new Vector3D(1, v, -1, dV)),
802                 orbit.getFrame(), orbit.getDate(), orbit.getMu());
803         CircularOrbit oP1h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, +1, dP), new Vector3D(1, v, +1, dV)),
804                 orbit.getFrame(), orbit.getDate(), orbit.getMu());
805         CircularOrbit oP2h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, +2, dP), new Vector3D(1, v, +2, dV)),
806                 orbit.getFrame(), orbit.getDate(), orbit.getMu());
807         CircularOrbit oP3h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, +3, dP), new Vector3D(1, v, +3, dV)),
808                 orbit.getFrame(), orbit.getDate(), orbit.getMu());
809         CircularOrbit oP4h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, +4, dP), new Vector3D(1, v, +4, dV)),
810                 orbit.getFrame(), orbit.getDate(), orbit.getMu());
811 
812         jacobian[0][i] = (-3 * (oP4h.getA()                             - oM4h.getA()) +
813                 32 * (oP3h.getA()                             - oM3h.getA()) -
814                 168 * (oP2h.getA()                             - oM2h.getA()) +
815                 672 * (oP1h.getA()                             - oM1h.getA())) / (840 * h);
816         jacobian[1][i] = (-3 * (oP4h.getCircularEx()                    - oM4h.getCircularEx()) +
817                 32 * (oP3h.getCircularEx()                    - oM3h.getCircularEx()) -
818                 168 * (oP2h.getCircularEx()                    - oM2h.getCircularEx()) +
819                 672 * (oP1h.getCircularEx()                    - oM1h.getCircularEx())) / (840 * h);
820         jacobian[2][i] = (-3 * (oP4h.getCircularEy()                    - oM4h.getCircularEy()) +
821                 32 * (oP3h.getCircularEy()                    - oM3h.getCircularEy()) -
822                 168 * (oP2h.getCircularEy()                    - oM2h.getCircularEy()) +
823                 672 * (oP1h.getCircularEy()                    - oM1h.getCircularEy())) / (840 * h);
824         jacobian[3][i] = (-3 * (oP4h.getI()                             - oM4h.getI()) +
825                 32 * (oP3h.getI()                             - oM3h.getI()) -
826                 168 * (oP2h.getI()                             - oM2h.getI()) +
827                 672 * (oP1h.getI()                             - oM1h.getI())) / (840 * h);
828         jacobian[4][i] = (-3 * (oP4h.getRightAscensionOfAscendingNode() - oM4h.getRightAscensionOfAscendingNode()) +
829                 32 * (oP3h.getRightAscensionOfAscendingNode() - oM3h.getRightAscensionOfAscendingNode()) -
830                 168 * (oP2h.getRightAscensionOfAscendingNode() - oM2h.getRightAscensionOfAscendingNode()) +
831                 672 * (oP1h.getRightAscensionOfAscendingNode() - oM1h.getRightAscensionOfAscendingNode())) / (840 * h);
832         jacobian[5][i] = (-3 * (oP4h.getAlpha(type)                     - oM4h.getAlpha(type)) +
833                 32 * (oP3h.getAlpha(type)                     - oM3h.getAlpha(type)) -
834                 168 * (oP2h.getAlpha(type)                     - oM2h.getAlpha(type)) +
835                 672 * (oP1h.getAlpha(type)                     - oM1h.getAlpha(type))) / (840 * h);
836 
837     }
838 
839     @Test
840     void testNonKeplerianDerivatives() {
841         final AbsoluteDate date         = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
842         final Vector3D     position     = new Vector3D(6896874.444705,  1956581.072644,  -147476.245054);
843         final Vector3D     velocity     = new Vector3D(166.816407662, -1106.783301861, -7372.745712770);
844         final Vector3D     acceleration = new Vector3D(-7.466182457944, -2.118153357345,  0.160004048437);
845         final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity, acceleration);
846         final Frame frame = FramesFactory.getEME2000();
847         final double mu   = Constants.EIGEN5C_EARTH_MU;
848         final CircularOrbit orbit = new CircularOrbit(pv, frame, mu);
849 
850         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getA),
851                 orbit.getADot(),
852                 4.3e-8);
853         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getEquinoctialEx),
854                 orbit.getEquinoctialExDot(),
855                 2.1e-15);
856         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getEquinoctialEy),
857                 orbit.getEquinoctialEyDot(),
858                 5.4e-16);
859         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getHx),
860                 orbit.getHxDot(),
861                 1.6e-15);
862         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getHy),
863                 orbit.getHyDot(),
864                 7.3e-17);
865         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getLv),
866                 orbit.getLvDot(),
867                 3.4e-16);
868         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getLE),
869                 orbit.getLEDot(),
870                 3.5e-15);
871         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getLM),
872                 orbit.getLMDot(),
873                 5.3e-15);
874         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getE),
875                 orbit.getEDot(),
876                 6.8e-16);
877         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getI),
878                 orbit.getIDot(),
879                 5.7e-16);
880         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getCircularEx),
881                 orbit.getCircularExDot(),
882                 2.2e-15);
883         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getCircularEy),
884                 orbit.getCircularEyDot(),
885                 5.3e-17);
886         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getAlphaV),
887                 orbit.getAlphaVDot(),
888                 4.3e-15);
889         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getAlphaE),
890                 orbit.getAlphaEDot(),
891                 1.2e-15);
892         Assertions.assertEquals(differentiate(pv, frame, mu, CircularOrbit::getAlphaM),
893                 orbit.getAlphaMDot(),
894                 3.7e-15);
895         Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAlpha(PositionAngleType.TRUE)),
896                 orbit.getAlphaDot(PositionAngleType.TRUE),
897                 4.3e-15);
898         Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAlpha(PositionAngleType.ECCENTRIC)),
899                 orbit.getAlphaDot(PositionAngleType.ECCENTRIC),
900                 1.2e-15);
901         Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAlpha(PositionAngleType.MEAN)),
902                 orbit.getAlphaDot(PositionAngleType.MEAN),
903                 3.7e-15);
904 
905     }
906 
907     private <S extends Function<CircularOrbit, Double>>
908     double differentiate(TimeStampedPVCoordinates pv, Frame frame, double mu, S picker) {
909         final DSFactory factory = new DSFactory(1, 1);
910         FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(8, 0.1);
911         UnivariateDifferentiableFunction diff =
912                 differentiator.differentiate((UnivariateFunction) dt -> picker.apply(new CircularOrbit(pv.shiftedBy(dt), frame, mu)));
913         return diff.value(factory.variable(0, 0.0)).getPartialDerivative(1);
914     }
915 
916     @Test
917     void testPositionAngleDerivatives() {
918         final AbsoluteDate date         = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
919         final Vector3D     position     = new Vector3D(6896874.444705,  1956581.072644,  -147476.245054);
920         final Vector3D     velocity     = new Vector3D(166.816407662, -1106.783301861, -7372.745712770);
921         final Vector3D     acceleration = new Vector3D(-7.466182457944, -2.118153357345,  0.160004048437);
922         final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity, acceleration);
923         final Frame frame = FramesFactory.getEME2000();
924         final double mu   = Constants.EIGEN5C_EARTH_MU;
925         final CircularOrbit orbit = new CircularOrbit(pv, frame, mu);
926 
927         for (PositionAngleType type : PositionAngleType.values()) {
928             final CircularOrbit rebuilt = new CircularOrbit(orbit.getA(),
929                     orbit.getCircularEx(),
930                     orbit.getCircularEy(),
931                     orbit.getI(),
932                     orbit.getRightAscensionOfAscendingNode(),
933                     orbit.getAlpha(type),
934                     orbit.getADot(),
935                     orbit.getCircularExDot(),
936                     orbit.getCircularEyDot(),
937                     orbit.getIDot(),
938                     orbit.getRightAscensionOfAscendingNodeDot(),
939                     orbit.getAlphaDot(type),
940                     type, orbit.getFrame(), orbit.getDate(), orbit.getMu());
941             MatcherAssert.assertThat(rebuilt.getA(),                                relativelyCloseTo(orbit.getA(),                                1));
942             MatcherAssert.assertThat(rebuilt.getCircularEx(),                       relativelyCloseTo(orbit.getCircularEx(),                       1));
943             MatcherAssert.assertThat(rebuilt.getCircularEy(),                       relativelyCloseTo(orbit.getCircularEy(),                       1));
944             MatcherAssert.assertThat(rebuilt.getE(),                                relativelyCloseTo(orbit.getE(),                                1));
945             MatcherAssert.assertThat(rebuilt.getI(),                                relativelyCloseTo(orbit.getI(),                                1));
946             MatcherAssert.assertThat(rebuilt.getRightAscensionOfAscendingNode(),    relativelyCloseTo(orbit.getRightAscensionOfAscendingNode(),    1));
947             MatcherAssert.assertThat(rebuilt.getADot(),                             relativelyCloseTo(orbit.getADot(),                             1));
948             MatcherAssert.assertThat(rebuilt.getCircularExDot(),                    relativelyCloseTo(orbit.getCircularExDot(),                    1));
949             MatcherAssert.assertThat(rebuilt.getCircularEyDot(),                    relativelyCloseTo(orbit.getCircularEyDot(),                    1));
950             MatcherAssert.assertThat(rebuilt.getEDot(),                             relativelyCloseTo(orbit.getEDot(),                             1));
951             MatcherAssert.assertThat(rebuilt.getIDot(),                             relativelyCloseTo(orbit.getIDot(),                             1));
952             MatcherAssert.assertThat(rebuilt.getRightAscensionOfAscendingNodeDot(), relativelyCloseTo(orbit.getRightAscensionOfAscendingNodeDot(), 1));
953             for (PositionAngleType type2 : PositionAngleType.values()) {
954                 MatcherAssert.assertThat(rebuilt.getAlpha(type2),    relativelyCloseTo(orbit.getAlpha(type2),    1));
955                 MatcherAssert.assertThat(rebuilt.getAlphaDot(type2), relativelyCloseTo(orbit.getAlphaDot(type2), 1));
956             }
957         }
958 
959     }
960 
961     @Test
962     void testEquatorialRetrograde() {
963         Vector3D position = new Vector3D(10000000.0, 0.0, 0.0);
964         Vector3D velocity = new Vector3D(0.0, -6500.0, 1.0e-10);
965         double r2 = position.getNormSq();
966         double r  = FastMath.sqrt(r2);
967         Vector3D acceleration = new Vector3D(-mu / (r * r2), position,
968                 1, new Vector3D(-0.1, 0.2, 0.3));
969         PVCoordinates pvCoordinates = new PVCoordinates(position, velocity, acceleration);
970         CircularOrbit orbit = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
971         Assertions.assertEquals(10637829.465, orbit.getA(), 1.0e-3);
972         Assertions.assertEquals(-738.145, orbit.getADot(), 1.0e-3);
973         Assertions.assertEquals(0.05995861, orbit.getE(), 1.0e-8);
974         Assertions.assertEquals(-6.523e-5, orbit.getEDot(), 1.0e-8);
975         Assertions.assertEquals(FastMath.PI, orbit.getI(), 2.0e-14);
976         Assertions.assertEquals(-4.615e-5, orbit.getIDot(), 1.0e-8);
977         Assertions.assertTrue(Double.isNaN(orbit.getHx()));
978         Assertions.assertTrue(Double.isNaN(orbit.getHxDot()));
979         Assertions.assertTrue(Double.isNaN(orbit.getHy()));
980         Assertions.assertTrue(Double.isNaN(orbit.getHyDot()));
981     }
982 
983     @Test
984     void testDerivativesConversionSymmetry() {
985         final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:01:20.000", TimeScalesFactory.getUTC());
986         Vector3D position     = new Vector3D(6893443.400234382, 1886406.1073757345, -589265.1150359757);
987         Vector3D velocity     = new Vector3D(-281.1261461082365, -1231.6165642450928, -7348.756363469432);
988         Vector3D acceleration = new Vector3D(-7.460341170581685, -2.0415957334584527, 0.6393322823627762);
989         PVCoordinates pvCoordinates = new PVCoordinates( position, velocity, acceleration);
990         CircularOrbit orbit = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(),
991                 date, Constants.EIGEN5C_EARTH_MU);
992         Assertions.assertTrue(orbit.hasNonKeplerianAcceleration());
993         double r2 = position.getNormSq();
994         double r  = FastMath.sqrt(r2);
995         Vector3D keplerianAcceleration = new Vector3D(-orbit.getMu() / (r2 * r), position);
996         Assertions.assertEquals(0.0101, Vector3D.distance(keplerianAcceleration, acceleration), 1.0e-4);
997 
998         for (OrbitType type : OrbitType.values()) {
999             Orbit converted = type.convertType(orbit);
1000             Assertions.assertTrue(converted.hasNonKeplerianAcceleration());
1001             CircularOrbit rebuilt = (CircularOrbit) OrbitType.CIRCULAR.convertType(converted);
1002             Assertions.assertTrue(rebuilt.hasNonKeplerianAcceleration());
1003             Assertions.assertEquals(orbit.getADot(),                             rebuilt.getADot(),                             3.0e-13);
1004             Assertions.assertEquals(orbit.getCircularExDot(),                    rebuilt.getCircularExDot(),                    1.0e-15);
1005             Assertions.assertEquals(orbit.getCircularEyDot(),                    rebuilt.getCircularEyDot(),                    1.0e-15);
1006             Assertions.assertEquals(orbit.getIDot(),                             rebuilt.getIDot(),                             1.0e-15);
1007             Assertions.assertEquals(orbit.getRightAscensionOfAscendingNodeDot(), rebuilt.getRightAscensionOfAscendingNodeDot(), 1.0e-15);
1008             Assertions.assertEquals(orbit.getAlphaVDot(),                        rebuilt.getAlphaVDot(),                        1.0e-15);
1009         }
1010 
1011     }
1012 
1013     @Test
1014     void testToString() {
1015         Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
1016         Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
1017         PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
1018         CircularOrbit orbit = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
1019         Assertions.assertEquals("circular parameters: {a: 4.225517000282565E7, ex: 0.002082917137146049, ey: 5.173980074371024E-4, i: 0.20189257051515358, raan: -87.91788415673473, alphaV: -137.84099636616548;}",
1020                 orbit.toString());
1021     }
1022 
1023     @Test
1024     void testWithKeplerianRates() {
1025         // GIVEN
1026         final Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
1027         final Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
1028         final PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
1029         final CircularOrbit orbit = new CircularOrbit(pvCoordinates, FramesFactory.getGCRF(), date, mu);
1030         // WHEN
1031         final CircularOrbit orbitWithoutRates = orbit.withKeplerianRates();
1032         // THEN
1033         Assertions.assertFalse(orbitWithoutRates.hasNonKeplerianRates());
1034         Assertions.assertEquals(0., orbitWithoutRates.getADot());
1035         Assertions.assertEquals(0., orbitWithoutRates.getCircularExDot());
1036         Assertions.assertEquals(0., orbitWithoutRates.getCircularEyDot());
1037         Assertions.assertEquals(0., orbitWithoutRates.getRightAscensionOfAscendingNodeDot());
1038         Assertions.assertEquals(0., orbitWithoutRates.getIDot());
1039         Assertions.assertEquals(orbit.getMu(), orbitWithoutRates.getMu());
1040         Assertions.assertEquals(orbit.getDate(), orbitWithoutRates.getDate());
1041         Assertions.assertEquals(orbit.getFrame(), orbitWithoutRates.getFrame());
1042         Assertions.assertEquals(orbit.getA(), orbitWithoutRates.getA());
1043         Assertions.assertEquals(orbit.getCircularEx(), orbitWithoutRates.getCircularEx());
1044         Assertions.assertEquals(orbit.getCircularEy(), orbitWithoutRates.getCircularEy());
1045         Assertions.assertEquals(orbit.getRightAscensionOfAscendingNode(),
1046                 orbitWithoutRates.getRightAscensionOfAscendingNode());
1047         Assertions.assertEquals(orbit.getI(), orbitWithoutRates.getI());
1048         Assertions.assertEquals(orbit.getAlphaV(), orbitWithoutRates.getAlphaV());
1049     }
1050 
1051     @Test
1052     void testCopyNonKeplerianAcceleration() {
1053 
1054         final Frame eme2000     = FramesFactory.getEME2000();
1055 
1056         // Define GEO satellite position
1057         final Vector3D position = new Vector3D(42164140, 0, 0);
1058         // Build PVCoodrinates starting from its position and computing the corresponding circular velocity
1059         final PVCoordinates pv  = new PVCoordinates(position,
1060                 new Vector3D(0, FastMath.sqrt(mu / position.getNorm()), 0));
1061         // Build a KeplerianOrbit in eme2000
1062         final Orbit orbit = new CircularOrbit(pv, eme2000, date, mu);
1063 
1064         // Build another KeplerianOrbit as a copy of the first one
1065         final Orbit orbitCopy = new CircularOrbit(orbit);
1066 
1067         // Shift the orbit of a time-interval
1068         final Orbit shiftedOrbit = orbit.shiftedBy(10); // This works good
1069         final Orbit shiftedOrbitCopy = orbitCopy.shiftedBy(10); // This does not work
1070 
1071         Assertions.assertEquals(0.0,
1072                 Vector3D.distance(shiftedOrbit.getPosition(),
1073                         shiftedOrbitCopy.getPosition()),
1074                 1.0e-10);
1075         Assertions.assertEquals(0.0,
1076                 Vector3D.distance(shiftedOrbit.getVelocity(),
1077                         shiftedOrbitCopy.getVelocity()),
1078                 1.0e-10);
1079 
1080     }
1081 
1082     @Test
1083     void testNormalize() {
1084         CircularOrbit withoutDerivatives =
1085                 new CircularOrbit(42166712.0, 0.005, -0.025, 1.6,
1086                         1.25, 0.4, PositionAngleType.MEAN,
1087                         FramesFactory.getEME2000(), date, mu);
1088         CircularOrbit ref =
1089                 new CircularOrbit(24000000.0, -0.012, 0.01, 0.2,
1090                         -6.28, 6.28, PositionAngleType.MEAN,
1091                         FramesFactory.getEME2000(), date, mu);
1092 
1093         CircularOrbit normalized1 = (CircularOrbit) OrbitType.CIRCULAR.normalize(withoutDerivatives, ref);
1094         Assertions.assertFalse(normalized1.hasNonKeplerianAcceleration());
1095         Assertions.assertEquals(0.0, normalized1.getA()          - withoutDerivatives.getA(),          1.0e-6);
1096         Assertions.assertEquals(0.0, normalized1.getCircularEx() - withoutDerivatives.getCircularEx(), 1.0e-10);
1097         Assertions.assertEquals(0.0, normalized1.getCircularEy() - withoutDerivatives.getCircularEy(), 1.0e-10);
1098         Assertions.assertEquals(0.0, normalized1.getI()          - withoutDerivatives.getI(),          1.0e-10);
1099         Assertions.assertEquals(-MathUtils.TWO_PI, normalized1.getRightAscensionOfAscendingNode() - withoutDerivatives.getRightAscensionOfAscendingNode(), 1.0e-10);
1100         Assertions.assertEquals(+MathUtils.TWO_PI, normalized1.getAlphaV() - withoutDerivatives.getAlphaV(), 1.0e-10);
1101         Assertions.assertEquals(withoutDerivatives.getADot(), normalized1.getADot());
1102         Assertions.assertEquals(withoutDerivatives.getCircularExDot(), normalized1.getCircularExDot());
1103         Assertions.assertEquals(withoutDerivatives.getCircularEyDot(), normalized1.getCircularEyDot());
1104         Assertions.assertEquals(withoutDerivatives.getIDot(), normalized1.getIDot());
1105         Assertions.assertEquals(withoutDerivatives.getRightAscensionOfAscendingNodeDot(), normalized1.getRightAscensionOfAscendingNodeDot());
1106         Assertions.assertEquals(withoutDerivatives.getAlphaVDot(), normalized1.getAlphaVDot());
1107 
1108         double[] p = new double[6];
1109         OrbitType.CIRCULAR.mapOrbitToArray(withoutDerivatives, PositionAngleType.TRUE, p, null);
1110         CircularOrbit withDerivatives = (CircularOrbit) OrbitType.CIRCULAR.mapArrayToOrbit(p,
1111                 new double[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 },
1112                 PositionAngleType.TRUE,
1113                 withoutDerivatives.getDate(),
1114                 withoutDerivatives.getMu(),
1115                 withoutDerivatives.getFrame());
1116         CircularOrbit normalized2 = (CircularOrbit) OrbitType.CIRCULAR.normalize(withDerivatives, ref);
1117         Assertions.assertTrue(normalized2.hasNonKeplerianAcceleration());
1118         Assertions.assertEquals(0.0, normalized2.getA()          - withDerivatives.getA(),          1.0e-6);
1119         Assertions.assertEquals(0.0, normalized2.getCircularEx() - withDerivatives.getCircularEx(), 1.0e-10);
1120         Assertions.assertEquals(0.0, normalized2.getCircularEy() - withDerivatives.getCircularEy(), 1.0e-10);
1121         Assertions.assertEquals(0.0, normalized2.getI()          - withDerivatives.getI(),          1.0e-10);
1122         Assertions.assertEquals(-MathUtils.TWO_PI, normalized2.getRightAscensionOfAscendingNode() - withDerivatives.getRightAscensionOfAscendingNode(), 1.0e-10);
1123         Assertions.assertEquals(+MathUtils.TWO_PI, normalized2.getAlphaV() - withDerivatives.getAlphaV(), 1.0e-10);
1124         Assertions.assertEquals(0.0, normalized2.getADot()          - withDerivatives.getADot(),          1.0e-10);
1125         Assertions.assertEquals(0.0, normalized2.getCircularExDot() - withDerivatives.getCircularExDot(), 1.0e-10);
1126         Assertions.assertEquals(0.0, normalized2.getCircularEyDot() - withDerivatives.getCircularEyDot(), 1.0e-10);
1127         Assertions.assertEquals(0.0, normalized2.getIDot()          - withDerivatives.getIDot(),          1.0e-10);
1128         Assertions.assertEquals(0.0, normalized2.getRightAscensionOfAscendingNodeDot() - withDerivatives.getRightAscensionOfAscendingNodeDot(), 1.0e-10);
1129         Assertions.assertEquals(0.0, normalized2.getAlphaVDot() - withDerivatives.getAlphaVDot(), 1.0e-10);
1130 
1131     }
1132 
1133     @Test
1134     void testCoverageCachedPositionAngleTypeWithRates() {
1135         // GIVEN
1136         final double semiMajorAxis = 1e4;
1137         final double eccentricity = 0.;
1138         final double expectedAlpha = 0.;
1139         final double expectedAlphaDot = 0.;
1140         // WHEN & THEN
1141         for (final PositionAngleType inputPositionAngleType: PositionAngleType.values()) {
1142             for (final PositionAngleType cachedPositionAngleType: PositionAngleType.values()) {
1143                 final CircularOrbit circularOrbit = new CircularOrbit(semiMajorAxis, eccentricity, 0., 0., 0.,
1144                         expectedAlpha, 0., 0., 0., 0., 0., expectedAlphaDot,
1145                         inputPositionAngleType, cachedPositionAngleType, FramesFactory.getGCRF(), date, mu);
1146                 Assertions.assertEquals(expectedAlpha, circularOrbit.getAlphaV());
1147                 Assertions.assertEquals(expectedAlpha, circularOrbit.getAlphaM());
1148                 Assertions.assertEquals(expectedAlpha, circularOrbit.getAlphaE());
1149                 Assertions.assertEquals(expectedAlphaDot, circularOrbit.getAlphaVDot());
1150                 Assertions.assertEquals(expectedAlphaDot, circularOrbit.getAlphaMDot());
1151                 Assertions.assertEquals(expectedAlphaDot, circularOrbit.getAlphaEDot());
1152             }
1153         }
1154     }
1155 
1156     @Test
1157     void testCoverageCachedPositionAngleTypeWithoutRates() {
1158         // GIVEN
1159         final double semiMajorAxis = 1e5;
1160         final double eccentricity = 0.;
1161         final double expectedAlpha = 0.;
1162         // WHEN & THEN
1163         for (final PositionAngleType inputPositionAngleType: PositionAngleType.values()) {
1164             for (final PositionAngleType cachedPositionAngleType: PositionAngleType.values()) {
1165                 final CircularOrbit circularOrbit = new CircularOrbit(semiMajorAxis, eccentricity, 0., 0., 0.,
1166                         expectedAlpha, inputPositionAngleType, cachedPositionAngleType,
1167                         FramesFactory.getGCRF(), date, mu);
1168                 Assertions.assertEquals(expectedAlpha, circularOrbit.getAlphaV());
1169                 Assertions.assertEquals(expectedAlpha, circularOrbit.getAlphaM());
1170                 Assertions.assertEquals(expectedAlpha, circularOrbit.getAlphaE());
1171             }
1172         }
1173     }
1174 
1175     @BeforeEach
1176     public void setUp() {
1177 
1178         Utils.setDataRoot("regular-data");
1179 
1180         // Computation date
1181         date = AbsoluteDate.J2000_EPOCH;
1182 
1183         // Body mu
1184         mu = 3.9860047e14;
1185     }
1186 
1187     @AfterEach
1188     public void tearDown() {
1189         date = null;
1190     }
1191 
1192 }