1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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.MatcherAssert;
31 import org.hipparchus.analysis.UnivariateFunction;
32 import org.hipparchus.analysis.differentiation.DSFactory;
33 import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
34 import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
35 import org.hipparchus.geometry.euclidean.threed.Vector3D;
36 import org.hipparchus.linear.MatrixUtils;
37 import org.hipparchus.linear.RealMatrixPreservingVisitor;
38 import org.hipparchus.util.FastMath;
39 import org.hipparchus.util.MathUtils;
40 import org.junit.After;
41 import org.junit.Assert;
42 import org.junit.Before;
43 import org.junit.Test;
44 import org.orekit.Utils;
45 import org.orekit.errors.OrekitIllegalArgumentException;
46 import org.orekit.errors.OrekitMessages;
47 import org.orekit.frames.Frame;
48 import org.orekit.frames.FramesFactory;
49 import org.orekit.frames.Transform;
50 import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
51 import org.orekit.time.AbsoluteDate;
52 import org.orekit.time.TimeScalesFactory;
53 import org.orekit.utils.Constants;
54 import org.orekit.utils.PVCoordinates;
55 import org.orekit.utils.TimeStampedPVCoordinates;
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 Assert.assertEquals(param.getA(), circ.getA(), Utils.epsilonTest * circ.getA());
86 Assert.assertEquals(param.getEquinoctialEx(), circ.getEquinoctialEx(), Utils.epsilonE * FastMath.abs(circ.getE()));
87 Assert.assertEquals(param.getEquinoctialEy(), circ.getEquinoctialEy(), Utils.epsilonE * FastMath.abs(circ.getE()));
88 Assert.assertEquals(param.getHx(), circ.getHx(), Utils.epsilonAngle * FastMath.abs(circ.getI()));
89 Assert.assertEquals(param.getHy(), circ.getHy(), Utils.epsilonAngle * FastMath.abs(circ.getI()));
90 Assert.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 Assert.assertEquals(paramCir.getA(), circCir.getA(), Utils.epsilonTest * circCir.getA());
114 Assert.assertEquals(paramCir.getEquinoctialEx(), circCir.getEquinoctialEx(), Utils.epsilonEcir * FastMath.abs(circCir.getE()));
115 Assert.assertEquals(paramCir.getEquinoctialEy(), circCir.getEquinoctialEy(), Utils.epsilonEcir * FastMath.abs(circCir.getE()));
116 Assert.assertEquals(paramCir.getHx(), circCir.getHx(), Utils.epsilonAngle * FastMath.abs(circCir.getI()));
117 Assert.assertEquals(paramCir.getHy(), circCir.getHy(), Utils.epsilonAngle * FastMath.abs(circCir.getI()));
118 Assert.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 Assert.assertEquals(2 / r - v * v / mu, 1 / circ.getA(), 1.0e-7);
147
148 Assert.assertEquals( 0.233745668678733e+08, pos.getX(), Utils.epsilonTest * r);
149 Assert.assertEquals(-0.350998914352669e+08, pos.getY(), Utils.epsilonTest * r);
150 Assert.assertEquals(-0.150053723123334e+04, pos.getZ(), Utils.epsilonTest * r);
151
152 Assert.assertEquals(2558.7096558809967, vel.getX(), Utils.epsilonTest * v);
153 Assert.assertEquals(1704.1586039092576, vel.getY(), Utils.epsilonTest * v);
154 Assert.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 Assert.assertEquals(42166712.000, circ.getA(), Utils.epsilonTest * kep.getA());
179 Assert.assertEquals(0.110283316961361e-03, kep.getE(), Utils.epsilonE * FastMath.abs(kep.getE()));
180 Assert.assertEquals(0.166901168553917e-03, kep.getI(),
181 Utils.epsilonAngle * FastMath.abs(kep.getI()));
182 Assert.assertEquals(MathUtils.normalizeAngle(-3.87224326008837, kep.getPerigeeArgument()),
183 kep.getPerigeeArgument(),
184 Utils.epsilonTest * FastMath.abs(kep.getPerigeeArgument()));
185 Assert.assertEquals(MathUtils.normalizeAngle(5.51473467358854, kep.getRightAscensionOfAscendingNode()),
186 kep.getRightAscensionOfAscendingNode(),
187 Utils.epsilonTest * FastMath.abs(kep.getRightAscensionOfAscendingNode()));
188 Assert.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 Assert.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 Assert.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 Assert.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
241 Assert.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
242 Assert.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 Assert.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
252 Assert.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
253 Assert.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 Assert.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
262 Assert.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
263 Assert.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 Assert.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
288 Assert.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
289 Assert.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 Assert.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
299 Assert.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
300 Assert.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 Assert.assertEquals(p.getAlphaV() + raan, lv, Utils.epsilonAngle * FastMath.abs(lv));
309 Assert.assertEquals(p.getAlphaE() + raan, lE, Utils.epsilonAngle * FastMath.abs(lE));
310 Assert.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 Assert.assertEquals(a * epsilon * epsilon / ksi,
338 p.getPVCoordinates().getPosition().getNorm(),
339 Utils.epsilonTest * FastMath.abs(p.getPVCoordinates().getPosition().getNorm()));
340 Assert.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 Assert.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 Assert.assertEquals(0.0, orbit.getI(), 2.0e-14);
368 Assert.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 Assert.assertEquals(a * epsilon * epsilon / ksi,
395 pCirEqua.getPVCoordinates().getPosition().getNorm(),
396 Utils.epsilonTest * FastMath.abs(pCirEqua.getPVCoordinates().getPosition().getNorm()));
397 Assert.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 Assert.assertTrue((position.getNorm() - apogeeRadius) <= ( apogeeRadius * Utils.epsilonTest));
430 Assert.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 Assert.assertTrue(FastMath.abs(Vector3D.dotProduct(position, momentum)) < Utils.epsilonTest);
440
441 Assert.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 Assert.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 Assert.assertTrue((position.getNorm() - apogeeRadius) <= ( apogeeRadius * Utils.epsilonTest));
476 Assert.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 Assert.assertTrue(FastMath.abs(Vector3D.dotProduct(position, momentum)) < Utils.epsilonTest);
486
487 Assert.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 Assert.assertEquals(0.0, positionOffset.getNorm(), position.getNorm() * Utils.epsilonTest);
508 Assert.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 Assert.assertEquals(0.0, positionOffset.getNorm(), position.getNorm() * Utils.epsilonTest);
525 Assert.assertEquals(0.0, velocityOffset.getNorm(), velocity.getNorm() * Utils.epsilonTest);
526
527 }
528
529 @Test(expected=IllegalArgumentException.class)
530 public void testNonInertialFrame() throws IllegalArgumentException {
531
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 @Test
541 public void testJacobianReference() {
542
543 AbsoluteDate dateTca = new AbsoluteDate(2000, 04, 01, 0, 0, 0.000, TimeScalesFactory.getUTC());
544 double mu = 3.986004415e+14;
545 CircularOrbit orbCir = new CircularOrbit(7000000.0, 0.01, -0.02, 1.2, 2.1,
546 0.7, PositionAngle.MEAN,
547 FramesFactory.getEME2000(), dateTca, mu);
548
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 Vector3D pRef = new Vector3D(-4106905.105389204807580, 3603162.539798960555345, 4439730.167038885876536);
594 Vector3D vRef = new Vector3D(740.132407342422994, -5308.773280141396754, 5250.338353483879473);
595 double[][] jRef = {
596 { -1.1535467596325562, 1.0120556393573172, 1.2470306024626943, 181.96913090864561, -1305.2162699469984, 1290.8494448855752 },
597 { -5.07367368325471104E-008, -1.27870567070456834E-008, 1.31544531338558113E-007, -3.09332106417043592E-005, -9.60781276304445404E-005, 1.91506964883791605E-004 },
598 { -6.59428471712402018E-008, 1.24561703203882533E-007, -1.41907027322388158E-008, 7.63442601186485441E-005, -1.77446722746170009E-004, 5.99464401287846734E-005 },
599 { 7.55079920652274275E-008, 4.41606835295069131E-008, 3.40079310688458225E-008, 7.89724635377817962E-005, 4.61868720707717372E-005, 3.55682891687782599E-005 },
600 { -9.20788748896973282E-008, -5.38521280004949642E-008, -4.14712660805579618E-008, 7.78626692360739821E-005, 4.55378113077967091E-005, 3.50684505810897702E-005 },
601 { 1.85082436324531617E-008, 1.20506219457886855E-007, -8.31277842285972640E-008, 1.27364008345789645E-004, -1.54770720974742483E-004, -1.78589436862677754E-004 }
602 };
603
604 PVCoordinates pv = orbCir.getPVCoordinates();
605 Assert.assertEquals(0, pv.getPosition().subtract(pRef).getNorm(), 3.0e-16 * pRef.getNorm());
606 Assert.assertEquals(0, pv.getVelocity().subtract(vRef).getNorm(), 2.0e-16 * vRef.getNorm());
607
608 double[][] jacobian = new double[6][6];
609 orbCir.getJacobianWrtCartesian(PositionAngle.MEAN, jacobian);
610
611 for (int i = 0; i < jacobian.length; i++) {
612 double[] row = jacobian[i];
613 double[] rowRef = jRef[i];
614 for (int j = 0; j < row.length; j++) {
615 Assert.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 5.0e-15);
616 }
617 }
618
619 }
620
621 @Test
622 public void testJacobianFinitedifferences() {
623
624 AbsoluteDate dateTca = new AbsoluteDate(2000, 04, 01, 0, 0, 0.000, TimeScalesFactory.getUTC());
625 double mu = 3.986004415e+14;
626 CircularOrbit orbCir = new CircularOrbit(7000000.0, 0.01, -0.02, 1.2, 2.1,
627 0.7, PositionAngle.MEAN,
628 FramesFactory.getEME2000(), dateTca, mu);
629
630 for (PositionAngle type : PositionAngle.values()) {
631 double hP = 2.0;
632 double[][] finiteDiffJacobian = finiteDifferencesJacobian(type, orbCir, hP);
633 double[][] jacobian = new double[6][6];
634 orbCir.getJacobianWrtCartesian(type, jacobian);
635
636 for (int i = 0; i < jacobian.length; i++) {
637 double[] row = jacobian[i];
638 double[] rowRef = finiteDiffJacobian[i];
639 for (int j = 0; j < row.length; j++) {
640 Assert.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 8.0e-9);
641 }
642 }
643
644 double[][] invJacobian = new double[6][6];
645 orbCir.getJacobianWrtParameters(type, invJacobian);
646 MatrixUtils.createRealMatrix(jacobian).
647 multiply(MatrixUtils.createRealMatrix(invJacobian)).
648 walkInRowOrder(new RealMatrixPreservingVisitor() {
649 public void start(int rows, int columns,
650 int startRow, int endRow, int startColumn, int endColumn) {
651 }
652
653 public void visit(int row, int column, double value) {
654 Assert.assertEquals(row == column ? 1.0 : 0.0, value, 4.0e-9);
655 }
656
657 public double end() {
658 return Double.NaN;
659 }
660 });
661
662 }
663
664 }
665
666 private double[][] finiteDifferencesJacobian(PositionAngle type, CircularOrbit orbit, double hP)
667 {
668 double[][] jacobian = new double[6][6];
669 for (int i = 0; i < 6; ++i) {
670 fillColumn(type, i, orbit, hP, jacobian);
671 }
672 return jacobian;
673 }
674
675 private void fillColumn(PositionAngle type, int i, CircularOrbit orbit, double hP, double[][] jacobian) {
676
677
678
679 Vector3D p = orbit.getPVCoordinates().getPosition();
680 Vector3D v = orbit.getPVCoordinates().getVelocity();
681 double hV = orbit.getMu() * hP / (v.getNorm() * p.getNormSq());
682
683 double h;
684 Vector3D dP = Vector3D.ZERO;
685 Vector3D dV = Vector3D.ZERO;
686 switch (i) {
687 case 0:
688 h = hP;
689 dP = new Vector3D(hP, 0, 0);
690 break;
691 case 1:
692 h = hP;
693 dP = new Vector3D(0, hP, 0);
694 break;
695 case 2:
696 h = hP;
697 dP = new Vector3D(0, 0, hP);
698 break;
699 case 3:
700 h = hV;
701 dV = new Vector3D(hV, 0, 0);
702 break;
703 case 4:
704 h = hV;
705 dV = new Vector3D(0, hV, 0);
706 break;
707 default:
708 h = hV;
709 dV = new Vector3D(0, 0, hV);
710 break;
711 }
712
713 CircularOrbit oM4h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, -4, dP), new Vector3D(1, v, -4, dV)),
714 orbit.getFrame(), orbit.getDate(), orbit.getMu());
715 CircularOrbit oM3h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, -3, dP), new Vector3D(1, v, -3, dV)),
716 orbit.getFrame(), orbit.getDate(), orbit.getMu());
717 CircularOrbit oM2h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, -2, dP), new Vector3D(1, v, -2, dV)),
718 orbit.getFrame(), orbit.getDate(), orbit.getMu());
719 CircularOrbit oM1h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, -1, dP), new Vector3D(1, v, -1, dV)),
720 orbit.getFrame(), orbit.getDate(), orbit.getMu());
721 CircularOrbit oP1h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, +1, dP), new Vector3D(1, v, +1, dV)),
722 orbit.getFrame(), orbit.getDate(), orbit.getMu());
723 CircularOrbit oP2h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, +2, dP), new Vector3D(1, v, +2, dV)),
724 orbit.getFrame(), orbit.getDate(), orbit.getMu());
725 CircularOrbit oP3h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, +3, dP), new Vector3D(1, v, +3, dV)),
726 orbit.getFrame(), orbit.getDate(), orbit.getMu());
727 CircularOrbit oP4h = new CircularOrbit(new PVCoordinates(new Vector3D(1, p, +4, dP), new Vector3D(1, v, +4, dV)),
728 orbit.getFrame(), orbit.getDate(), orbit.getMu());
729
730 jacobian[0][i] = (-3 * (oP4h.getA() - oM4h.getA()) +
731 32 * (oP3h.getA() - oM3h.getA()) -
732 168 * (oP2h.getA() - oM2h.getA()) +
733 672 * (oP1h.getA() - oM1h.getA())) / (840 * h);
734 jacobian[1][i] = (-3 * (oP4h.getCircularEx() - oM4h.getCircularEx()) +
735 32 * (oP3h.getCircularEx() - oM3h.getCircularEx()) -
736 168 * (oP2h.getCircularEx() - oM2h.getCircularEx()) +
737 672 * (oP1h.getCircularEx() - oM1h.getCircularEx())) / (840 * h);
738 jacobian[2][i] = (-3 * (oP4h.getCircularEy() - oM4h.getCircularEy()) +
739 32 * (oP3h.getCircularEy() - oM3h.getCircularEy()) -
740 168 * (oP2h.getCircularEy() - oM2h.getCircularEy()) +
741 672 * (oP1h.getCircularEy() - oM1h.getCircularEy())) / (840 * h);
742 jacobian[3][i] = (-3 * (oP4h.getI() - oM4h.getI()) +
743 32 * (oP3h.getI() - oM3h.getI()) -
744 168 * (oP2h.getI() - oM2h.getI()) +
745 672 * (oP1h.getI() - oM1h.getI())) / (840 * h);
746 jacobian[4][i] = (-3 * (oP4h.getRightAscensionOfAscendingNode() - oM4h.getRightAscensionOfAscendingNode()) +
747 32 * (oP3h.getRightAscensionOfAscendingNode() - oM3h.getRightAscensionOfAscendingNode()) -
748 168 * (oP2h.getRightAscensionOfAscendingNode() - oM2h.getRightAscensionOfAscendingNode()) +
749 672 * (oP1h.getRightAscensionOfAscendingNode() - oM1h.getRightAscensionOfAscendingNode())) / (840 * h);
750 jacobian[5][i] = (-3 * (oP4h.getAlpha(type) - oM4h.getAlpha(type)) +
751 32 * (oP3h.getAlpha(type) - oM3h.getAlpha(type)) -
752 168 * (oP2h.getAlpha(type) - oM2h.getAlpha(type)) +
753 672 * (oP1h.getAlpha(type) - oM1h.getAlpha(type))) / (840 * h);
754
755 }
756
757 @Test
758 public void testInterpolationWithDerivatives() {
759 doTestInterpolation(true,
760 397, 1.88e-8,
761 610, 3.52e-6,
762 4870, 115);
763 }
764
765 @Test
766 public void testInterpolationWithoutDerivatives() {
767 doTestInterpolation(false,
768 397, 0.0372,
769 610.0, 1.23,
770 4870, 8869);
771 }
772
773 private void doTestInterpolation(boolean useDerivatives,
774 double shiftErrorWithin, double interpolationErrorWithin,
775 double shiftErrorSlightlyPast, double interpolationErrorSlightlyPast,
776 double shiftErrorFarPast, double interpolationErrorFarPast)
777 {
778
779 final double ehMu = 3.9860047e14;
780 final double ae = 6.378137e6;
781 final double c20 = -1.08263e-3;
782 final double c30 = 2.54e-6;
783 final double c40 = 1.62e-6;
784 final double c50 = 2.3e-7;
785 final double c60 = -5.5e-7;
786
787 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
788 final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
789 final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
790 final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity),
791 FramesFactory.getEME2000(), date, ehMu);
792
793 EcksteinHechlerPropagator propagator =
794 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
795
796
797 List<Orbit> sample = new ArrayList<Orbit>();
798 for (double dt = 0; dt < 300.0; dt += 60.0) {
799 Orbit orbit = propagator.propagate(date.shiftedBy(dt)).getOrbit();
800 if (!useDerivatives) {
801
802 double[] stateVector = new double[6];
803 orbit.getType().mapOrbitToArray(orbit, PositionAngle.TRUE, stateVector, null);
804 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngle.TRUE,
805 orbit.getDate(), orbit.getMu(), orbit.getFrame());
806 }
807 sample.add(orbit);
808 }
809
810
811 double maxShiftError = 0;
812 double maxInterpolationError = 0;
813 for (double dt = 0; dt < 241.0; dt += 1.0) {
814 AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
815 Vector3D shifted = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
816 Vector3D interpolated = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
817 Vector3D propagated = propagator.propagate(t).getPVCoordinates().getPosition();
818 maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
819 maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
820 }
821 Assert.assertEquals(shiftErrorWithin, maxShiftError, 0.01 * shiftErrorWithin);
822 Assert.assertEquals(interpolationErrorWithin, maxInterpolationError, 0.01 * interpolationErrorWithin);
823
824
825 maxShiftError = 0;
826 maxInterpolationError = 0;
827 for (double dt = 240; dt < 300.0; dt += 1.0) {
828 AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
829 Vector3D shifted = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
830 Vector3D interpolated = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
831 Vector3D propagated = propagator.propagate(t).getPVCoordinates().getPosition();
832 maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
833 maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
834 }
835 Assert.assertEquals(shiftErrorSlightlyPast, maxShiftError, 0.01 * shiftErrorSlightlyPast);
836 Assert.assertEquals(interpolationErrorSlightlyPast, maxInterpolationError, 0.01 * interpolationErrorSlightlyPast);
837
838
839 maxShiftError = 0;
840 maxInterpolationError = 0;
841 for (double dt = 300; dt < 1000; dt += 1.0) {
842 AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
843 Vector3D shifted = initialOrbit.shiftedBy(dt).getPVCoordinates().getPosition();
844 Vector3D interpolated = initialOrbit.interpolate(t, sample).getPVCoordinates().getPosition();
845 Vector3D propagated = propagator.propagate(t).getPVCoordinates().getPosition();
846 maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
847 maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
848 }
849 Assert.assertEquals(shiftErrorFarPast, maxShiftError, 0.01 * shiftErrorFarPast);
850 Assert.assertEquals(interpolationErrorFarPast, maxInterpolationError, 0.01 * interpolationErrorFarPast);
851
852 }
853
854 @Test
855 public void testSerialization()
856 throws IOException, ClassNotFoundException {
857 Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
858 Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
859 PVCoordinates pvCoordinates = new PVCoordinates( position, velocity);
860 CircularOrbit orbit = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
861 Assert.assertEquals(42255170.003, orbit.getA(), 1.0e-3);
862
863 ByteArrayOutputStream bos = new ByteArrayOutputStream();
864 ObjectOutputStream oos = new ObjectOutputStream(bos);
865 oos.writeObject(orbit);
866
867 Assert.assertTrue(bos.size() > 350);
868 Assert.assertTrue(bos.size() < 400);
869
870 ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
871 ObjectInputStream ois = new ObjectInputStream(bis);
872 CircularOrbit deserialized = (CircularOrbit) ois.readObject();
873 Assert.assertEquals(orbit.getA(), deserialized.getA(), 1.0e-10);
874 Assert.assertEquals(orbit.getCircularEx(), deserialized.getCircularEx(), 1.0e-10);
875 Assert.assertEquals(orbit.getCircularEy(), deserialized.getCircularEy(), 1.0e-10);
876 Assert.assertEquals(orbit.getI(), deserialized.getI(), 1.0e-10);
877 Assert.assertEquals(orbit.getRightAscensionOfAscendingNode(), deserialized.getRightAscensionOfAscendingNode(), 1.0e-10);
878 Assert.assertEquals(orbit.getAlphaV(), deserialized.getAlphaV(), 1.0e-10);
879 Assert.assertTrue(Double.isNaN(orbit.getADot()) && Double.isNaN(deserialized.getADot()));
880 Assert.assertTrue(Double.isNaN(orbit.getCircularExDot()) && Double.isNaN(deserialized.getCircularExDot()));
881 Assert.assertTrue(Double.isNaN(orbit.getCircularEyDot()) && Double.isNaN(deserialized.getCircularEyDot()));
882 Assert.assertTrue(Double.isNaN(orbit.getIDot()) && Double.isNaN(deserialized.getIDot()));
883 Assert.assertTrue(Double.isNaN(orbit.getRightAscensionOfAscendingNodeDot()) && Double.isNaN(deserialized.getRightAscensionOfAscendingNodeDot()));
884 Assert.assertTrue(Double.isNaN(orbit.getAlphaVDot()) && Double.isNaN(deserialized.getAlphaVDot()));
885 Assert.assertEquals(orbit.getDate(), deserialized.getDate());
886 Assert.assertEquals(orbit.getMu(), deserialized.getMu(), 1.0e-10);
887 Assert.assertEquals(orbit.getFrame().getName(), deserialized.getFrame().getName());
888
889 }
890
891 @Test
892 public void testSerializationWithDerivatives()
893 throws IOException, ClassNotFoundException {
894 Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
895 Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
896 double r2 = position.getNormSq();
897 double r = FastMath.sqrt(r2);
898 Vector3D acceleration = new Vector3D(-mu / (r * r2), position,
899 1, new Vector3D(-0.1, 0.2, 0.3));
900 PVCoordinates pvCoordinates = new PVCoordinates( position, velocity, acceleration);
901 CircularOrbit orbit = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
902 Assert.assertEquals(42255170.003, orbit.getA(), 1.0e-3);
903
904 ByteArrayOutputStream bos = new ByteArrayOutputStream();
905 ObjectOutputStream oos = new ObjectOutputStream(bos);
906 oos.writeObject(orbit);
907
908 Assert.assertTrue(bos.size() > 400);
909 Assert.assertTrue(bos.size() < 450);
910
911 ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
912 ObjectInputStream ois = new ObjectInputStream(bis);
913 CircularOrbit deserialized = (CircularOrbit) ois.readObject();
914 Assert.assertEquals(orbit.getA(), deserialized.getA(), 1.0e-10);
915 Assert.assertEquals(orbit.getCircularEx(), deserialized.getCircularEx(), 1.0e-10);
916 Assert.assertEquals(orbit.getCircularEy(), deserialized.getCircularEy(), 1.0e-10);
917 Assert.assertEquals(orbit.getI(), deserialized.getI(), 1.0e-10);
918 Assert.assertEquals(orbit.getRightAscensionOfAscendingNode(), deserialized.getRightAscensionOfAscendingNode(), 1.0e-10);
919 Assert.assertEquals(orbit.getAlphaV(), deserialized.getAlphaV(), 1.0e-10);
920 Assert.assertEquals(orbit.getADot(), deserialized.getADot(), 1.0e-10);
921 Assert.assertEquals(orbit.getCircularExDot(), deserialized.getCircularExDot(), 1.0e-10);
922 Assert.assertEquals(orbit.getCircularEyDot(), deserialized.getCircularEyDot(), 1.0e-10);
923 Assert.assertEquals(orbit.getIDot(), deserialized.getIDot(), 1.0e-10);
924 Assert.assertEquals(orbit.getRightAscensionOfAscendingNodeDot(), deserialized.getRightAscensionOfAscendingNodeDot(), 1.0e-10);
925 Assert.assertEquals(orbit.getAlphaVDot(), deserialized.getAlphaVDot(), 1.0e-10);
926 Assert.assertEquals(orbit.getDate(), deserialized.getDate());
927 Assert.assertEquals(orbit.getMu(), deserialized.getMu(), 1.0e-10);
928 Assert.assertEquals(orbit.getFrame().getName(), deserialized.getFrame().getName());
929
930 }
931
932 @Test
933 public void testSerializationNoPVWithDerivatives()
934 throws IOException, ClassNotFoundException {
935 Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
936 Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
937 double r2 = position.getNormSq();
938 double r = FastMath.sqrt(r2);
939 Vector3D acceleration = new Vector3D(-mu / (r * r2), position,
940 1, new Vector3D(-0.1, 0.2, 0.3));
941 PVCoordinates pvCoordinates = new PVCoordinates( position, velocity, acceleration);
942 CircularOrbit original = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
943
944
945
946 CircularOrbit orbit = new CircularOrbit(original.getA(), original.getCircularEx(), original.getCircularEy(),
947 original.getI(), original.getRightAscensionOfAscendingNode(),
948 original.getAlphaV(),
949 original.getADot(), original.getCircularExDot(), original.getCircularEyDot(),
950 original.getIDot(), original.getRightAscensionOfAscendingNodeDot(),
951 original.getAlphaVDot(),
952 PositionAngle.TRUE, original.getFrame(),
953 original.getDate(), original.getMu());
954 Assert.assertEquals(42255170.003, orbit.getA(), 1.0e-3);
955
956 ByteArrayOutputStream bos = new ByteArrayOutputStream();
957 ObjectOutputStream oos = new ObjectOutputStream(bos);
958 oos.writeObject(orbit);
959
960 Assert.assertTrue(bos.size() > 330);
961 Assert.assertTrue(bos.size() < 380);
962
963 ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
964 ObjectInputStream ois = new ObjectInputStream(bis);
965 CircularOrbit deserialized = (CircularOrbit) ois.readObject();
966 Assert.assertEquals(orbit.getA(), deserialized.getA(), 1.0e-10);
967 Assert.assertEquals(orbit.getCircularEx(), deserialized.getCircularEx(), 1.0e-10);
968 Assert.assertEquals(orbit.getCircularEy(), deserialized.getCircularEy(), 1.0e-10);
969 Assert.assertEquals(orbit.getI(), deserialized.getI(), 1.0e-10);
970 Assert.assertEquals(orbit.getRightAscensionOfAscendingNode(), deserialized.getRightAscensionOfAscendingNode(), 1.0e-10);
971 Assert.assertEquals(orbit.getAlphaV(), deserialized.getAlphaV(), 1.0e-10);
972 Assert.assertEquals(orbit.getADot(), deserialized.getADot(), 1.0e-10);
973 Assert.assertEquals(orbit.getCircularExDot(), deserialized.getCircularExDot(), 1.0e-10);
974 Assert.assertEquals(orbit.getCircularEyDot(), deserialized.getCircularEyDot(), 1.0e-10);
975 Assert.assertEquals(orbit.getIDot(), deserialized.getIDot(), 1.0e-10);
976 Assert.assertEquals(orbit.getRightAscensionOfAscendingNodeDot(), deserialized.getRightAscensionOfAscendingNodeDot(), 1.0e-10);
977 Assert.assertEquals(orbit.getAlphaVDot(), deserialized.getAlphaVDot(), 1.0e-10);
978 Assert.assertEquals(orbit.getDate(), deserialized.getDate());
979 Assert.assertEquals(orbit.getMu(), deserialized.getMu(), 1.0e-10);
980 Assert.assertEquals(orbit.getFrame().getName(), deserialized.getFrame().getName());
981
982 }
983
984 @Test
985 public void testNonKeplerianDerivatives() {
986 final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
987 final Vector3D position = new Vector3D(6896874.444705, 1956581.072644, -147476.245054);
988 final Vector3D velocity = new Vector3D(166.816407662, -1106.783301861, -7372.745712770);
989 final Vector3D acceleration = new Vector3D(-7.466182457944, -2.118153357345, 0.160004048437);
990 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity, acceleration);
991 final Frame frame = FramesFactory.getEME2000();
992 final double mu = Constants.EIGEN5C_EARTH_MU;
993 final CircularOrbit orbit = new CircularOrbit(pv, frame, mu);
994
995 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getA()),
996 orbit.getADot(),
997 4.3e-8);
998 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getEquinoctialEx()),
999 orbit.getEquinoctialExDot(),
1000 2.1e-15);
1001 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getEquinoctialEy()),
1002 orbit.getEquinoctialEyDot(),
1003 5.4e-16);
1004 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getHx()),
1005 orbit.getHxDot(),
1006 1.6e-15);
1007 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getHy()),
1008 orbit.getHyDot(),
1009 7.3e-17);
1010 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLv()),
1011 orbit.getLvDot(),
1012 3.4e-16);
1013 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLE()),
1014 orbit.getLEDot(),
1015 3.5e-15);
1016 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLM()),
1017 orbit.getLMDot(),
1018 5.3e-15);
1019 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getE()),
1020 orbit.getEDot(),
1021 6.8e-16);
1022 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getI()),
1023 orbit.getIDot(),
1024 5.7e-16);
1025 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getCircularEx()),
1026 orbit.getCircularExDot(),
1027 2.2e-15);
1028 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getCircularEy()),
1029 orbit.getCircularEyDot(),
1030 5.3e-17);
1031 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAlphaV()),
1032 orbit.getAlphaVDot(),
1033 4.3e-15);
1034 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAlphaE()),
1035 orbit.getAlphaEDot(),
1036 1.2e-15);
1037 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAlphaM()),
1038 orbit.getAlphaMDot(),
1039 3.7e-15);
1040 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAlpha(PositionAngle.TRUE)),
1041 orbit.getAlphaDot(PositionAngle.TRUE),
1042 4.3e-15);
1043 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAlpha(PositionAngle.ECCENTRIC)),
1044 orbit.getAlphaDot(PositionAngle.ECCENTRIC),
1045 1.2e-15);
1046 Assert.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getAlpha(PositionAngle.MEAN)),
1047 orbit.getAlphaDot(PositionAngle.MEAN),
1048 3.7e-15);
1049
1050 }
1051
1052 private <S extends Function<CircularOrbit, Double>>
1053 double differentiate(TimeStampedPVCoordinates pv, Frame frame, double mu, S picker) {
1054 final DSFactory factory = new DSFactory(1, 1);
1055 FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(8, 0.1);
1056 UnivariateDifferentiableFunction diff = differentiator.differentiate(new UnivariateFunction() {
1057 public double value(double dt) {
1058 return picker.apply(new CircularOrbit(pv.shiftedBy(dt), frame, mu));
1059 }
1060 });
1061 return diff.value(factory.variable(0, 0.0)).getPartialDerivative(1);
1062 }
1063
1064 @Test
1065 public void testPositionAngleDerivatives() {
1066 final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1067 final Vector3D position = new Vector3D(6896874.444705, 1956581.072644, -147476.245054);
1068 final Vector3D velocity = new Vector3D(166.816407662, -1106.783301861, -7372.745712770);
1069 final Vector3D acceleration = new Vector3D(-7.466182457944, -2.118153357345, 0.160004048437);
1070 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity, acceleration);
1071 final Frame frame = FramesFactory.getEME2000();
1072 final double mu = Constants.EIGEN5C_EARTH_MU;
1073 final CircularOrbit orbit = new CircularOrbit(pv, frame, mu);
1074
1075 for (PositionAngle type : PositionAngle.values()) {
1076 final CircularOrbit rebuilt = new CircularOrbit(orbit.getA(),
1077 orbit.getCircularEx(),
1078 orbit.getCircularEy(),
1079 orbit.getI(),
1080 orbit.getRightAscensionOfAscendingNode(),
1081 orbit.getAlpha(type),
1082 orbit.getADot(),
1083 orbit.getCircularExDot(),
1084 orbit.getCircularEyDot(),
1085 orbit.getIDot(),
1086 orbit.getRightAscensionOfAscendingNodeDot(),
1087 orbit.getAlphaDot(type),
1088 type, orbit.getFrame(), orbit.getDate(), orbit.getMu());
1089 MatcherAssert.assertThat(rebuilt.getA(), relativelyCloseTo(orbit.getA(), 1));
1090 MatcherAssert.assertThat(rebuilt.getCircularEx(), relativelyCloseTo(orbit.getCircularEx(), 1));
1091 MatcherAssert.assertThat(rebuilt.getCircularEy(), relativelyCloseTo(orbit.getCircularEy(), 1));
1092 MatcherAssert.assertThat(rebuilt.getE(), relativelyCloseTo(orbit.getE(), 1));
1093 MatcherAssert.assertThat(rebuilt.getI(), relativelyCloseTo(orbit.getI(), 1));
1094 MatcherAssert.assertThat(rebuilt.getRightAscensionOfAscendingNode(), relativelyCloseTo(orbit.getRightAscensionOfAscendingNode(), 1));
1095 MatcherAssert.assertThat(rebuilt.getADot(), relativelyCloseTo(orbit.getADot(), 1));
1096 MatcherAssert.assertThat(rebuilt.getCircularExDot(), relativelyCloseTo(orbit.getCircularExDot(), 1));
1097 MatcherAssert.assertThat(rebuilt.getCircularEyDot(), relativelyCloseTo(orbit.getCircularEyDot(), 1));
1098 MatcherAssert.assertThat(rebuilt.getEDot(), relativelyCloseTo(orbit.getEDot(), 1));
1099 MatcherAssert.assertThat(rebuilt.getIDot(), relativelyCloseTo(orbit.getIDot(), 1));
1100 MatcherAssert.assertThat(rebuilt.getRightAscensionOfAscendingNodeDot(), relativelyCloseTo(orbit.getRightAscensionOfAscendingNodeDot(), 1));
1101 for (PositionAngle type2 : PositionAngle.values()) {
1102 MatcherAssert.assertThat(rebuilt.getAlpha(type2), relativelyCloseTo(orbit.getAlpha(type2), 1));
1103 MatcherAssert.assertThat(rebuilt.getAlphaDot(type2), relativelyCloseTo(orbit.getAlphaDot(type2), 1));
1104 }
1105 }
1106
1107 }
1108
1109 @Test
1110 public void testEquatorialRetrograde() {
1111 Vector3D position = new Vector3D(10000000.0, 0.0, 0.0);
1112 Vector3D velocity = new Vector3D(0.0, -6500.0, 1.0e-10);
1113 double r2 = position.getNormSq();
1114 double r = FastMath.sqrt(r2);
1115 Vector3D acceleration = new Vector3D(-mu / (r * r2), position,
1116 1, new Vector3D(-0.1, 0.2, 0.3));
1117 PVCoordinates pvCoordinates = new PVCoordinates(position, velocity, acceleration);
1118 CircularOrbit orbit = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
1119 Assert.assertEquals(10637829.465, orbit.getA(), 1.0e-3);
1120 Assert.assertEquals(-738.145, orbit.getADot(), 1.0e-3);
1121 Assert.assertEquals(0.05995861, orbit.getE(), 1.0e-8);
1122 Assert.assertEquals(-6.523e-5, orbit.getEDot(), 1.0e-8);
1123 Assert.assertEquals(FastMath.PI, orbit.getI(), 2.0e-14);
1124 Assert.assertEquals(-4.615e-5, orbit.getIDot(), 1.0e-8);
1125 Assert.assertTrue(Double.isNaN(orbit.getHx()));
1126 Assert.assertTrue(Double.isNaN(orbit.getHxDot()));
1127 Assert.assertTrue(Double.isNaN(orbit.getHy()));
1128 Assert.assertTrue(Double.isNaN(orbit.getHyDot()));
1129 }
1130
1131 @Test
1132 public void testDerivativesConversionSymmetry() {
1133 final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:01:20.000", TimeScalesFactory.getUTC());
1134 Vector3D position = new Vector3D(6893443.400234382, 1886406.1073757345, -589265.1150359757);
1135 Vector3D velocity = new Vector3D(-281.1261461082365, -1231.6165642450928, -7348.756363469432);
1136 Vector3D acceleration = new Vector3D(-7.460341170581685, -2.0415957334584527, 0.6393322823627762);
1137 PVCoordinates pvCoordinates = new PVCoordinates( position, velocity, acceleration);
1138 CircularOrbit orbit = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(),
1139 date, Constants.EIGEN5C_EARTH_MU);
1140 Assert.assertTrue(orbit.hasDerivatives());
1141 double r2 = position.getNormSq();
1142 double r = FastMath.sqrt(r2);
1143 Vector3D keplerianAcceleration = new Vector3D(-orbit.getMu() / (r2 * r), position);
1144 Assert.assertEquals(0.0101, Vector3D.distance(keplerianAcceleration, acceleration), 1.0e-4);
1145
1146 for (OrbitType type : OrbitType.values()) {
1147 Orbit converted = type.convertType(orbit);
1148 Assert.assertTrue(converted.hasDerivatives());
1149 CircularOrbit rebuilt = (CircularOrbit) OrbitType.CIRCULAR.convertType(converted);
1150 Assert.assertTrue(rebuilt.hasDerivatives());
1151 Assert.assertEquals(orbit.getADot(), rebuilt.getADot(), 3.0e-13);
1152 Assert.assertEquals(orbit.getCircularExDot(), rebuilt.getCircularExDot(), 1.0e-15);
1153 Assert.assertEquals(orbit.getCircularEyDot(), rebuilt.getCircularEyDot(), 1.0e-15);
1154 Assert.assertEquals(orbit.getIDot(), rebuilt.getIDot(), 1.0e-15);
1155 Assert.assertEquals(orbit.getRightAscensionOfAscendingNodeDot(), rebuilt.getRightAscensionOfAscendingNodeDot(), 1.0e-15);
1156 Assert.assertEquals(orbit.getAlphaVDot(), rebuilt.getAlphaVDot(), 1.0e-15);
1157 }
1158
1159 }
1160
1161 @Test
1162 public void testToString() {
1163 Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
1164 Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
1165 PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
1166 CircularOrbit orbit = new CircularOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
1167 Assert.assertEquals("circular parameters: {a: 4.225517000282565E7, ex: 0.002082917137146049, ey: 5.173980074371024E-4, i: 0.20189257051515358, raan: -87.91788415673473, alphaV: -137.84099636616548;}",
1168 orbit.toString());
1169 }
1170
1171 @Test
1172 public void testCopyNonKeplerianAcceleration() {
1173
1174 final Frame eme2000 = FramesFactory.getEME2000();
1175
1176
1177 final Vector3D position = new Vector3D(42164140, 0, 0);
1178
1179 final PVCoordinates pv = new PVCoordinates(position,
1180 new Vector3D(0, FastMath.sqrt(mu / position.getNorm()), 0));
1181
1182 final Orbit orbit = new CircularOrbit(pv, eme2000, date, mu);
1183
1184
1185 final Orbit orbitCopy = new CircularOrbit(orbit);
1186
1187
1188 final Orbit shiftedOrbit = orbit.shiftedBy(10);
1189 final Orbit shiftedOrbitCopy = orbitCopy.shiftedBy(10);
1190
1191 Assert.assertEquals(0.0,
1192 Vector3D.distance(shiftedOrbit.getPVCoordinates().getPosition(),
1193 shiftedOrbitCopy.getPVCoordinates().getPosition()),
1194 1.0e-10);
1195 Assert.assertEquals(0.0,
1196 Vector3D.distance(shiftedOrbit.getPVCoordinates().getVelocity(),
1197 shiftedOrbitCopy.getPVCoordinates().getVelocity()),
1198 1.0e-10);
1199
1200 }
1201
1202 @Test
1203 public void testNormalize() {
1204 CircularOrbit withoutDerivatives =
1205 new CircularOrbit(42166712.0, 0.005, -0.025, 1.6,
1206 1.25, 0.4, PositionAngle.MEAN,
1207 FramesFactory.getEME2000(), date, mu);
1208 CircularOrbit ref =
1209 new CircularOrbit(24000000.0, -0.012, 0.01, 0.2,
1210 -6.28, 6.28, PositionAngle.MEAN,
1211 FramesFactory.getEME2000(), date, mu);
1212
1213 CircularOrbit normalized1 = (CircularOrbit) OrbitType.CIRCULAR.normalize(withoutDerivatives, ref);
1214 Assert.assertFalse(normalized1.hasDerivatives());
1215 Assert.assertEquals(0.0, normalized1.getA() - withoutDerivatives.getA(), 1.0e-6);
1216 Assert.assertEquals(0.0, normalized1.getCircularEx() - withoutDerivatives.getCircularEx(), 1.0e-10);
1217 Assert.assertEquals(0.0, normalized1.getCircularEy() - withoutDerivatives.getCircularEy(), 1.0e-10);
1218 Assert.assertEquals(0.0, normalized1.getI() - withoutDerivatives.getI(), 1.0e-10);
1219 Assert.assertEquals(-MathUtils.TWO_PI, normalized1.getRightAscensionOfAscendingNode() - withoutDerivatives.getRightAscensionOfAscendingNode(), 1.0e-10);
1220 Assert.assertEquals(+MathUtils.TWO_PI, normalized1.getAlphaV() - withoutDerivatives.getAlphaV(), 1.0e-10);
1221 Assert.assertTrue(Double.isNaN(normalized1.getADot()));
1222 Assert.assertTrue(Double.isNaN(normalized1.getCircularExDot()));
1223 Assert.assertTrue(Double.isNaN(normalized1.getCircularEyDot()));
1224 Assert.assertTrue(Double.isNaN(normalized1.getIDot()));
1225 Assert.assertTrue(Double.isNaN(normalized1.getRightAscensionOfAscendingNodeDot()));
1226 Assert.assertTrue(Double.isNaN(normalized1.getAlphaVDot()));
1227
1228 double[] p = new double[6];
1229 OrbitType.CIRCULAR.mapOrbitToArray(withoutDerivatives, PositionAngle.TRUE, p, null);
1230 CircularOrbit withDerivatives = (CircularOrbit) OrbitType.CIRCULAR.mapArrayToOrbit(p,
1231 new double[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 },
1232 PositionAngle.TRUE,
1233 withoutDerivatives.getDate(),
1234 withoutDerivatives.getMu(),
1235 withoutDerivatives.getFrame());
1236 CircularOrbit normalized2 = (CircularOrbit) OrbitType.CIRCULAR.normalize(withDerivatives, ref);
1237 Assert.assertTrue(normalized2.hasDerivatives());
1238 Assert.assertEquals(0.0, normalized2.getA() - withDerivatives.getA(), 1.0e-6);
1239 Assert.assertEquals(0.0, normalized2.getCircularEx() - withDerivatives.getCircularEx(), 1.0e-10);
1240 Assert.assertEquals(0.0, normalized2.getCircularEy() - withDerivatives.getCircularEy(), 1.0e-10);
1241 Assert.assertEquals(0.0, normalized2.getI() - withDerivatives.getI(), 1.0e-10);
1242 Assert.assertEquals(-MathUtils.TWO_PI, normalized2.getRightAscensionOfAscendingNode() - withDerivatives.getRightAscensionOfAscendingNode(), 1.0e-10);
1243 Assert.assertEquals(+MathUtils.TWO_PI, normalized2.getAlphaV() - withDerivatives.getAlphaV(), 1.0e-10);
1244 Assert.assertEquals(0.0, normalized2.getADot() - withDerivatives.getADot(), 1.0e-10);
1245 Assert.assertEquals(0.0, normalized2.getCircularExDot() - withDerivatives.getCircularExDot(), 1.0e-10);
1246 Assert.assertEquals(0.0, normalized2.getCircularEyDot() - withDerivatives.getCircularEyDot(), 1.0e-10);
1247 Assert.assertEquals(0.0, normalized2.getIDot() - withDerivatives.getIDot(), 1.0e-10);
1248 Assert.assertEquals(0.0, normalized2.getRightAscensionOfAscendingNodeDot() - withDerivatives.getRightAscensionOfAscendingNodeDot(), 1.0e-10);
1249 Assert.assertEquals(0.0, normalized2.getAlphaVDot() - withDerivatives.getAlphaVDot(), 1.0e-10);
1250
1251 }
1252
1253 @Before
1254 public void setUp() {
1255
1256 Utils.setDataRoot("regular-data");
1257
1258
1259 date = AbsoluteDate.J2000_EPOCH;
1260
1261
1262 mu = 3.9860047e14;
1263 }
1264
1265 @After
1266 public void tearDown() {
1267 date = null;
1268 }
1269
1270 }