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