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