1   /* Copyright 2002-2022 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 static org.orekit.OrekitMatchers.relativelyCloseTo;
20  
21  import java.io.ByteArrayInputStream;
22  import java.io.ByteArrayOutputStream;
23  import java.io.IOException;
24  import java.io.ObjectInputStream;
25  import java.io.ObjectOutputStream;
26  import java.util.ArrayList;
27  import java.util.List;
28  import java.util.function.Function;
29  
30  import org.hamcrest.CoreMatchers;
31  import org.hamcrest.MatcherAssert;
32  import org.hipparchus.analysis.UnivariateFunction;
33  import org.hipparchus.analysis.differentiation.DSFactory;
34  import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
35  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
36  import org.hipparchus.geometry.euclidean.threed.Vector3D;
37  import org.hipparchus.linear.MatrixUtils;
38  import org.hipparchus.linear.RealMatrixPreservingVisitor;
39  import org.hipparchus.util.FastMath;
40  import org.hipparchus.util.MathUtils;
41  import org.junit.After;
42  import org.junit.Assert;
43  import org.junit.Before;
44  import org.junit.Test;
45  import org.orekit.Utils;
46  import org.orekit.errors.OrekitException;
47  import org.orekit.errors.OrekitIllegalArgumentException;
48  import org.orekit.errors.OrekitMessages;
49  import org.orekit.frames.Frame;
50  import org.orekit.frames.FramesFactory;
51  import org.orekit.frames.Transform;
52  import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
53  import org.orekit.time.AbsoluteDate;
54  import org.orekit.time.TimeScalesFactory;
55  import org.orekit.utils.Constants;
56  import org.orekit.utils.PVCoordinates;
57  import org.orekit.utils.TimeStampedPVCoordinates;
58  
59  public class KeplerianOrbitTest {
60  
61      // Computation date
62      private AbsoluteDate date;
63  
64     // Body mu
65      private double mu;
66  
67      @Test
68      public void testKeplerianToKeplerian() {
69  
70          // elliptic orbit
71          KeplerianOrbit kep =
72              new KeplerianOrbit(24464560.0, 0.7311, 0.122138, 3.10686, 1.00681,
73                                      0.048363, PositionAngle.MEAN,
74                                      FramesFactory.getEME2000(), date, mu);
75  
76          Vector3D pos = kep.getPVCoordinates().getPosition();
77          Vector3D vit = kep.getPVCoordinates().getVelocity();
78  
79          KeplerianOrbit param = new KeplerianOrbit(new PVCoordinates(pos, vit),
80                                                    FramesFactory.getEME2000(), date, mu);
81          Assert.assertEquals(param.getA(), kep.getA(), Utils.epsilonTest * kep.getA());
82          Assert.assertEquals(param.getE(), kep.getE(), Utils.epsilonE * FastMath.abs(kep.getE()));
83          Assert.assertEquals(MathUtils.normalizeAngle(param.getI(), kep.getI()), kep.getI(), Utils.epsilonAngle * FastMath.abs(kep.getI()));
84          Assert.assertEquals(MathUtils.normalizeAngle(param.getPerigeeArgument(), kep.getPerigeeArgument()), kep.getPerigeeArgument(), Utils.epsilonAngle * FastMath.abs(kep.getPerigeeArgument()));
85          Assert.assertEquals(MathUtils.normalizeAngle(param.getRightAscensionOfAscendingNode(), kep.getRightAscensionOfAscendingNode()), kep.getRightAscensionOfAscendingNode(), Utils.epsilonAngle * FastMath.abs(kep.getRightAscensionOfAscendingNode()));
86          Assert.assertEquals(MathUtils.normalizeAngle(param.getMeanAnomaly(), kep.getMeanAnomaly()), kep.getMeanAnomaly(), Utils.epsilonAngle * FastMath.abs(kep.getMeanAnomaly()));
87  
88          // circular orbit
89          KeplerianOrbit kepCir =
90              new KeplerianOrbit(24464560.0, 0.0, 0.122138, 3.10686, 1.00681,
91                                      0.048363, PositionAngle.MEAN,
92                                      FramesFactory.getEME2000(), date, mu);
93  
94          Vector3D posCir = kepCir.getPVCoordinates().getPosition();
95          Vector3D vitCir = kepCir.getPVCoordinates().getVelocity();
96  
97          KeplerianOrbit paramCir = new KeplerianOrbit(new PVCoordinates(posCir, vitCir),
98                                                       FramesFactory.getEME2000(), date, mu);
99          Assert.assertEquals(paramCir.getA(), kepCir.getA(), Utils.epsilonTest * kepCir.getA());
100         Assert.assertEquals(paramCir.getE(), kepCir.getE(), Utils.epsilonE * FastMath.max(1., FastMath.abs(kepCir.getE())));
101         Assert.assertEquals(MathUtils.normalizeAngle(paramCir.getI(), kepCir.getI()), kepCir.getI(), Utils.epsilonAngle * FastMath.abs(kepCir.getI()));
102         Assert.assertEquals(MathUtils.normalizeAngle(paramCir.getLM(), kepCir.getLM()), kepCir.getLM(), Utils.epsilonAngle * FastMath.abs(kepCir.getLM()));
103         Assert.assertEquals(MathUtils.normalizeAngle(paramCir.getLE(), kepCir.getLE()), kepCir.getLE(), Utils.epsilonAngle * FastMath.abs(kepCir.getLE()));
104         Assert.assertEquals(MathUtils.normalizeAngle(paramCir.getLv(), kepCir.getLv()), kepCir.getLv(), Utils.epsilonAngle * FastMath.abs(kepCir.getLv()));
105 
106         // hyperbolic orbit
107         KeplerianOrbit kepHyp =
108             new KeplerianOrbit(-24464560.0, 1.7311, 0.122138, 3.10686, 1.00681,
109                                     0.048363, PositionAngle.MEAN,
110                                     FramesFactory.getEME2000(), date, mu);
111 
112         Vector3D posHyp = kepHyp.getPVCoordinates().getPosition();
113         Vector3D vitHyp = kepHyp.getPVCoordinates().getVelocity();
114 
115         KeplerianOrbit paramHyp = new KeplerianOrbit(new PVCoordinates(posHyp, vitHyp),
116                                                   FramesFactory.getEME2000(), date, mu);
117         Assert.assertEquals(paramHyp.getA(), kepHyp.getA(), Utils.epsilonTest * FastMath.abs(kepHyp.getA()));
118         Assert.assertEquals(paramHyp.getE(), kepHyp.getE(), Utils.epsilonE * FastMath.abs(kepHyp.getE()));
119         Assert.assertEquals(MathUtils.normalizeAngle(paramHyp.getI(), kepHyp.getI()), kepHyp.getI(), Utils.epsilonAngle * FastMath.abs(kepHyp.getI()));
120         Assert.assertEquals(MathUtils.normalizeAngle(paramHyp.getPerigeeArgument(), kepHyp.getPerigeeArgument()), kepHyp.getPerigeeArgument(), Utils.epsilonAngle * FastMath.abs(kepHyp.getPerigeeArgument()));
121         Assert.assertEquals(MathUtils.normalizeAngle(paramHyp.getRightAscensionOfAscendingNode(), kepHyp.getRightAscensionOfAscendingNode()), kepHyp.getRightAscensionOfAscendingNode(), Utils.epsilonAngle * FastMath.abs(kepHyp.getRightAscensionOfAscendingNode()));
122         Assert.assertEquals(MathUtils.normalizeAngle(paramHyp.getMeanAnomaly(), kepHyp.getMeanAnomaly()), kepHyp.getMeanAnomaly(), Utils.epsilonAngle * FastMath.abs(kepHyp.getMeanAnomaly()));
123 
124     }
125 
126     public void testKeplerianToCartesian() {
127 
128         KeplerianOrbit kep =
129             new KeplerianOrbit(24464560.0, 0.7311, 0.122138, 3.10686, 1.00681,
130                                     0.048363, PositionAngle.MEAN,
131                                     FramesFactory.getEME2000(), date, mu);
132 
133         Vector3D pos = kep.getPVCoordinates().getPosition();
134         Vector3D vit = kep.getPVCoordinates().getVelocity();
135         Assert.assertEquals(-0.107622532467967e+07, pos.getX(), Utils.epsilonTest * FastMath.abs(pos.getX()));
136         Assert.assertEquals(-0.676589636432773e+07, pos.getY(), Utils.epsilonTest * FastMath.abs(pos.getY()));
137         Assert.assertEquals(-0.332308783350379e+06, pos.getZ(), Utils.epsilonTest * FastMath.abs(pos.getZ()));
138 
139         Assert.assertEquals( 0.935685775154103e+04, vit.getX(), Utils.epsilonTest * FastMath.abs(vit.getX()));
140         Assert.assertEquals(-0.331234775037644e+04, vit.getY(), Utils.epsilonTest * FastMath.abs(vit.getY()));
141         Assert.assertEquals(-0.118801577532701e+04, vit.getZ(), Utils.epsilonTest * FastMath.abs(vit.getZ()));
142     }
143 
144     @Test
145     public void testKeplerianToEquinoctial() {
146 
147         KeplerianOrbit kep =
148             new KeplerianOrbit(24464560.0, 0.7311, 0.122138, 3.10686, 1.00681,
149                                     0.048363, PositionAngle.MEAN,
150                                     FramesFactory.getEME2000(), date, mu);
151 
152         Assert.assertEquals(24464560.0, kep.getA(), Utils.epsilonTest * kep.getA());
153         Assert.assertEquals(-0.412036802887626, kep.getEquinoctialEx(), Utils.epsilonE * FastMath.abs(kep.getE()));
154         Assert.assertEquals(-0.603931190671706, kep.getEquinoctialEy(), Utils.epsilonE * FastMath.abs(kep.getE()));
155         Assert.assertEquals(MathUtils.normalizeAngle(2*FastMath.asin(FastMath.sqrt((FastMath.pow(0.652494417368829e-01, 2)+FastMath.pow(0.103158450084864, 2))/4.)), kep.getI()), kep.getI(), Utils.epsilonAngle * FastMath.abs(kep.getI()));
156         Assert.assertEquals(MathUtils.normalizeAngle(0.416203300000000e+01, kep.getLM()), kep.getLM(), Utils.epsilonAngle * FastMath.abs(kep.getLM()));
157 
158     }
159 
160     @Test
161     public void testAnomaly() {
162 
163         Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
164         Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
165         double mu = 3.9860047e14;
166 
167         KeplerianOrbit p = new KeplerianOrbit(new PVCoordinates(position, velocity),
168                                               FramesFactory.getEME2000(), date, mu);
169 
170         // elliptic orbit
171         double e = p.getE();
172         double eRatio = FastMath.sqrt((1 - e) / (1 + e));
173 
174         double v = 1.1;
175         // formulations for elliptic case
176         double E = 2 * FastMath.atan(eRatio * FastMath.tan(v / 2));
177         double M = E - e * FastMath.sin(E);
178 
179         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
180                                     p.getRightAscensionOfAscendingNode(), v , PositionAngle.TRUE,
181                                     p.getFrame(), p.getDate(), p.getMu());
182         Assert.assertEquals(p.getTrueAnomaly(), v, Utils.epsilonAngle * FastMath.abs(v));
183         Assert.assertEquals(p.getEccentricAnomaly(), E, Utils.epsilonAngle * FastMath.abs(E));
184         Assert.assertEquals(p.getMeanAnomaly(), M, Utils.epsilonAngle * FastMath.abs(M));
185         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
186                                     p.getRightAscensionOfAscendingNode(), 0 , PositionAngle.TRUE,
187                                     p.getFrame(), p.getDate(), p.getMu());
188 
189         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
190                                     p.getRightAscensionOfAscendingNode(), E , PositionAngle.ECCENTRIC,
191                                     p.getFrame(), p.getDate(), p.getMu());
192         Assert.assertEquals(p.getTrueAnomaly(), v, Utils.epsilonAngle * FastMath.abs(v));
193         Assert.assertEquals(p.getEccentricAnomaly(), E, Utils.epsilonAngle * FastMath.abs(E));
194         Assert.assertEquals(p.getMeanAnomaly(), M, Utils.epsilonAngle * FastMath.abs(M));
195         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
196                                     p.getRightAscensionOfAscendingNode(), 0 , PositionAngle.TRUE,
197                                     p.getFrame(), p.getDate(), p.getMu());
198 
199         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
200                                     p.getRightAscensionOfAscendingNode(), M, PositionAngle.MEAN,
201                                     p.getFrame(), p.getDate(), p.getMu());
202         Assert.assertEquals(p.getTrueAnomaly(), v, Utils.epsilonAngle * FastMath.abs(v));
203         Assert.assertEquals(p.getEccentricAnomaly(), E, Utils.epsilonAngle * FastMath.abs(E));
204         Assert.assertEquals(p.getMeanAnomaly(), M, Utils.epsilonAngle * FastMath.abs(M));
205 
206         // circular orbit
207         p = new KeplerianOrbit(p.getA(), 0, p.getI(), p.getPerigeeArgument(),
208                                     p.getRightAscensionOfAscendingNode(), p.getLv() , PositionAngle.TRUE,
209                                     p.getFrame(), p.getDate(), p.getMu());
210 
211         E = v;
212         M = E;
213 
214         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
215                                     p.getRightAscensionOfAscendingNode(), v , PositionAngle.TRUE,
216                                     p.getFrame(), p.getDate(), p.getMu());
217         Assert.assertEquals(p.getTrueAnomaly(), v, Utils.epsilonAngle * FastMath.abs(v));
218         Assert.assertEquals(p.getEccentricAnomaly(), E, Utils.epsilonAngle * FastMath.abs(E));
219         Assert.assertEquals(p.getMeanAnomaly(), M, Utils.epsilonAngle * FastMath.abs(M));
220         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
221                                     p.getRightAscensionOfAscendingNode(), 0 , PositionAngle.TRUE,
222                                     p.getFrame(), p.getDate(), p.getMu());
223 
224         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
225                                     p.getRightAscensionOfAscendingNode(), E , PositionAngle.ECCENTRIC, p.getFrame(), p.getDate(), p.getMu());
226         Assert.assertEquals(p.getTrueAnomaly(), v, Utils.epsilonAngle * FastMath.abs(v));
227         Assert.assertEquals(p.getEccentricAnomaly(), E, Utils.epsilonAngle * FastMath.abs(E));
228         Assert.assertEquals(p.getMeanAnomaly(), M, Utils.epsilonAngle * FastMath.abs(M));
229         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
230                                     p.getRightAscensionOfAscendingNode(), 0 , PositionAngle.TRUE,
231                                     p.getFrame(), p.getDate(), p.getMu());
232 
233         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
234                                     p.getRightAscensionOfAscendingNode(), M, PositionAngle.MEAN,
235                                     p.getFrame(), p.getDate(), p.getMu());
236         Assert.assertEquals(p.getTrueAnomaly(), v, Utils.epsilonAngle * FastMath.abs(v));
237         Assert.assertEquals(p.getEccentricAnomaly(), E, Utils.epsilonAngle * FastMath.abs(E));
238         Assert.assertEquals(p.getMeanAnomaly(), M, Utils.epsilonAngle * FastMath.abs(M));
239 
240     }
241 
242     @Test
243     public void testPositionVelocityNorms() {
244         double mu = 3.9860047e14;
245 
246         // elliptic and non equatorial orbit
247         KeplerianOrbit p =
248             new KeplerianOrbit(24464560.0, 0.7311, 2.1, 3.10686, 1.00681,
249                                     0.67, PositionAngle.TRUE,
250                                     FramesFactory.getEME2000(), date, mu);
251 
252         double e       = p.getE();
253         double v       = p.getTrueAnomaly();
254         double ksi     = 1 + e * FastMath.cos(v);
255         double nu      = e * FastMath.sin(v);
256         double epsilon = FastMath.sqrt((1 - e) * (1 + e));
257 
258         double a  = p.getA();
259         double na = FastMath.sqrt(mu / a);
260 
261         // validation of: r = a .(1 - e2) / (1 + e.cos(v))
262         Assert.assertEquals(a * epsilon * epsilon / ksi,
263                      p.getPVCoordinates().getPosition().getNorm(),
264                      Utils.epsilonTest * FastMath.abs(p.getPVCoordinates().getPosition().getNorm()));
265 
266         // validation of: V = sqrt(mu.(1+2e.cos(v)+e2)/a.(1-e2) )
267         Assert.assertEquals(na * FastMath.sqrt(ksi * ksi + nu * nu) / epsilon,
268                      p.getPVCoordinates().getVelocity().getNorm(),
269                      Utils.epsilonTest * FastMath.abs(p.getPVCoordinates().getVelocity().getNorm()));
270 
271 
272         //  circular and equatorial orbit
273         KeplerianOrbit pCirEqua =
274             new KeplerianOrbit(24464560.0, 0.1e-10, 0.1e-8, 3.10686, 1.00681,
275                                     0.67, PositionAngle.TRUE,
276                                     FramesFactory.getEME2000(), date, mu);
277 
278         e       = pCirEqua.getE();
279         v       = pCirEqua.getTrueAnomaly();
280         ksi     = 1 + e * FastMath.cos(v);
281         nu      = e * FastMath.sin(v);
282         epsilon = FastMath.sqrt((1 - e) * (1 + e));
283 
284         a  = pCirEqua.getA();
285         na = FastMath.sqrt(mu / a);
286 
287         // validation of: r = a .(1 - e2) / (1 + e.cos(v))
288         Assert.assertEquals(a * epsilon * epsilon / ksi,
289                      pCirEqua.getPVCoordinates().getPosition().getNorm(),
290                      Utils.epsilonTest * FastMath.abs(pCirEqua.getPVCoordinates().getPosition().getNorm()));
291 
292         // validation of: V = sqrt(mu.(1+2e.cos(v)+e2)/a.(1-e2) )
293         Assert.assertEquals(na * FastMath.sqrt(ksi * ksi + nu * nu) / epsilon,
294                      pCirEqua.getPVCoordinates().getVelocity().getNorm(),
295                      Utils.epsilonTest * FastMath.abs(pCirEqua.getPVCoordinates().getVelocity().getNorm()));
296     }
297 
298     @Test
299     public void testGeometry() {
300         double mu = 3.9860047e14;
301 
302         // elliptic and non equatorial orbit
303         KeplerianOrbit p =
304             new KeplerianOrbit(24464560.0, 0.7311, 2.1, 3.10686, 1.00681,
305                                     0.67, PositionAngle.TRUE,
306                                     FramesFactory.getEME2000(), date, mu);
307 
308         Vector3D position = p.getPVCoordinates().getPosition();
309         Vector3D velocity = p.getPVCoordinates().getVelocity();
310         Vector3D momentum = p.getPVCoordinates().getMomentum().normalize();
311 
312         double apogeeRadius  = p.getA() * (1 + p.getE());
313         double perigeeRadius = p.getA() * (1 - p.getE());
314 
315         for (double lv = 0; lv <= 2 * FastMath.PI; lv += 2 * FastMath.PI/100.) {
316             p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
317                                    p.getRightAscensionOfAscendingNode(), lv , PositionAngle.TRUE,
318                                    p.getFrame(), p.getDate(), p.getMu());
319             position = p.getPVCoordinates().getPosition();
320 
321             // test if the norm of the position is in the range [perigee radius, apogee radius]
322             Assert.assertTrue((position.getNorm() - apogeeRadius)  <= (  apogeeRadius * Utils.epsilonTest));
323             Assert.assertTrue((position.getNorm() - perigeeRadius) >= (- perigeeRadius * Utils.epsilonTest));
324 
325             position = position.normalize();
326             velocity = p.getPVCoordinates().getVelocity();
327             velocity = velocity.normalize();
328 
329             // at this stage of computation, all the vectors (position, velocity and momemtum) are normalized here
330 
331             // test of orthogonality between position and momentum
332             Assert.assertTrue(FastMath.abs(Vector3D.dotProduct(position, momentum)) < Utils.epsilonTest);
333             // test of orthogonality between velocity and momentum
334             Assert.assertTrue(FastMath.abs(Vector3D.dotProduct(velocity, momentum)) < Utils.epsilonTest);
335 
336         }
337 
338         // apsides
339         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
340                                p.getRightAscensionOfAscendingNode(), 0 , PositionAngle.TRUE, p.getFrame(), p.getDate(), p.getMu());
341         Assert.assertEquals(p.getPVCoordinates().getPosition().getNorm(), perigeeRadius, perigeeRadius * Utils.epsilonTest);
342 
343         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
344                                p.getRightAscensionOfAscendingNode(), FastMath.PI , PositionAngle.TRUE, p.getFrame(), p.getDate(), p.getMu());
345         Assert.assertEquals(p.getPVCoordinates().getPosition().getNorm(), apogeeRadius, apogeeRadius * Utils.epsilonTest);
346 
347         // nodes
348         // descending node
349         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
350                                p.getRightAscensionOfAscendingNode(), FastMath.PI - p.getPerigeeArgument() , PositionAngle.TRUE,
351                                p.getFrame(), p.getDate(), p.getMu());
352         Assert.assertTrue(FastMath.abs(p.getPVCoordinates().getPosition().getZ()) < p.getPVCoordinates().getPosition().getNorm() * Utils.epsilonTest);
353         Assert.assertTrue(p.getPVCoordinates().getVelocity().getZ() < 0);
354 
355         // ascending node
356         p = new KeplerianOrbit(p.getA(), p.getE(), p.getI(), p.getPerigeeArgument(),
357                                p.getRightAscensionOfAscendingNode(), 2.0 * FastMath.PI - p.getPerigeeArgument() , PositionAngle.TRUE,
358                                p.getFrame(), p.getDate(), p.getMu());
359         Assert.assertTrue(FastMath.abs(p.getPVCoordinates().getPosition().getZ()) < p.getPVCoordinates().getPosition().getNorm() * Utils.epsilonTest);
360         Assert.assertTrue(p.getPVCoordinates().getVelocity().getZ() > 0);
361 
362 
363         //  circular and equatorial orbit
364         KeplerianOrbit pCirEqua =
365             new KeplerianOrbit(24464560.0, 0.1e-10, 0.1e-8, 3.10686, 1.00681,
366                                0.67, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu);
367 
368         position = pCirEqua.getPVCoordinates().getPosition();
369         velocity = pCirEqua.getPVCoordinates().getVelocity();
370         momentum = Vector3D.crossProduct(position, velocity).normalize();
371 
372         apogeeRadius  = pCirEqua.getA() * (1 + pCirEqua.getE());
373         perigeeRadius = pCirEqua.getA() * (1 - pCirEqua.getE());
374         // test if apogee equals perigee
375         Assert.assertEquals(perigeeRadius, apogeeRadius, 1.e+4 * Utils.epsilonTest * apogeeRadius);
376 
377         for (double lv = 0; lv <= 2 * FastMath.PI; lv += 2 * FastMath.PI/100.) {
378             pCirEqua = new KeplerianOrbit(pCirEqua.getA(), pCirEqua.getE(), pCirEqua.getI(), pCirEqua.getPerigeeArgument(),
379                                           pCirEqua.getRightAscensionOfAscendingNode(), lv, PositionAngle.TRUE,
380                                           pCirEqua.getFrame(), pCirEqua.getDate(), pCirEqua.getMu());
381             position = pCirEqua.getPVCoordinates().getPosition();
382 
383             // test if the norm pf the position is in the range [perigee radius, apogee radius]
384             // Warning: these tests are without absolute value by choice
385             Assert.assertTrue((position.getNorm() - apogeeRadius)  <= (  apogeeRadius * Utils.epsilonTest));
386             Assert.assertTrue((position.getNorm() - perigeeRadius) >= (- perigeeRadius * Utils.epsilonTest));
387 
388             position = position.normalize();
389             velocity = pCirEqua.getPVCoordinates().getVelocity();
390             velocity = velocity.normalize();
391 
392             // at this stage of computation, all the vectors (position, velocity and momemtum) are normalized here
393 
394             // test of orthogonality between position and momentum
395             Assert.assertTrue(FastMath.abs(Vector3D.dotProduct(position, momentum)) < Utils.epsilonTest);
396             // test of orthogonality between velocity and momentum
397             Assert.assertTrue(FastMath.abs(Vector3D.dotProduct(velocity, momentum)) < Utils.epsilonTest);
398 
399         }
400     }
401 
402     @Test
403     public void testSymmetry() {
404 
405         // elliptic and non equatorial orbit
406         Vector3D position = new Vector3D(-4947831., -3765382., -3708221.);
407         Vector3D velocity = new Vector3D(-2079., 5291., -7842.);
408         double mu = 3.9860047e14;
409 
410         KeplerianOrbit p = new KeplerianOrbit(new PVCoordinates(position, velocity),
411                                               FramesFactory.getEME2000(), date, mu);
412         Vector3D positionOffset = p.getPVCoordinates().getPosition().subtract(position);
413         Vector3D velocityOffset = p.getPVCoordinates().getVelocity().subtract(velocity);
414 
415         Assert.assertTrue(positionOffset.getNorm() < Utils.epsilonTest);
416         Assert.assertTrue(velocityOffset.getNorm() < Utils.epsilonTest);
417 
418         // circular and equatorial orbit
419         position = new Vector3D(1742382., -2.440243e7, -0.014517);
420         velocity = new Vector3D(4026.2, 287.479, -3.e-6);
421 
422 
423         p = new KeplerianOrbit(new PVCoordinates(position, velocity),
424                                FramesFactory.getEME2000(), date, mu);
425         positionOffset = p.getPVCoordinates().getPosition().subtract(position);
426         velocityOffset = p.getPVCoordinates().getVelocity().subtract(velocity);
427 
428         Assert.assertTrue(positionOffset.getNorm() < Utils.epsilonTest);
429         Assert.assertTrue(velocityOffset.getNorm() < Utils.epsilonTest);
430 
431     }
432 
433     @Test(expected=IllegalArgumentException.class)
434     public void testNonInertialFrame() throws IllegalArgumentException {
435 
436         Vector3D position = new Vector3D(-4947831., -3765382., -3708221.);
437         Vector3D velocity = new Vector3D(-2079., 5291., -7842.);
438         PVCoordinates pvCoordinates = new PVCoordinates( position, velocity);
439         new KeplerianOrbit(pvCoordinates,
440                            new Frame(FramesFactory.getEME2000(), Transform.IDENTITY, "non-inertial", false),
441                            date, mu);
442     }
443 
444     @Test
445     public void testPeriod() {
446         KeplerianOrbit orbit = new KeplerianOrbit(7654321.0, 0.1, 0.2, 0, 0, 0,
447                                                   PositionAngle.TRUE,
448                                                   FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
449                                                   mu);
450         Assert.assertEquals(6664.5521723383589487, orbit.getKeplerianPeriod(), 1.0e-12);
451         Assert.assertEquals(0.00094277682051291315229, orbit.getKeplerianMeanMotion(), 1.0e-16);
452     }
453 
454     @Test
455     public void testHyperbola1() {
456         KeplerianOrbit orbit = new KeplerianOrbit(-10000000.0, 2.5, 0.3, 0, 0, 0.0,
457                                                   PositionAngle.TRUE,
458                                                   FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
459                                                   mu);
460         Vector3D perigeeP  = orbit.getPVCoordinates().getPosition();
461         Vector3D u = perigeeP.normalize();
462         Vector3D focus1 = Vector3D.ZERO;
463         Vector3D focus2 = new Vector3D(-2 * orbit.getA() * orbit.getE(), u);
464         for (double dt = -5000; dt < 5000; dt += 60) {
465             PVCoordinates pv = orbit.shiftedBy(dt).getPVCoordinates();
466             double d1 = Vector3D.distance(pv.getPosition(), focus1);
467             double d2 = Vector3D.distance(pv.getPosition(), focus2);
468             Assert.assertEquals(-2 * orbit.getA(), FastMath.abs(d1 - d2), 1.0e-6);
469             KeplerianOrbit rebuilt =
470                 new KeplerianOrbit(pv, orbit.getFrame(), orbit.getDate().shiftedBy(dt), mu);
471             Assert.assertEquals(-10000000.0, rebuilt.getA(), 1.0e-6);
472             Assert.assertEquals(2.5, rebuilt.getE(), 1.0e-13);
473         }
474     }
475 
476     @Test
477     public void testHyperbola2() {
478         KeplerianOrbit orbit = new KeplerianOrbit(-10000000.0, 1.2, 0.3, 0, 0, -1.75,
479                                                   PositionAngle.MEAN,
480                                                   FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
481                                                   mu);
482         Vector3D perigeeP  = new KeplerianOrbit(orbit.getA(), orbit.getE(), orbit.getI(),
483                                                 orbit.getPerigeeArgument(), orbit.getRightAscensionOfAscendingNode(),
484                                                 0.0, PositionAngle.TRUE, orbit.getFrame(),
485                                                 orbit.getDate(), orbit.getMu()).getPVCoordinates().getPosition();
486         Vector3D u = perigeeP.normalize();
487         Vector3D focus1 = Vector3D.ZERO;
488         Vector3D focus2 = new Vector3D(-2 * orbit.getA() * orbit.getE(), u);
489         for (double dt = -5000; dt < 5000; dt += 60) {
490             PVCoordinates pv = orbit.shiftedBy(dt).getPVCoordinates();
491             double d1 = Vector3D.distance(pv.getPosition(), focus1);
492             double d2 = Vector3D.distance(pv.getPosition(), focus2);
493             Assert.assertEquals(-2 * orbit.getA(), FastMath.abs(d1 - d2), 1.0e-6);
494             KeplerianOrbit rebuilt =
495                 new KeplerianOrbit(pv, orbit.getFrame(), orbit.getDate().shiftedBy(dt), mu);
496             Assert.assertEquals(-10000000.0, rebuilt.getA(), 1.0e-6);
497             Assert.assertEquals(1.2, rebuilt.getE(), 1.0e-13);
498         }
499     }
500 
501     @Test
502     public void testInconsistentHyperbola() {
503         try {
504             new KeplerianOrbit(+10000000.0, 2.5, 0.3, 0, 0, 0.0,
505                                                   PositionAngle.TRUE,
506                                                   FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
507                                                   mu);
508             Assert.fail("an exception should have been thrown");
509         } catch (OrekitIllegalArgumentException oe) {
510             Assert.assertEquals(OrekitMessages.ORBIT_A_E_MISMATCH_WITH_CONIC_TYPE, oe.getSpecifier());
511             Assert.assertEquals(+10000000.0, ((Double) oe.getParts()[0]).doubleValue(), 1.0e-3);
512             Assert.assertEquals(2.5,         ((Double) oe.getParts()[1]).doubleValue(), 1.0e-15);
513         }
514     }
515 
516     @Test
517     public void testVeryLargeEccentricity() {
518 
519         final Frame eme2000 = FramesFactory.getEME2000();
520         final double meanAnomaly = 1.;
521         final KeplerianOrbit orb0 = new KeplerianOrbit(42600e3, 0.9, 0.00001, 0, 0,
522                                                        FastMath.toRadians(meanAnomaly),
523                                                        PositionAngle.MEAN, eme2000, date, mu);
524 
525         // big dV along Y
526         final Vector3D deltaV = new Vector3D(0.0, 110000.0, 0.0);
527         final PVCoordinates pv1 = new PVCoordinates(orb0.getPVCoordinates().getPosition(),
528                                                     orb0.getPVCoordinates().getVelocity().add(deltaV));
529         final KeplerianOrbit orb1 = new KeplerianOrbit(pv1, eme2000, date, mu);
530 
531         // Despite large eccentricity, the conversion of mean anomaly to hyperbolic eccentric anomaly
532         // converges in less than 50 iterations (issue #114)
533         final PVCoordinates pvTested    = orb1.shiftedBy(0).getPVCoordinates();
534         final Vector3D      pTested     = pvTested.getPosition();
535         final Vector3D      vTested     = pvTested.getVelocity();
536 
537         final PVCoordinates pvReference = orb1.getPVCoordinates();
538         final Vector3D      pReference  = pvReference.getPosition();
539         final Vector3D      vReference  = pvReference.getVelocity();
540 
541         final double threshold = 1.e-15;
542         Assert.assertEquals(0, pTested.subtract(pReference).getNorm(), threshold * pReference.getNorm());
543         Assert.assertEquals(0, vTested.subtract(vReference).getNorm(), threshold * vReference.getNorm());
544 
545     }
546 
547     @Test
548     public void testKeplerEquation() {
549 
550         for (double M = -6 * FastMath.PI; M < 6 * FastMath.PI; M += 0.01) {
551             KeplerianOrbit pElliptic =
552                 new KeplerianOrbit(24464560.0, 0.7311, 2.1, 3.10686, 1.00681,
553                                    M, PositionAngle.MEAN,
554                                    FramesFactory.getEME2000(), date, mu);
555             double E = pElliptic.getEccentricAnomaly();
556             double e = pElliptic.getE();
557             Assert.assertEquals(M, E - e * FastMath.sin(E), 2.0e-14);
558         }
559 
560         for (double M = -6 * FastMath.PI; M < 6 * FastMath.PI; M += 0.01) {
561             KeplerianOrbit pAlmostParabolic =
562                 new KeplerianOrbit(24464560.0, 0.9999, 2.1, 3.10686, 1.00681,
563                                    M, PositionAngle.MEAN,
564                                    FramesFactory.getEME2000(), date, mu);
565             double E = pAlmostParabolic.getEccentricAnomaly();
566             double e = pAlmostParabolic.getE();
567             Assert.assertEquals(M, E - e * FastMath.sin(E), 4.0e-13);
568         }
569 
570     }
571 
572     @Test(expected=IllegalArgumentException.class)
573     public void testOutOfRangeV() {
574         new KeplerianOrbit(-7000434.460140012, 1.1999785407363386, 1.3962787004479158,
575                            1.3962320168955138, 0.3490728321331678, -2.55593407037698,
576                            PositionAngle.TRUE, FramesFactory.getEME2000(),
577                            new AbsoluteDate("2000-01-01T12:00:00.391", TimeScalesFactory.getUTC()),
578                            3.986004415E14);
579     }
580 
581     @Test
582     public void testNumericalIssue25() {
583         Vector3D position = new Vector3D(3782116.14107698, 416663.11924914, 5875541.62103057);
584         Vector3D velocity = new Vector3D(-6349.7848910501, 288.4061811651, 4066.9366759691);
585         KeplerianOrbit orbit = new KeplerianOrbit(new PVCoordinates(position, velocity),
586                                                   FramesFactory.getEME2000(),
587                                                   new AbsoluteDate("2004-01-01T23:00:00.000",
588                                                                    TimeScalesFactory.getUTC()),
589                                                   3.986004415E14);
590         Assert.assertEquals(0.0, orbit.getE(), 2.0e-14);
591     }
592 
593     @Test
594     public void testPerfectlyEquatorial() {
595         Vector3D position = new Vector3D(6957904.3624652653594, 766529.11411558074507, 0);
596         Vector3D velocity = new Vector3D(-7538.2817012412102845, 342.38751001881413381, 0.);
597         KeplerianOrbit orbit = new KeplerianOrbit(new PVCoordinates(position, velocity),
598                                                   FramesFactory.getEME2000(),
599                                                   new AbsoluteDate("2004-01-01T23:00:00.000",
600                                                                    TimeScalesFactory.getUTC()),
601                                                   3.986004415E14);
602         Assert.assertEquals(0.0, orbit.getI(), 2.0e-14);
603         Assert.assertEquals(0.0, orbit.getRightAscensionOfAscendingNode(), 2.0e-14);
604     }
605 
606     @Test
607     public void testJacobianReferenceEllipse() {
608 
609         AbsoluteDate dateTca = new AbsoluteDate(2000, 04, 01, 0, 0, 0.000, TimeScalesFactory.getUTC());
610         double mu =  3.986004415e+14;
611         KeplerianOrbit orbKep = new KeplerianOrbit(7000000.0, 0.01, FastMath.toRadians(80.), FastMath.toRadians(80.), FastMath.toRadians(20.),
612                                           FastMath.toRadians(40.), PositionAngle.MEAN,
613                                           FramesFactory.getEME2000(), dateTca, mu);
614 
615         // the following reference values have been computed using the free software
616         // version 6.2 of the MSLIB fortran library by the following program:
617         //        program kep_jacobian
618         //
619         //        use mslib
620         //        implicit none
621         //
622         //        integer, parameter :: nb = 11
623         //        integer :: i,j
624         //        type(tm_code_retour)      ::  code_retour
625         //
626         //        real(pm_reel), parameter :: mu= 3.986004415e+14_pm_reel
627         //        real(pm_reel),dimension(3)::vit_car,pos_car
628         //        type(tm_orb_kep)::kep
629         //        real(pm_reel), dimension(6,6)::jacob
630         //        real(pm_reel)::norme
631         //
632         //        kep%a=7000000_pm_reel
633         //        kep%e=0.01_pm_reel
634         //        kep%i=80_pm_reel*pm_deg_rad
635         //        kep%pom=80_pm_reel*pm_deg_rad
636         //        kep%gom=20_pm_reel*pm_deg_rad
637         //        kep%M=40_pm_reel*pm_deg_rad
638         //
639         //        call mv_kep_car(mu,kep,pos_car,vit_car,code_retour)
640         //        write(*,*)code_retour%valeur
641         //        write(*,1000)pos_car,vit_car
642         //
643         //
644         //        call mu_norme(pos_car,norme,code_retour)
645         //        write(*,*)norme
646         //
647         //        call mv_car_kep (mu, pos_car, vit_car, kep, code_retour, jacob)
648         //        write(*,*)code_retour%valeur
649         //
650         //        write(*,*)"kep = ", kep%a, kep%e, kep%i*pm_rad_deg,&
651         //                            kep%pom*pm_rad_deg, kep%gom*pm_rad_deg, kep%M*pm_rad_deg
652         //
653         //        do i = 1,6
654         //           write(*,*) " ",(jacob(i,j),j=1,6)
655         //        end do
656         //
657         //        1000 format (6(f24.15,1x))
658         //        end program kep_jacobian
659         Vector3D pRef = new Vector3D(-3691555.569874833337963, -240330.253992714860942, 5879700.285850423388183);
660         Vector3D vRef = new Vector3D(-5936.229884450408463, -2871.067660163344044, -3786.209549192726627);
661         double[][] jRef = {
662             { -1.0792090588217809,       -7.02594292049818631E-002,  1.7189029642216496,       -1459.4829009393857,       -705.88138246206040,       -930.87838644776593       },
663             { -1.31195762636625214E-007, -3.90087231593959271E-008,  4.65917592901869866E-008, -2.02467187867647177E-004, -7.89767994436215424E-005, -2.81639203329454407E-005 },
664             {  4.18334478744371316E-008, -1.14936453412947957E-007,  2.15670500707930151E-008, -2.26450325965329431E-005,  6.22167157217876380E-005, -1.16745469637130306E-005 },
665             {  3.52735168061691945E-006,  3.82555734454450974E-006,  1.34715077236557634E-005, -8.06586262922115264E-003, -6.13725651685311825E-003, -1.71765290503914092E-002 },
666             {  2.48948022169790885E-008, -6.83979069529389238E-008,  1.28344057971888544E-008,  3.86597661353874888E-005, -1.06216834498373629E-004,  1.99308724078785540E-005 },
667             { -3.41911705254704525E-006, -3.75913623359912437E-006, -1.34013845492518465E-005,  8.19851888816422458E-003,  6.16449264680494959E-003,  1.69495878276556648E-002 }
668         };
669 
670         PVCoordinates pv = orbKep.getPVCoordinates();
671         Assert.assertEquals(0, pv.getPosition().subtract(pRef).getNorm(), 1.0e-15 * pRef.getNorm());
672         Assert.assertEquals(0, pv.getVelocity().subtract(vRef).getNorm(), 1.0e-16 * vRef.getNorm());
673 
674         double[][] jacobian = new double[6][6];
675         orbKep.getJacobianWrtCartesian(PositionAngle.MEAN, jacobian);
676 
677         for (int i = 0; i < jacobian.length; i++) {
678             double[] row    = jacobian[i];
679             double[] rowRef = jRef[i];
680             for (int j = 0; j < row.length; j++) {
681                 Assert.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 2.0e-12);
682             }
683         }
684 
685     }
686 
687     @Test
688     public void testJacobianFinitedifferencesEllipse() {
689 
690         AbsoluteDate dateTca = new AbsoluteDate(2000, 04, 01, 0, 0, 0.000, TimeScalesFactory.getUTC());
691         double mu =  3.986004415e+14;
692         KeplerianOrbit orbKep = new KeplerianOrbit(7000000.0, 0.01, FastMath.toRadians(80.), FastMath.toRadians(80.), FastMath.toRadians(20.),
693                                           FastMath.toRadians(40.), PositionAngle.MEAN,
694                                           FramesFactory.getEME2000(), dateTca, mu);
695 
696         for (PositionAngle type : PositionAngle.values()) {
697             double hP = 2.0;
698             double[][] finiteDiffJacobian = finiteDifferencesJacobian(type, orbKep, hP);
699             double[][] jacobian = new double[6][6];
700             orbKep.getJacobianWrtCartesian(type, jacobian);
701 
702             for (int i = 0; i < jacobian.length; i++) {
703                 double[] row    = jacobian[i];
704                 double[] rowRef = finiteDiffJacobian[i];
705                 for (int j = 0; j < row.length; j++) {
706                     Assert.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 2.0e-7);
707                 }
708             }
709 
710             double[][] invJacobian = new double[6][6];
711             orbKep.getJacobianWrtParameters(type, invJacobian);
712             MatrixUtils.createRealMatrix(jacobian).
713                             multiply(MatrixUtils.createRealMatrix(invJacobian)).
714             walkInRowOrder(new RealMatrixPreservingVisitor() {
715                 public void start(int rows, int columns,
716                                   int startRow, int endRow, int startColumn, int endColumn) {
717                 }
718 
719                 public void visit(int row, int column, double value) {
720                     Assert.assertEquals(row == column ? 1.0 : 0.0, value, 5.0e-9);
721                 }
722 
723                 public double end() {
724                     return Double.NaN;
725                 }
726             });
727 
728         }
729 
730     }
731 
732     @Test
733     public void testJacobianReferenceHyperbola() {
734 
735         AbsoluteDate dateTca = new AbsoluteDate(2000, 04, 01, 0, 0, 0.000, TimeScalesFactory.getUTC());
736         double mu =  3.986004415e+14;
737         KeplerianOrbit orbKep = new KeplerianOrbit(-7000000.0, 1.2, FastMath.toRadians(80.), FastMath.toRadians(80.), FastMath.toRadians(20.),
738                                           FastMath.toRadians(40.), PositionAngle.MEAN,
739                                           FramesFactory.getEME2000(), dateTca, mu);
740 
741         // the following reference values have been computed using the free software
742         // version 6.2 of the MSLIB fortran library by the following program:
743         //        program kep_hyperb_jacobian
744         //
745         //        use mslib
746         //        implicit none
747         //
748         //        integer, parameter :: nb = 11
749         //        integer :: i,j
750         //        type(tm_code_retour)      ::  code_retour
751         //
752         //        real(pm_reel), parameter :: mu= 3.986004415e+14_pm_reel
753         //        real(pm_reel),dimension(3)::vit_car,pos_car
754         //        type(tm_orb_kep)::kep
755         //        real(pm_reel), dimension(6,6)::jacob
756         //        real(pm_reel)::norme
757         //
758         //        kep%a=7000000_pm_reel
759         //        kep%e=1.2_pm_reel
760         //        kep%i=80_pm_reel*pm_deg_rad
761         //        kep%pom=80_pm_reel*pm_deg_rad
762         //        kep%gom=20_pm_reel*pm_deg_rad
763         //        kep%M=40_pm_reel*pm_deg_rad
764         //
765         //        call mv_kep_car(mu,kep,pos_car,vit_car,code_retour)
766         //        write(*,*)code_retour%valeur
767         //        write(*,1000)pos_car,vit_car
768         //
769         //
770         //        call mu_norme(pos_car,norme,code_retour)
771         //        write(*,*)norme
772         //
773         //        call mv_car_kep (mu, pos_car, vit_car, kep, code_retour, jacob)
774         //        write(*,*)code_retour%valeur
775         //
776         //        write(*,*)"kep = ", kep%a, kep%e, kep%i*pm_rad_deg,&
777         //                            kep%pom*pm_rad_deg, kep%gom*pm_rad_deg, kep%M*pm_rad_deg
778         //
779         //        ! convert the sign of da row since mslib uses a > 0 for all orbits
780         //        ! whereas we use a < 0 for hyperbolic orbits
781         //        write(*,*) " ",(-jacob(1,j),j=1,6)
782         //        do i = 2,6
783         //           write(*,*) " ",(jacob(i,j),j=1,6)
784         //        end do
785         //
786         //        1000 format (6(f24.15,1x))
787         //        end program kep_hyperb_jacobian
788         Vector3D pRef = new Vector3D(-7654711.206549182534218, -3460171.872979687992483, -3592374.514463655184954);
789         Vector3D vRef = new Vector3D(   -7886.368091820805603,    -4359.739012331759113,    -7937.060044548694350);
790         double[][] jRef = {
791             {  -0.98364725131848019,      -0.44463970750901238,      -0.46162803814668391,       -1938.9443476028839,       -1071.8864775981751,       -1951.4074832397598      },
792             {  -1.10548813242982574E-007, -2.52906747183730431E-008,  7.96500937398593591E-008, -9.70479823470940108E-006, -2.93209076428001017E-005, -1.37434463892791042E-004 },
793             {   8.55737680891616672E-008, -2.35111995522618220E-007,  4.41171797903162743E-008, -8.05235180390949802E-005,  2.21236547547460423E-004, -4.15135455876865407E-005 },
794             {  -1.52641427784095578E-007,  1.10250447958827901E-008,  1.21265251605359894E-007,  7.63347077200903542E-005, -3.54738331412232378E-005, -2.31400737283033359E-004 },
795             {   7.86711766048035274E-008, -2.16147281283624453E-007,  4.05585791077187359E-008, -3.56071805267582894E-005,  9.78299244677127374E-005, -1.83571253224293247E-005 },
796             {  -2.41488884881911384E-007, -1.00119615610276537E-007, -6.51494225096757969E-008, -2.43295075073248163E-004, -1.43273725071890463E-004, -2.91625510452094873E-004 }
797         };
798 
799         PVCoordinates pv = orbKep.getPVCoordinates();
800         Assert.assertEquals(0, pv.getPosition().subtract(pRef).getNorm() / pRef.getNorm(), 1.0e-16);
801 //        Assert.assertEquals(0, pv.getPosition().subtract(pRef).getNorm() / pRef.getNorm(), 2.0e-15);
802         Assert.assertEquals(0, pv.getVelocity().subtract(vRef).getNorm() / vRef.getNorm(), 3.0e-16);
803 
804         double[][] jacobian = new double[6][6];
805         orbKep.getJacobianWrtCartesian(PositionAngle.MEAN, jacobian);
806 
807         for (int i = 0; i < jacobian.length; i++) {
808             double[] row    = jacobian[i];
809             double[] rowRef = jRef[i];
810             for (int j = 0; j < row.length; j++) {
811                 Assert.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 1.0e-14);
812             }
813         }
814 
815     }
816 
817     @Test
818     public void testJacobianFinitedifferencesHyperbola() {
819 
820         AbsoluteDate dateTca = new AbsoluteDate(2000, 04, 01, 0, 0, 0.000, TimeScalesFactory.getUTC());
821         double mu =  3.986004415e+14;
822         KeplerianOrbit orbKep = new KeplerianOrbit(-7000000.0, 1.2, FastMath.toRadians(80.), FastMath.toRadians(80.), FastMath.toRadians(20.),
823                                                    FastMath.toRadians(40.), PositionAngle.MEAN,
824                                                    FramesFactory.getEME2000(), dateTca, mu);
825 
826         for (PositionAngle type : PositionAngle.values()) {
827             double hP = 2.0;
828             double[][] finiteDiffJacobian = finiteDifferencesJacobian(type, orbKep, hP);
829             double[][] jacobian = new double[6][6];
830             orbKep.getJacobianWrtCartesian(type, jacobian);
831             for (int i = 0; i < jacobian.length; i++) {
832                 double[] row    = jacobian[i];
833                 double[] rowRef = finiteDiffJacobian[i];
834                 for (int j = 0; j < row.length; j++) {
835                     Assert.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 3.0e-8);
836                 }
837             }
838 
839             double[][] invJacobian = new double[6][6];
840             orbKep.getJacobianWrtParameters(type, invJacobian);
841             MatrixUtils.createRealMatrix(jacobian).
842                             multiply(MatrixUtils.createRealMatrix(invJacobian)).
843             walkInRowOrder(new RealMatrixPreservingVisitor() {
844                 public void start(int rows, int columns,
845                                   int startRow, int endRow, int startColumn, int endColumn) {
846                 }
847 
848                 public void visit(int row, int column, double value) {
849                     Assert.assertEquals(row == column ? 1.0 : 0.0, value, 2.0e-8);
850                 }
851 
852                 public double end() {
853                     return Double.NaN;
854                 }
855             });
856 
857         }
858 
859     }
860 
861     private double[][] finiteDifferencesJacobian(PositionAngle type, KeplerianOrbit orbit, double hP)
862         {
863         double[][] jacobian = new double[6][6];
864         for (int i = 0; i < 6; ++i) {
865             fillColumn(type, i, orbit, hP, jacobian);
866         }
867         return jacobian;
868     }
869 
870     private void fillColumn(PositionAngle type, int i, KeplerianOrbit orbit, double hP, double[][] jacobian) {
871 
872         // at constant energy (i.e. constant semi major axis), we have dV = -mu dP / (V * r^2)
873         // we use this to compute a velocity step size from the position step size
874         Vector3D p = orbit.getPVCoordinates().getPosition();
875         Vector3D v = orbit.getPVCoordinates().getVelocity();
876         double hV = orbit.getMu() * hP / (v.getNorm() * p.getNormSq());
877 
878         double h;
879         Vector3D dP = Vector3D.ZERO;
880         Vector3D dV = Vector3D.ZERO;
881         switch (i) {
882         case 0:
883             h = hP;
884             dP = new Vector3D(hP, 0, 0);
885             break;
886         case 1:
887             h = hP;
888             dP = new Vector3D(0, hP, 0);
889             break;
890         case 2:
891             h = hP;
892             dP = new Vector3D(0, 0, hP);
893             break;
894         case 3:
895             h = hV;
896             dV = new Vector3D(hV, 0, 0);
897             break;
898         case 4:
899             h = hV;
900             dV = new Vector3D(0, hV, 0);
901             break;
902         default:
903             h = hV;
904             dV = new Vector3D(0, 0, hV);
905             break;
906         }
907 
908         KeplerianOrbit oM4h = new KeplerianOrbit(new PVCoordinates(new Vector3D(1, p, -4, dP), new Vector3D(1, v, -4, dV)),
909                                                  orbit.getFrame(), orbit.getDate(), orbit.getMu());
910         KeplerianOrbit oM3h = new KeplerianOrbit(new PVCoordinates(new Vector3D(1, p, -3, dP), new Vector3D(1, v, -3, dV)),
911                                                  orbit.getFrame(), orbit.getDate(), orbit.getMu());
912         KeplerianOrbit oM2h = new KeplerianOrbit(new PVCoordinates(new Vector3D(1, p, -2, dP), new Vector3D(1, v, -2, dV)),
913                                                  orbit.getFrame(), orbit.getDate(), orbit.getMu());
914         KeplerianOrbit oM1h = new KeplerianOrbit(new PVCoordinates(new Vector3D(1, p, -1, dP), new Vector3D(1, v, -1, dV)),
915                                                  orbit.getFrame(), orbit.getDate(), orbit.getMu());
916         KeplerianOrbit oP1h = new KeplerianOrbit(new PVCoordinates(new Vector3D(1, p, +1, dP), new Vector3D(1, v, +1, dV)),
917                                                  orbit.getFrame(), orbit.getDate(), orbit.getMu());
918         KeplerianOrbit oP2h = new KeplerianOrbit(new PVCoordinates(new Vector3D(1, p, +2, dP), new Vector3D(1, v, +2, dV)),
919                                                  orbit.getFrame(), orbit.getDate(), orbit.getMu());
920         KeplerianOrbit oP3h = new KeplerianOrbit(new PVCoordinates(new Vector3D(1, p, +3, dP), new Vector3D(1, v, +3, dV)),
921                                                  orbit.getFrame(), orbit.getDate(), orbit.getMu());
922         KeplerianOrbit oP4h = new KeplerianOrbit(new PVCoordinates(new Vector3D(1, p, +4, dP), new Vector3D(1, v, +4, dV)),
923                                                  orbit.getFrame(), orbit.getDate(), orbit.getMu());
924 
925         jacobian[0][i] = (-3 * (oP4h.getA()                             - oM4h.getA()) +
926                           32 * (oP3h.getA()                             - oM3h.getA()) -
927                          168 * (oP2h.getA()                             - oM2h.getA()) +
928                          672 * (oP1h.getA()                             - oM1h.getA())) / (840 * h);
929         jacobian[1][i] = (-3 * (oP4h.getE()                             - oM4h.getE()) +
930                           32 * (oP3h.getE()                             - oM3h.getE()) -
931                          168 * (oP2h.getE()                             - oM2h.getE()) +
932                          672 * (oP1h.getE()                             - oM1h.getE())) / (840 * h);
933         jacobian[2][i] = (-3 * (oP4h.getI()                             - oM4h.getI()) +
934                           32 * (oP3h.getI()                             - oM3h.getI()) -
935                          168 * (oP2h.getI()                             - oM2h.getI()) +
936                          672 * (oP1h.getI()                             - oM1h.getI())) / (840 * h);
937         jacobian[3][i] = (-3 * (oP4h.getPerigeeArgument()               - oM4h.getPerigeeArgument()) +
938                           32 * (oP3h.getPerigeeArgument()               - oM3h.getPerigeeArgument()) -
939                          168 * (oP2h.getPerigeeArgument()               - oM2h.getPerigeeArgument()) +
940                          672 * (oP1h.getPerigeeArgument()               - oM1h.getPerigeeArgument())) / (840 * h);
941         jacobian[4][i] = (-3 * (oP4h.getRightAscensionOfAscendingNode() - oM4h.getRightAscensionOfAscendingNode()) +
942                           32 * (oP3h.getRightAscensionOfAscendingNode() - oM3h.getRightAscensionOfAscendingNode()) -
943                          168 * (oP2h.getRightAscensionOfAscendingNode() - oM2h.getRightAscensionOfAscendingNode()) +
944                          672 * (oP1h.getRightAscensionOfAscendingNode() - oM1h.getRightAscensionOfAscendingNode())) / (840 * h);
945         jacobian[5][i] = (-3 * (oP4h.getAnomaly(type)                   - oM4h.getAnomaly(type)) +
946                           32 * (oP3h.getAnomaly(type)                   - oM3h.getAnomaly(type)) -
947                          168 * (oP2h.getAnomaly(type)                   - oM2h.getAnomaly(type)) +
948                          672 * (oP1h.getAnomaly(type)                   - oM1h.getAnomaly(type))) / (840 * h);
949 
950     }
951 
952     @Test
953     public void testInterpolationWithDerivatives() {
954         doTestInterpolation(true,
955                             397, 4.01, 4.75e-4, 1.28e-7,
956                             2159, 1.05e7, 1.19e-3, 0.773);
957     }
958 
959     @Test
960     public void testInterpolationWithoutDerivatives() {
961         doTestInterpolation(false,
962                             397, 62.0, 4.75e-4, 2.87e-6,
963                             2159, 79365, 1.19e-3, 3.89e-3);
964     }
965 
966     private void doTestInterpolation(boolean useDerivatives,
967                                      double shiftPositionErrorWithin, double interpolationPositionErrorWithin,
968                                      double shiftEccentricityErrorWithin, double interpolationEccentricityErrorWithin,
969                                      double shiftPositionErrorSlightlyPast, double interpolationPositionErrorSlightlyPast,
970                                      double shiftEccentricityErrorSlightlyPast, double interpolationEccentricityErrorSlightlyPast)
971         {
972 
973         final double ehMu  = 3.9860047e14;
974         final double ae  = 6.378137e6;
975         final double c20 = -1.08263e-3;
976         final double c30 = 2.54e-6;
977         final double c40 = 1.62e-6;
978         final double c50 = 2.3e-7;
979         final double c60 = -5.5e-7;
980 
981         final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
982         final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
983         final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
984         final KeplerianOrbit initialOrbit = new KeplerianOrbit(new PVCoordinates(position, velocity),
985                                                                FramesFactory.getEME2000(), date, ehMu);
986 
987         EcksteinHechlerPropagator propagator =
988                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
989 
990         // set up a 5 points sample
991         List<Orbit> sample = new ArrayList<Orbit>();
992         for (double dt = 0; dt < 300.0; dt += 60.0) {
993             Orbit orbit = propagator.propagate(date.shiftedBy(dt)).getOrbit();
994             if (!useDerivatives) {
995                 // remove derivatives
996                 double[] stateVector = new double[6];
997                 orbit.getType().mapOrbitToArray(orbit, PositionAngle.TRUE, stateVector, null);
998                 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngle.TRUE,
999                                                         orbit.getDate(), orbit.getMu(), orbit.getFrame());
1000             }
1001             sample.add(orbit);
1002         }
1003 
1004         // well inside the sample, interpolation should be slightly better than Keplerian shift
1005         // the relative bad behaviour here is due to eccentricity, which cannot be
1006         // accurately interpolated with a polynomial in this case
1007         double maxShiftPositionError = 0;
1008         double maxInterpolationPositionError = 0;
1009         double maxShiftEccentricityError = 0;
1010         double maxInterpolationEccentricityError = 0;
1011         for (double dt = 0; dt < 241.0; dt += 1.0) {
1012             AbsoluteDate t         = initialOrbit.getDate().shiftedBy(dt);
1013             Vector3D shiftedP      = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
1014             Vector3D interpolatedP = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
1015             Vector3D propagatedP   = propagator.propagate(t).getPVCoordinates().getPosition();
1016             double shiftedE        = initialOrbit.shiftedBy(dt).getE();
1017             double interpolatedE   = initialOrbit.interpolate(t, sample).getE();
1018             double propagatedE     = propagator.propagate(t).getE();
1019             maxShiftPositionError = FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm());
1020             maxInterpolationPositionError = FastMath.max(maxInterpolationPositionError, interpolatedP.subtract(propagatedP).getNorm());
1021             maxShiftEccentricityError = FastMath.max(maxShiftEccentricityError, FastMath.abs(shiftedE - propagatedE));
1022             maxInterpolationEccentricityError = FastMath.max(maxInterpolationEccentricityError, FastMath.abs(interpolatedE - propagatedE));
1023         }
1024         Assert.assertEquals(shiftPositionErrorWithin,             maxShiftPositionError,             0.01 * shiftPositionErrorWithin);
1025         Assert.assertEquals(interpolationPositionErrorWithin,     maxInterpolationPositionError,     0.01 * interpolationPositionErrorWithin);
1026         Assert.assertEquals(shiftEccentricityErrorWithin,         maxShiftEccentricityError,         0.01 * shiftEccentricityErrorWithin);
1027         Assert.assertEquals(interpolationEccentricityErrorWithin, maxInterpolationEccentricityError, 0.01 * interpolationEccentricityErrorWithin);
1028 
1029         // slightly past sample end, bad eccentricity interpolation shows up
1030         // (in this case, interpolated eccentricity exceeds 1.0 btween 1900
1031         // and 1910s, while semi-majaxis remains positive, so this is not
1032         // even a proper hyperbolic orbit...)
1033         maxShiftPositionError = 0;
1034         maxInterpolationPositionError = 0;
1035         maxShiftEccentricityError = 0;
1036         maxInterpolationEccentricityError = 0;
1037         for (double dt = 240; dt < 600; dt += 1.0) {
1038             AbsoluteDate t         = initialOrbit.getDate().shiftedBy(dt);
1039             Vector3D shiftedP      = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
1040             Vector3D interpolatedP = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
1041             Vector3D propagatedP   = propagator.propagate(t).getPVCoordinates().getPosition();
1042             double shiftedE        = initialOrbit.shiftedBy(dt).getE();
1043             double interpolatedE   = initialOrbit.interpolate(t, sample).getE();
1044             double propagatedE     = propagator.propagate(t).getE();
1045             maxShiftPositionError = FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm());
1046             maxInterpolationPositionError = FastMath.max(maxInterpolationPositionError, interpolatedP.subtract(propagatedP).getNorm());
1047             maxShiftEccentricityError = FastMath.max(maxShiftEccentricityError, FastMath.abs(shiftedE - propagatedE));
1048             maxInterpolationEccentricityError = FastMath.max(maxInterpolationEccentricityError, FastMath.abs(interpolatedE - propagatedE));
1049         }
1050         Assert.assertEquals(shiftPositionErrorSlightlyPast,             maxShiftPositionError,             0.01 * shiftPositionErrorSlightlyPast);
1051         Assert.assertEquals(interpolationPositionErrorSlightlyPast,     maxInterpolationPositionError,     0.01 * interpolationPositionErrorSlightlyPast);
1052         Assert.assertEquals(shiftEccentricityErrorSlightlyPast,         maxShiftEccentricityError,         0.01 * shiftEccentricityErrorSlightlyPast);
1053         Assert.assertEquals(interpolationEccentricityErrorSlightlyPast, maxInterpolationEccentricityError, 0.01 * interpolationEccentricityErrorSlightlyPast);
1054 
1055     }
1056 
1057     @Test
1058     public void testPerfectlyEquatorialConversion() {
1059         KeplerianOrbit initial = new KeplerianOrbit(13378000.0, 0.05, 0.0, 0.0, FastMath.PI,
1060                                                     0.0, PositionAngle.MEAN,
1061                                                     FramesFactory.getEME2000(), date,
1062                                                     Constants.EIGEN5C_EARTH_MU);
1063         EquinoctialOrbit equ = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(initial);
1064         KeplerianOrbit converted = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(equ);
1065         Assert.assertEquals(FastMath.PI,
1066                             MathUtils.normalizeAngle(converted.getRightAscensionOfAscendingNode() +
1067                                                      converted.getPerigeeArgument(), FastMath.PI),
1068                             1.0e-10);
1069     }
1070 
1071     @Test
1072     public void testKeplerianDerivatives() {
1073 
1074         final KeplerianOrbit orbit = new KeplerianOrbit(new PVCoordinates(new Vector3D(-4947831., -3765382., -3708221.),
1075                                                                           new Vector3D(-2079., 5291., -7842.)),
1076                                                         FramesFactory.getEME2000(), date, 3.9860047e14);
1077         final Vector3D p = orbit.getPVCoordinates().getPosition();
1078         final Vector3D v = orbit.getPVCoordinates().getVelocity();
1079         final Vector3D a = orbit.getPVCoordinates().getAcceleration();
1080 
1081         // check that despite we did not provide acceleration, it got recomputed
1082         Assert.assertEquals(7.605422, a.getNorm(), 1.0e-6);
1083 
1084         // check velocity is the derivative of position
1085         Assert.assertEquals(differentiate(orbit, shifted -> shifted.getPVCoordinates().getPosition().getX()),
1086                             orbit.getPVCoordinates().getVelocity().getX(),
1087                             3.0e-12 * v.getNorm());
1088         Assert.assertEquals(differentiate(orbit, shifted -> shifted.getPVCoordinates().getPosition().getY()),
1089                             orbit.getPVCoordinates().getVelocity().getY(),
1090                             3.0e-12 * v.getNorm());
1091         Assert.assertEquals(differentiate(orbit, shifted -> shifted.getPVCoordinates().getPosition().getZ()),
1092                             orbit.getPVCoordinates().getVelocity().getZ(),
1093                             3.0e-12 * v.getNorm());
1094 
1095         // check acceleration is the derivative of velocity
1096         Assert.assertEquals(differentiate(orbit, shifted -> shifted.getPVCoordinates().getVelocity().getX()),
1097                             orbit.getPVCoordinates().getAcceleration().getX(),
1098                             3.0e-12 * a.getNorm());
1099         Assert.assertEquals(differentiate(orbit, shifted -> shifted.getPVCoordinates().getVelocity().getY()),
1100                             orbit.getPVCoordinates().getAcceleration().getY(),
1101                             3.0e-12 * a.getNorm());
1102         Assert.assertEquals(differentiate(orbit, shifted -> shifted.getPVCoordinates().getVelocity().getZ()),
1103                             orbit.getPVCoordinates().getAcceleration().getZ(),
1104                             3.0e-12 * a.getNorm());
1105 
1106         // check jerk is the derivative of acceleration
1107         final double r2 = p.getNormSq();
1108         final double r  = FastMath.sqrt(r2);
1109         Vector3D keplerianJerk = new Vector3D(-3 * Vector3D.dotProduct(p, v) / r2, a, -a.getNorm() / r, v);
1110         Assert.assertEquals(differentiate(orbit, shifted -> shifted.getPVCoordinates().getAcceleration().getX()),
1111                             keplerianJerk.getX(),
1112                             3.0e-12 * keplerianJerk.getNorm());
1113         Assert.assertEquals(differentiate(orbit, shifted -> shifted.getPVCoordinates().getAcceleration().getY()),
1114                             keplerianJerk.getY(),
1115                             3.0e-12 * keplerianJerk.getNorm());
1116         Assert.assertEquals(differentiate(orbit, shifted -> shifted.getPVCoordinates().getAcceleration().getZ()),
1117                             keplerianJerk.getZ(),
1118                             3.0e-12 * keplerianJerk.getNorm());
1119 
1120         Assert.assertTrue(Double.isNaN(orbit.getADot()));
1121         Assert.assertTrue(Double.isNaN(orbit.getEquinoctialExDot()));
1122         Assert.assertTrue(Double.isNaN(orbit.getEquinoctialEyDot()));
1123         Assert.assertTrue(Double.isNaN(orbit.getHxDot()));
1124         Assert.assertTrue(Double.isNaN(orbit.getHyDot()));
1125         Assert.assertTrue(Double.isNaN(orbit.getLvDot()));
1126         Assert.assertTrue(Double.isNaN(orbit.getLEDot()));
1127         Assert.assertTrue(Double.isNaN(orbit.getLMDot()));
1128         Assert.assertTrue(Double.isNaN(orbit.getEDot()));
1129         Assert.assertTrue(Double.isNaN(orbit.getIDot()));
1130         Assert.assertTrue(Double.isNaN(orbit.getPerigeeArgumentDot()));
1131         Assert.assertTrue(Double.isNaN(orbit.getRightAscensionOfAscendingNodeDot()));
1132         Assert.assertTrue(Double.isNaN(orbit.getTrueAnomalyDot()));
1133         Assert.assertTrue(Double.isNaN(orbit.getEccentricAnomalyDot()));
1134         Assert.assertTrue(Double.isNaN(orbit.getMeanAnomalyDot()));
1135         Assert.assertTrue(Double.isNaN(orbit.getAnomalyDot(PositionAngle.TRUE)));
1136         Assert.assertTrue(Double.isNaN(orbit.getAnomalyDot(PositionAngle.ECCENTRIC)));
1137         Assert.assertTrue(Double.isNaN(orbit.getAnomalyDot(PositionAngle.MEAN)));
1138 
1139     }
1140 
1141     private <S extends Function<KeplerianOrbit, Double>>
1142     double differentiate(KeplerianOrbit orbit, S picker) {
1143         final DSFactory factory = new DSFactory(1, 1);
1144         FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(8, 0.1);
1145         UnivariateDifferentiableFunction diff = differentiator.differentiate(new UnivariateFunction() {
1146             public double value(double dt) {
1147                 return picker.apply(orbit.shiftedBy(dt));
1148             }
1149         });
1150         return diff.value(factory.variable(0, 0.0)).getPartialDerivative(1);
1151      }
1152 
1153     @Test
1154     public void testNonKeplerianEllipticDerivatives() {
1155         final AbsoluteDate date         = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1156         final Vector3D     position     = new Vector3D(6896874.444705,  1956581.072644,  -147476.245054);
1157         final Vector3D     velocity     = new Vector3D(166.816407662, -1106.783301861, -7372.745712770);
1158         final Vector3D     acceleration = new Vector3D(-7.466182457944, -2.118153357345,  0.160004048437);
1159         final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity, acceleration);
1160         final Frame frame = FramesFactory.getEME2000();
1161         final double mu   = Constants.EIGEN5C_EARTH_MU;
1162         final KeplerianOrbit orbit = new KeplerianOrbit(pv, frame, mu);
1163 
1164         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getA()),
1165                             orbit.getADot(),
1166                             4.3e-8);
1167         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getEquinoctialEx()),
1168                             orbit.getEquinoctialExDot(),
1169                             2.1e-15);
1170         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getEquinoctialEy()),
1171                             orbit.getEquinoctialEyDot(),
1172                             5.3e-16);
1173         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getHx()),
1174                             orbit.getHxDot(),
1175                             1.6e-15);
1176         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getHy()),
1177                             orbit.getHyDot(),
1178                             7.3e-17);
1179         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLv()),
1180                             orbit.getLvDot(),
1181                             1.1e-14);
1182         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLE()),
1183                             orbit.getLEDot(),
1184                             7.2e-15);
1185         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLM()),
1186                             orbit.getLMDot(),
1187                             4.7e-15);
1188         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getE()),
1189                             orbit.getEDot(),
1190                             6.9e-16);
1191         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getI()),
1192                             orbit.getIDot(),
1193                             5.8e-16);
1194         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getPerigeeArgument()),
1195                             orbit.getPerigeeArgumentDot(),
1196                             1.5e-12);
1197         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getRightAscensionOfAscendingNode()),
1198                             orbit.getRightAscensionOfAscendingNodeDot(),
1199                             1.5e-15);
1200         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getTrueAnomaly()),
1201                             orbit.getTrueAnomalyDot(),
1202                             1.5e-12);
1203         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getEccentricAnomaly()),
1204                             orbit.getEccentricAnomalyDot(),
1205                             1.5e-12);
1206         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getMeanAnomaly()),
1207                             orbit.getMeanAnomalyDot(),
1208                             1.5e-12);
1209         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAnomaly(PositionAngle.TRUE)),
1210                             orbit.getAnomalyDot(PositionAngle.TRUE),
1211                             1.5e-12);
1212         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAnomaly(PositionAngle.ECCENTRIC)),
1213                             orbit.getAnomalyDot(PositionAngle.ECCENTRIC),
1214                             1.5e-12);
1215         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAnomaly(PositionAngle.MEAN)),
1216                             orbit.getAnomalyDot(PositionAngle.MEAN),
1217                             1.5e-12);
1218 
1219     }
1220 
1221     @Test
1222     public void testNonKeplerianHyperbolicDerivatives() {
1223         final AbsoluteDate date         = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1224         final Vector3D     position     = new Vector3D(224267911.905821, 290251613.109399, 45534292.777492);
1225         final Vector3D     velocity     = new Vector3D(-1494.068165293, 1124.771027677, 526.915286134);
1226         final Vector3D     acceleration = new Vector3D(-0.001295920501, -0.002233045187, -0.000349906292);
1227         final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity, acceleration);
1228         final Frame frame = FramesFactory.getEME2000();
1229         final double mu   = Constants.EIGEN5C_EARTH_MU;
1230         final KeplerianOrbit orbit = new KeplerianOrbit(pv, frame, mu);
1231 
1232         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getA()),
1233                             orbit.getADot(),
1234                             9.6e-8);
1235         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getEquinoctialEx()),
1236                             orbit.getEquinoctialExDot(),
1237                             2.8e-16);
1238         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getEquinoctialEy()),
1239                             orbit.getEquinoctialEyDot(),
1240                             3.6e-15);
1241         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getHx()),
1242                             orbit.getHxDot(),
1243                             1.4e-15);
1244         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getHy()),
1245                             orbit.getHyDot(),
1246                             9.4e-16);
1247         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLv()),
1248                             orbit.getLvDot(),
1249                             5.6e-16);
1250         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLE()),
1251                             orbit.getLEDot(),
1252                             9.0e-16);
1253         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLM()),
1254                             orbit.getLMDot(),
1255                             1.8e-15);
1256         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getE()),
1257                             orbit.getEDot(),
1258                             1.8e-15);
1259         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getI()),
1260                             orbit.getIDot(),
1261                             3.6e-15);
1262         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getPerigeeArgument()),
1263                             orbit.getPerigeeArgumentDot(),
1264                             9.4e-16);
1265         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getRightAscensionOfAscendingNode()),
1266                             orbit.getRightAscensionOfAscendingNodeDot(),
1267                             1.1e-15);
1268         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getTrueAnomaly()),
1269                             orbit.getTrueAnomalyDot(),
1270                             1.4e-15);
1271         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getEccentricAnomaly()),
1272                             orbit.getEccentricAnomalyDot(),
1273                             9.2e-16);
1274         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getMeanAnomaly()),
1275                             orbit.getMeanAnomalyDot(),
1276                             1.4e-15);
1277         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAnomaly(PositionAngle.TRUE)),
1278                             orbit.getAnomalyDot(PositionAngle.TRUE),
1279                             1.4e-15);
1280         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAnomaly(PositionAngle.ECCENTRIC)),
1281                             orbit.getAnomalyDot(PositionAngle.ECCENTRIC),
1282                             9.2e-16);
1283         Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAnomaly(PositionAngle.MEAN)),
1284                             orbit.getAnomalyDot(PositionAngle.MEAN),
1285                             1.4e-15);
1286 
1287     }
1288 
1289     private <S extends Function<KeplerianOrbit, Double>>
1290     double differentiate(TimeStampedPVCoordinates pv, Frame frame, double mu, S picker) {
1291         final DSFactory factory = new DSFactory(1, 1);
1292         FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(8, 0.1);
1293         UnivariateDifferentiableFunction diff = differentiator.differentiate(new UnivariateFunction() {
1294             public double value(double dt) {
1295                 return picker.apply(new KeplerianOrbit(pv.shiftedBy(dt), frame, mu));
1296             }
1297         });
1298         return diff.value(factory.variable(0, 0.0)).getPartialDerivative(1);
1299      }
1300 
1301     @Test
1302     public void testPositionAngleDerivatives() {
1303         final AbsoluteDate date         = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1304         final Vector3D     position     = new Vector3D(6896874.444705,  1956581.072644,  -147476.245054);
1305         final Vector3D     velocity     = new Vector3D(166.816407662, -1106.783301861, -7372.745712770);
1306         final Vector3D     acceleration = new Vector3D(-7.466182457944, -2.118153357345,  0.160004048437);
1307         final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity, acceleration);
1308         final Frame frame = FramesFactory.getEME2000();
1309         final double mu   = Constants.EIGEN5C_EARTH_MU;
1310         final KeplerianOrbit orbit = new KeplerianOrbit(pv, frame, mu);
1311 
1312         for (PositionAngle type : PositionAngle.values()) {
1313             final KeplerianOrbit rebuilt = new KeplerianOrbit(orbit.getA(),
1314                                                               orbit.getE(),
1315                                                               orbit.getI(),
1316                                                               orbit.getPerigeeArgument(),
1317                                                               orbit.getRightAscensionOfAscendingNode(),
1318                                                               orbit.getAnomaly(type),
1319                                                               orbit.getADot(),
1320                                                               orbit.getEDot(),
1321                                                               orbit.getIDot(),
1322                                                               orbit.getPerigeeArgumentDot(),
1323                                                               orbit.getRightAscensionOfAscendingNodeDot(),
1324                                                               orbit.getAnomalyDot(type),
1325                                                               type, orbit.getFrame(), orbit.getDate(), orbit.getMu());
1326             MatcherAssert.assertThat(rebuilt.getA(),                                relativelyCloseTo(orbit.getA(),                                1));
1327             MatcherAssert.assertThat(rebuilt.getE(),                                relativelyCloseTo(orbit.getE(),                                1));
1328             MatcherAssert.assertThat(rebuilt.getI(),                                relativelyCloseTo(orbit.getI(),                                1));
1329             MatcherAssert.assertThat(rebuilt.getPerigeeArgument(),                  relativelyCloseTo(orbit.getPerigeeArgument(),                  1));
1330             MatcherAssert.assertThat(rebuilt.getRightAscensionOfAscendingNode(),    relativelyCloseTo(orbit.getRightAscensionOfAscendingNode(),    1));
1331             MatcherAssert.assertThat(rebuilt.getADot(),                             relativelyCloseTo(orbit.getADot(),                             1));
1332             MatcherAssert.assertThat(rebuilt.getEDot(),                             relativelyCloseTo(orbit.getEDot(),                             1));
1333             MatcherAssert.assertThat(rebuilt.getIDot(),                             relativelyCloseTo(orbit.getIDot(),                             1));
1334             MatcherAssert.assertThat(rebuilt.getPerigeeArgumentDot(),               relativelyCloseTo(orbit.getPerigeeArgumentDot(),               1));
1335             MatcherAssert.assertThat(rebuilt.getRightAscensionOfAscendingNodeDot(), relativelyCloseTo(orbit.getRightAscensionOfAscendingNodeDot(), 1));
1336             for (PositionAngle type2 : PositionAngle.values()) {
1337                 MatcherAssert.assertThat(rebuilt.getAnomaly(type2),    relativelyCloseTo(orbit.getAnomaly(type2),    1));
1338                 MatcherAssert.assertThat(rebuilt.getAnomalyDot(type2), relativelyCloseTo(orbit.getAnomalyDot(type2), 1));
1339             }
1340         }
1341 
1342     }
1343 
1344     @Test
1345     public void testPositionAngleHyperbolicDerivatives() {
1346         final AbsoluteDate date         = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1347         final Vector3D     position     = new Vector3D(224267911.905821, 290251613.109399, 45534292.777492);
1348         final Vector3D     velocity     = new Vector3D(-1494.068165293, 1124.771027677, 526.915286134);
1349         final Vector3D     acceleration = new Vector3D(-0.001295920501, -0.002233045187, -0.000349906292);
1350         final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity, acceleration);
1351         final Frame frame = FramesFactory.getEME2000();
1352         final double mu   = Constants.EIGEN5C_EARTH_MU;
1353         final KeplerianOrbit orbit = new KeplerianOrbit(pv, frame, mu);
1354 
1355         for (PositionAngle type : PositionAngle.values()) {
1356             final KeplerianOrbit rebuilt = new KeplerianOrbit(orbit.getA(),
1357                                                               orbit.getE(),
1358                                                               orbit.getI(),
1359                                                               orbit.getPerigeeArgument(),
1360                                                               orbit.getRightAscensionOfAscendingNode(),
1361                                                               orbit.getAnomaly(type),
1362                                                               orbit.getADot(),
1363                                                               orbit.getEDot(),
1364                                                               orbit.getIDot(),
1365                                                               orbit.getPerigeeArgumentDot(),
1366                                                               orbit.getRightAscensionOfAscendingNodeDot(),
1367                                                               orbit.getAnomalyDot(type),
1368                                                               type, orbit.getFrame(), orbit.getDate(), orbit.getMu());
1369             MatcherAssert.assertThat(rebuilt.getA(),                                relativelyCloseTo(orbit.getA(),                                1));
1370             MatcherAssert.assertThat(rebuilt.getE(),                                relativelyCloseTo(orbit.getE(),                                1));
1371             MatcherAssert.assertThat(rebuilt.getI(),                                relativelyCloseTo(orbit.getI(),                                1));
1372             MatcherAssert.assertThat(rebuilt.getPerigeeArgument(),                  relativelyCloseTo(orbit.getPerigeeArgument(),                  1));
1373             MatcherAssert.assertThat(rebuilt.getRightAscensionOfAscendingNode(),    relativelyCloseTo(orbit.getRightAscensionOfAscendingNode(),    1));
1374             MatcherAssert.assertThat(rebuilt.getADot(),                             relativelyCloseTo(orbit.getADot(),                             1));
1375             MatcherAssert.assertThat(rebuilt.getEDot(),                             relativelyCloseTo(orbit.getEDot(),                             1));
1376             MatcherAssert.assertThat(rebuilt.getIDot(),                             relativelyCloseTo(orbit.getIDot(),                             1));
1377             MatcherAssert.assertThat(rebuilt.getPerigeeArgumentDot(),               relativelyCloseTo(orbit.getPerigeeArgumentDot(),               1));
1378             MatcherAssert.assertThat(rebuilt.getRightAscensionOfAscendingNodeDot(), relativelyCloseTo(orbit.getRightAscensionOfAscendingNodeDot(), 1));
1379             for (PositionAngle type2 : PositionAngle.values()) {
1380                 MatcherAssert.assertThat(rebuilt.getAnomaly(type2),    relativelyCloseTo(orbit.getAnomaly(type2),    2));
1381                 MatcherAssert.assertThat(rebuilt.getAnomalyDot(type2), relativelyCloseTo(orbit.getAnomalyDot(type2), 4));
1382             }
1383         }
1384 
1385     }
1386 
1387     @Test
1388     public void testSerialization()
1389             throws IOException, ClassNotFoundException {
1390         Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
1391         Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
1392         PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
1393         KeplerianOrbit orbit = new KeplerianOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
1394         Assert.assertEquals(42255170.003, orbit.getA(), 1.0e-3);
1395 
1396         ByteArrayOutputStream bos = new ByteArrayOutputStream();
1397         ObjectOutputStream    oos = new ObjectOutputStream(bos);
1398         oos.writeObject(orbit);
1399 
1400         Assert.assertTrue(bos.size() > 280);
1401         Assert.assertTrue(bos.size() < 330);
1402 
1403         ByteArrayInputStream  bis = new ByteArrayInputStream(bos.toByteArray());
1404         ObjectInputStream     ois = new ObjectInputStream(bis);
1405         KeplerianOrbit deserialized  = (KeplerianOrbit) ois.readObject();
1406         Assert.assertEquals(orbit.getA(), deserialized.getA(), 1.0e-10);
1407         Assert.assertEquals(orbit.getE(), deserialized.getE(), 1.0e-10);
1408         Assert.assertEquals(orbit.getI(), deserialized.getI(), 1.0e-10);
1409         Assert.assertEquals(orbit.getPerigeeArgument(), deserialized.getPerigeeArgument(), 1.0e-10);
1410         Assert.assertEquals(orbit.getRightAscensionOfAscendingNode(), deserialized.getRightAscensionOfAscendingNode(), 1.0e-10);
1411         Assert.assertEquals(orbit.getTrueAnomaly(), deserialized.getTrueAnomaly(), 1.0e-10);
1412         Assert.assertTrue(Double.isNaN(orbit.getADot()) && Double.isNaN(deserialized.getADot()));
1413         Assert.assertTrue(Double.isNaN(orbit.getEDot()) && Double.isNaN(deserialized.getEDot()));
1414         Assert.assertTrue(Double.isNaN(orbit.getIDot()) && Double.isNaN(deserialized.getIDot()));
1415         Assert.assertTrue(Double.isNaN(orbit.getPerigeeArgumentDot()) && Double.isNaN(deserialized.getPerigeeArgumentDot()));
1416         Assert.assertTrue(Double.isNaN(orbit.getRightAscensionOfAscendingNodeDot()) && Double.isNaN(deserialized.getRightAscensionOfAscendingNodeDot()));
1417         Assert.assertTrue(Double.isNaN(orbit.getTrueAnomalyDot()) && Double.isNaN(deserialized.getTrueAnomalyDot()));
1418         Assert.assertEquals(orbit.getDate(), deserialized.getDate());
1419         Assert.assertEquals(orbit.getMu(), deserialized.getMu(), 1.0e-10);
1420         Assert.assertEquals(orbit.getFrame().getName(), deserialized.getFrame().getName());
1421 
1422     }
1423 
1424     @Test
1425     public void testSerializationWithDerivatives()
1426             throws IOException, ClassNotFoundException {
1427         Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
1428         Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
1429         double r2 = position.getNormSq();
1430         double r  = FastMath.sqrt(r2);
1431         Vector3D acceleration = new Vector3D(-mu / (r * r2), position,
1432                                              1, new Vector3D(-0.1, 0.2, 0.3));
1433         PVCoordinates pvCoordinates = new PVCoordinates(position, velocity, acceleration);
1434         KeplerianOrbit orbit = new KeplerianOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
1435         Assert.assertEquals(42255170.003, orbit.getA(), 1.0e-3);
1436 
1437         ByteArrayOutputStream bos = new ByteArrayOutputStream();
1438         ObjectOutputStream    oos = new ObjectOutputStream(bos);
1439         oos.writeObject(orbit);
1440 
1441         Assert.assertTrue(bos.size() > 330);
1442         Assert.assertTrue(bos.size() < 380);
1443 
1444         ByteArrayInputStream  bis = new ByteArrayInputStream(bos.toByteArray());
1445         ObjectInputStream     ois = new ObjectInputStream(bis);
1446         KeplerianOrbit deserialized  = (KeplerianOrbit) ois.readObject();
1447         Assert.assertEquals(orbit.getA(), deserialized.getA(), 1.0e-10);
1448         Assert.assertEquals(orbit.getE(), deserialized.getE(), 1.0e-10);
1449         Assert.assertEquals(orbit.getI(), deserialized.getI(), 1.0e-10);
1450         Assert.assertEquals(orbit.getPerigeeArgument(), deserialized.getPerigeeArgument(), 1.0e-10);
1451         Assert.assertEquals(orbit.getRightAscensionOfAscendingNode(), deserialized.getRightAscensionOfAscendingNode(), 1.0e-10);
1452         Assert.assertEquals(orbit.getTrueAnomaly(), deserialized.getTrueAnomaly(), 1.0e-10);
1453         Assert.assertEquals(orbit.getADot(), deserialized.getADot(), 1.0e-10);
1454         Assert.assertEquals(orbit.getEDot(), deserialized.getEDot(), 1.0e-10);
1455         Assert.assertEquals(orbit.getIDot(), deserialized.getIDot(), 1.0e-10);
1456         Assert.assertEquals(orbit.getPerigeeArgumentDot(), deserialized.getPerigeeArgumentDot(), 1.0e-10);
1457         Assert.assertEquals(orbit.getRightAscensionOfAscendingNodeDot(), deserialized.getRightAscensionOfAscendingNodeDot(), 1.0e-10);
1458         Assert.assertEquals(orbit.getTrueAnomalyDot(), deserialized.getTrueAnomalyDot(), 1.0e-10);
1459         Assert.assertEquals(orbit.getDate(), deserialized.getDate());
1460         Assert.assertEquals(orbit.getMu(), deserialized.getMu(), 1.0e-10);
1461         Assert.assertEquals(orbit.getFrame().getName(), deserialized.getFrame().getName());
1462 
1463     }
1464 
1465     @Test
1466     public void testEquatorialRetrograde() {
1467         Vector3D position = new Vector3D(10000000.0, 0.0, 0.0);
1468         Vector3D velocity = new Vector3D(0.0, -6500.0, 1.0e-10);
1469         double r2 = position.getNormSq();
1470         double r  = FastMath.sqrt(r2);
1471         Vector3D acceleration = new Vector3D(-mu / (r * r2), position,
1472                                              1, new Vector3D(-0.1, 0.2, 0.3));
1473         PVCoordinates pvCoordinates = new PVCoordinates(position, velocity, acceleration);
1474         KeplerianOrbit orbit = new KeplerianOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
1475         Assert.assertEquals(10637829.465, orbit.getA(), 1.0e-3);
1476         Assert.assertEquals(-738.145, orbit.getADot(), 1.0e-3);
1477         Assert.assertEquals(0.05995861, orbit.getE(), 1.0e-8);
1478         Assert.assertEquals(-6.523e-5, orbit.getEDot(), 1.0e-8);
1479         Assert.assertEquals(FastMath.PI, orbit.getI(), 2.0e-14);
1480         Assert.assertEquals(-4.615e-5, orbit.getIDot(), 1.0e-8);
1481         Assert.assertTrue(Double.isNaN(orbit.getHx()));
1482         Assert.assertTrue(Double.isNaN(orbit.getHxDot()));
1483         Assert.assertTrue(Double.isNaN(orbit.getHy()));
1484         Assert.assertTrue(Double.isNaN(orbit.getHyDot()));
1485     }
1486 
1487     @Test
1488     public void testDerivativesConversionSymmetry() {
1489         final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:01:20.000", TimeScalesFactory.getUTC());
1490         Vector3D position     = new Vector3D(6893443.400234382, 1886406.1073757345, -589265.1150359757);
1491         Vector3D velocity     = new Vector3D(-281.1261461082365, -1231.6165642450928, -7348.756363469432);
1492         Vector3D acceleration = new Vector3D(-7.460341170581685, -2.0415957334584527, 0.6393322823627762);
1493         PVCoordinates pvCoordinates = new PVCoordinates( position, velocity, acceleration);
1494         KeplerianOrbit orbit = new KeplerianOrbit(pvCoordinates, FramesFactory.getEME2000(),
1495                                                   date, Constants.EIGEN5C_EARTH_MU);
1496         Assert.assertTrue(orbit.hasDerivatives());
1497         double r2 = position.getNormSq();
1498         double r  = FastMath.sqrt(r2);
1499         Vector3D keplerianAcceleration = new Vector3D(-orbit.getMu() / (r2 * r), position);
1500         Assert.assertEquals(0.0101, Vector3D.distance(keplerianAcceleration, acceleration), 1.0e-4);
1501 
1502         for (OrbitType type : OrbitType.values()) {
1503             Orbit converted = type.convertType(orbit);
1504             Assert.assertTrue(converted.hasDerivatives());
1505             KeplerianOrbit rebuilt = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(converted);
1506             Assert.assertTrue(rebuilt.hasDerivatives());
1507             Assert.assertEquals(orbit.getADot(),                             rebuilt.getADot(),                             3.0e-13);
1508             Assert.assertEquals(orbit.getEDot(),                             rebuilt.getEDot(),                             1.0e-15);
1509             Assert.assertEquals(orbit.getIDot(),                             rebuilt.getIDot(),                             1.0e-15);
1510             Assert.assertEquals(orbit.getPerigeeArgumentDot(),               rebuilt.getPerigeeArgumentDot(),               1.0e-15);
1511             Assert.assertEquals(orbit.getRightAscensionOfAscendingNodeDot(), rebuilt.getRightAscensionOfAscendingNodeDot(), 1.0e-15);
1512             Assert.assertEquals(orbit.getTrueAnomalyDot(),                   rebuilt.getTrueAnomalyDot(),                   1.0e-15);
1513         }
1514 
1515     }
1516 
1517     @Test
1518     public void testDerivativesConversionSymmetryHyperbolic() {
1519         final AbsoluteDate date         = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1520         final Vector3D     position     = new Vector3D(224267911.905821, 290251613.109399, 45534292.777492);
1521         final Vector3D     velocity     = new Vector3D(-1494.068165293, 1124.771027677, 526.915286134);
1522         final Vector3D     acceleration = new Vector3D(-0.001295920501, -0.002233045187, -0.000349906292);
1523         PVCoordinates pvCoordinates = new PVCoordinates( position, velocity, acceleration);
1524         KeplerianOrbit orbit = new KeplerianOrbit(pvCoordinates, FramesFactory.getEME2000(),
1525                                                   date, Constants.EIGEN5C_EARTH_MU);
1526         Assert.assertTrue(orbit.hasDerivatives());
1527         double r2 = position.getNormSq();
1528         double r  = FastMath.sqrt(r2);
1529         Vector3D keplerianAcceleration = new Vector3D(-orbit.getMu() / (r2 * r), position);
1530         Assert.assertEquals(4.78e-4, Vector3D.distance(keplerianAcceleration, acceleration), 1.0e-6);
1531 
1532         OrbitType type = OrbitType.CARTESIAN;
1533         Orbit converted = type.convertType(orbit);
1534         Assert.assertTrue(converted.hasDerivatives());
1535         KeplerianOrbit rebuilt = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(converted);
1536         Assert.assertTrue(rebuilt.hasDerivatives());
1537         Assert.assertEquals(orbit.getADot(),                             rebuilt.getADot(),                             3.0e-13);
1538         Assert.assertEquals(orbit.getEDot(),                             rebuilt.getEDot(),                             1.0e-15);
1539         Assert.assertEquals(orbit.getIDot(),                             rebuilt.getIDot(),                             1.0e-15);
1540         Assert.assertEquals(orbit.getPerigeeArgumentDot(),               rebuilt.getPerigeeArgumentDot(),               1.0e-15);
1541         Assert.assertEquals(orbit.getRightAscensionOfAscendingNodeDot(), rebuilt.getRightAscensionOfAscendingNodeDot(), 1.0e-15);
1542         Assert.assertEquals(orbit.getTrueAnomalyDot(),                   rebuilt.getTrueAnomalyDot(),                   1.0e-15);
1543 
1544     }
1545 
1546     @Test
1547     public void testToString() {
1548         Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
1549         Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
1550         PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
1551         KeplerianOrbit orbit = new KeplerianOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
1552         Assert.assertEquals("Keplerian parameters: {a: 4.225517000282565E7; e: 0.002146216321416967; i: 0.20189257051515358; pa: 13.949966363606599; raan: -87.91788415673473; v: -151.79096272977213;}",
1553                             orbit.toString());
1554     }
1555 
1556     @Test
1557     public void testCopyNonKeplerianAcceleration() {
1558 
1559         final Frame eme2000     = FramesFactory.getEME2000();
1560 
1561         // Define GEO satellite position
1562         final Vector3D position = new Vector3D(42164140, 0, 0);
1563         // Build PVCoodrinates starting from its position and computing the corresponding circular velocity
1564         final PVCoordinates pv  = new PVCoordinates(position,
1565                                        new Vector3D(0, FastMath.sqrt(mu / position.getNorm()), 0));
1566         // Build a KeplerianOrbit in eme2000
1567         final Orbit orbit = new KeplerianOrbit(pv, eme2000, date, mu);
1568 
1569         // Build another KeplerianOrbit as a copy of the first one
1570         final Orbit orbitCopy = new KeplerianOrbit(orbit);
1571 
1572         // Shift the orbit of a time-interval
1573         final Orbit shiftedOrbit = orbit.shiftedBy(10); // This works good
1574         final Orbit shiftedOrbitCopy = orbitCopy.shiftedBy(10); // This does not work
1575 
1576         Assert.assertEquals(0.0,
1577                             Vector3D.distance(shiftedOrbit.getPVCoordinates().getPosition(),
1578                                               shiftedOrbitCopy.getPVCoordinates().getPosition()),
1579                             1.0e-10);
1580         Assert.assertEquals(0.0,
1581                             Vector3D.distance(shiftedOrbit.getPVCoordinates().getVelocity(),
1582                                               shiftedOrbitCopy.getPVCoordinates().getVelocity()),
1583                             1.0e-10);
1584 
1585     }
1586 
1587     @Test
1588     public void testIssue544() {
1589         // Initial parameters
1590         // In order to test the issue, we volontary set the anomaly at Double.NaN.
1591         double e = 0.7311;
1592         double anomaly = Double.NaN;
1593         // Computes the elliptic eccentric anomaly 
1594         double E = KeplerianOrbit.meanToEllipticEccentric(anomaly, e);
1595         // Verify that an infinite loop did not occur
1596         Assert.assertTrue(Double.isNaN(E));  
1597     }
1598 
1599     @Test
1600     public void testIssue674() {
1601         try {
1602             new KeplerianOrbit(24464560.0, -0.7311, 0.122138, 3.10686, 1.00681,
1603                                0.048363, PositionAngle.MEAN,
1604                                FramesFactory.getEME2000(), date, mu);
1605             Assert.fail("an exception should have been thrown");
1606         } catch (OrekitException oe) {
1607             Assert.assertEquals(OrekitMessages.INVALID_PARAMETER_RANGE, oe.getSpecifier());
1608             Assert.assertEquals("eccentricity", oe.getParts()[0]);
1609             Assert.assertEquals(-0.7311, ((Double) oe.getParts()[1]).doubleValue(), Double.MIN_VALUE);
1610             Assert.assertEquals(0.0, ((Double) oe.getParts()[2]).doubleValue(), Double.MIN_VALUE);
1611             Assert.assertEquals(Double.POSITIVE_INFINITY, ((Double) oe.getParts()[3]).doubleValue(), Double.MIN_VALUE);
1612         }
1613     }
1614 
1615     @Test
1616     public void testNormalize() {
1617         KeplerianOrbit withoutDerivatives =
1618                         new KeplerianOrbit(42166712.0, 0.005, 1.6,
1619                                            -0.3, 1.25, 0.4, PositionAngle.MEAN,
1620                                            FramesFactory.getEME2000(), date, mu);
1621         KeplerianOrbit ref =
1622                         new KeplerianOrbit(24000000.0, 0.012, 1.9,
1623                                           -6.28, 6.28, 12.56, PositionAngle.MEAN,
1624                                           FramesFactory.getEME2000(), date, mu);
1625 
1626         KeplerianOrbit normalized1 = (KeplerianOrbit) OrbitType.KEPLERIAN.normalize(withoutDerivatives, ref);
1627         Assert.assertFalse(normalized1.hasDerivatives());
1628         Assert.assertEquals(0.0, normalized1.getA() - withoutDerivatives.getA(), 1.0e-6);
1629         Assert.assertEquals(0.0, normalized1.getE() - withoutDerivatives.getE(), 1.0e-10);
1630         Assert.assertEquals(0.0, normalized1.getI() - withoutDerivatives.getI(), 1.0e-10);
1631         Assert.assertEquals(-MathUtils.TWO_PI, normalized1.getPerigeeArgument()               - withoutDerivatives.getPerigeeArgument(),               1.0e-10);
1632         Assert.assertEquals(+MathUtils.TWO_PI, normalized1.getRightAscensionOfAscendingNode() - withoutDerivatives.getRightAscensionOfAscendingNode(), 1.0e-10);
1633         Assert.assertEquals(2 * MathUtils.TWO_PI, normalized1.getTrueAnomaly()                - withoutDerivatives.getTrueAnomaly(),                   1.0e-10);
1634         Assert.assertTrue(Double.isNaN(normalized1.getADot()));
1635         Assert.assertTrue(Double.isNaN(normalized1.getEDot()));
1636         Assert.assertTrue(Double.isNaN(normalized1.getIDot()));
1637         Assert.assertTrue(Double.isNaN(normalized1.getPerigeeArgumentDot()));
1638         Assert.assertTrue(Double.isNaN(normalized1.getRightAscensionOfAscendingNodeDot()));
1639         Assert.assertTrue(Double.isNaN(normalized1.getTrueAnomalyDot()));
1640 
1641         double[] p = new double[6];
1642         OrbitType.KEPLERIAN.mapOrbitToArray(withoutDerivatives, PositionAngle.TRUE, p, null);
1643         KeplerianOrbit withDerivatives = (KeplerianOrbit) OrbitType.KEPLERIAN.mapArrayToOrbit(p,
1644                                                                                               new double[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 },
1645                                                                                               PositionAngle.TRUE,
1646                                                                                               withoutDerivatives.getDate(),
1647                                                                                               withoutDerivatives.getMu(),
1648                                                                                               withoutDerivatives.getFrame());
1649         KeplerianOrbit normalized2 = (KeplerianOrbit) OrbitType.KEPLERIAN.normalize(withDerivatives, ref);
1650         Assert.assertTrue(normalized2.hasDerivatives());
1651         Assert.assertEquals(0.0, normalized2.getA() - withoutDerivatives.getA(), 1.0e-6);
1652         Assert.assertEquals(0.0, normalized2.getE() - withoutDerivatives.getE(), 1.0e-10);
1653         Assert.assertEquals(0.0, normalized2.getI() - withoutDerivatives.getI(), 1.0e-10);
1654         Assert.assertEquals(-MathUtils.TWO_PI, normalized2.getPerigeeArgument()               - withoutDerivatives.getPerigeeArgument(),               1.0e-10);
1655         Assert.assertEquals(+MathUtils.TWO_PI, normalized2.getRightAscensionOfAscendingNode() - withoutDerivatives.getRightAscensionOfAscendingNode(), 1.0e-10);
1656         Assert.assertEquals(2 * MathUtils.TWO_PI, normalized2.getTrueAnomaly()                - withoutDerivatives.getTrueAnomaly(),                   1.0e-10);
1657         Assert.assertEquals(0.0, normalized2.getADot() - withDerivatives.getADot(), 1.0e-10);
1658         Assert.assertEquals(0.0, normalized2.getEDot() - withDerivatives.getEDot(), 1.0e-10);
1659         Assert.assertEquals(0.0, normalized2.getIDot() - withDerivatives.getIDot(), 1.0e-10);
1660         Assert.assertEquals(0.0, normalized2.getPerigeeArgumentDot()               - withDerivatives.getPerigeeArgumentDot(),               1.0e-10);
1661         Assert.assertEquals(0.0, normalized2.getRightAscensionOfAscendingNodeDot() - withDerivatives.getRightAscensionOfAscendingNodeDot(), 1.0e-10);
1662         Assert.assertEquals(0.0, normalized2.getTrueAnomalyDot()                   - withDerivatives.getTrueAnomalyDot(),                   1.0e-10);
1663 
1664     }
1665     
1666     @Test
1667     public void testKeplerianToPvToKeplerian() {
1668         // setup
1669         Frame eci = FramesFactory.getGCRF();
1670         AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
1671         double mu = Constants.EGM96_EARTH_MU;
1672         KeplerianOrbit orbit = new KeplerianOrbit(6878e3, 0, 0, 0, 0, 0,
1673                 PositionAngle.TRUE, eci, date, mu);
1674 
1675         // action
1676         KeplerianOrbit pvOrbit = new KeplerianOrbit(orbit.getPVCoordinates(), eci, mu);
1677 
1678         // verify
1679         Assert.assertEquals(orbit.shiftedBy(1).getA(), pvOrbit.shiftedBy(1).getA(), 1e-6);
1680         Assert.assertEquals(orbit.shiftedBy(1).getE(), pvOrbit.shiftedBy(1).getE(), 1e-14);
1681         Assert.assertEquals(orbit.shiftedBy(1).getI(), pvOrbit.shiftedBy(1).getI(), 1e-10);
1682         Assert.assertEquals(orbit.shiftedBy(1).getPerigeeArgument(), pvOrbit.shiftedBy(1).getPerigeeArgument(), 1e-10);
1683         Assert.assertEquals(orbit.shiftedBy(1).getRightAscensionOfAscendingNode(), pvOrbit.shiftedBy(1).getRightAscensionOfAscendingNode(), 1e-10);
1684         Assert.assertEquals(orbit.shiftedBy(1).getTrueAnomaly(), pvOrbit.shiftedBy(1).getTrueAnomaly(), 1e-10);
1685 
1686     }
1687 
1688 
1689     @Before
1690     public void setUp() {
1691 
1692         Utils.setDataRoot("regular-data");
1693 
1694         // Computation date
1695         date = AbsoluteDate.J2000_EPOCH;
1696 
1697         // Body mu
1698         mu = 3.9860047e14;
1699     }
1700 
1701     @After
1702     public void tearDown() {
1703         date = null;
1704     }
1705 
1706 }