1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.orbits;
18
19 import java.io.Serializable;
20 import java.util.List;
21 import java.util.stream.Collectors;
22 import java.util.stream.Stream;
23
24 import org.hipparchus.analysis.differentiation.DSFactory;
25 import org.hipparchus.analysis.differentiation.DerivativeStructure;
26 import org.hipparchus.analysis.interpolation.HermiteInterpolator;
27 import org.hipparchus.exception.MathIllegalStateException;
28 import org.hipparchus.geometry.euclidean.threed.Vector3D;
29 import org.hipparchus.util.FastMath;
30 import org.hipparchus.util.MathUtils;
31 import org.orekit.errors.OrekitIllegalArgumentException;
32 import org.orekit.errors.OrekitInternalError;
33 import org.orekit.errors.OrekitMessages;
34 import org.orekit.frames.Frame;
35 import org.orekit.time.AbsoluteDate;
36 import org.orekit.utils.PVCoordinates;
37 import org.orekit.utils.TimeStampedPVCoordinates;
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79 public class KeplerianOrbit extends Orbit {
80
81
82 private static final long serialVersionUID = 20170414L;
83
84
85 private static final DSFactory FACTORY = new DSFactory(1, 1);
86
87
88 private static final double A;
89
90
91 private static final double B;
92
93 static {
94 final double k1 = 3 * FastMath.PI + 2;
95 final double k2 = FastMath.PI - 1;
96 final double k3 = 6 * FastMath.PI - 1;
97 A = 3 * k2 * k2 / k1;
98 B = k3 * k3 / (6 * k1);
99 }
100
101
102 private final double a;
103
104
105 private final double e;
106
107
108 private final double i;
109
110
111 private final double pa;
112
113
114 private final double raan;
115
116
117 private final double v;
118
119
120 private final double aDot;
121
122
123 private final double eDot;
124
125
126 private final double iDot;
127
128
129 private final double paDot;
130
131
132 private final double raanDot;
133
134
135 private final double vDot;
136
137
138 private transient PVCoordinates partialPV;
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156 public KeplerianOrbit(final double a, final double e, final double i,
157 final double pa, final double raan, final double anomaly,
158 final PositionAngle type,
159 final Frame frame, final AbsoluteDate date, final double mu)
160 throws IllegalArgumentException {
161 this(a, e, i, pa, raan, anomaly,
162 Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN,
163 type, frame, date, mu);
164 }
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189 public KeplerianOrbit(final double a, final double e, final double i,
190 final double pa, final double raan, final double anomaly,
191 final double aDot, final double eDot, final double iDot,
192 final double paDot, final double raanDot, final double anomalyDot,
193 final PositionAngle type,
194 final Frame frame, final AbsoluteDate date, final double mu)
195 throws IllegalArgumentException {
196 super(frame, date, mu);
197
198 if (a * (1 - e) < 0) {
199 throw new OrekitIllegalArgumentException(OrekitMessages.ORBIT_A_E_MISMATCH_WITH_CONIC_TYPE, a, e);
200 }
201
202 this.a = a;
203 this.aDot = aDot;
204 this.e = e;
205 this.eDot = eDot;
206 this.i = i;
207 this.iDot = iDot;
208 this.pa = pa;
209 this.paDot = paDot;
210 this.raan = raan;
211 this.raanDot = raanDot;
212
213 if (hasDerivatives()) {
214 final DerivativeStructure eDS = FACTORY.build(e, eDot);
215 final DerivativeStructure anomalyDS = FACTORY.build(anomaly, anomalyDot);
216 final DerivativeStructure vDS;
217 switch (type) {
218 case MEAN :
219 vDS = (a < 0) ?
220 FieldKeplerianOrbit.hyperbolicEccentricToTrue(FieldKeplerianOrbit.meanToHyperbolicEccentric(anomalyDS, eDS), eDS) :
221 FieldKeplerianOrbit.ellipticEccentricToTrue(FieldKeplerianOrbit.meanToEllipticEccentric(anomalyDS, eDS), eDS);
222 break;
223 case ECCENTRIC :
224 vDS = (a < 0) ?
225 FieldKeplerianOrbit.hyperbolicEccentricToTrue(anomalyDS, eDS) :
226 FieldKeplerianOrbit.ellipticEccentricToTrue(anomalyDS, eDS);
227 break;
228 case TRUE :
229 vDS = anomalyDS;
230 break;
231 default :
232 throw new OrekitInternalError(null);
233 }
234 this.v = vDS.getValue();
235 this.vDot = vDS.getPartialDerivative(1);
236 } else {
237 switch (type) {
238 case MEAN :
239 this.v = (a < 0) ?
240 hyperbolicEccentricToTrue(meanToHyperbolicEccentric(anomaly, e), e) :
241 ellipticEccentricToTrue(meanToEllipticEccentric(anomaly, e), e);
242 break;
243 case ECCENTRIC :
244 this.v = (a < 0) ?
245 hyperbolicEccentricToTrue(anomaly, e) :
246 ellipticEccentricToTrue(anomaly, e);
247 break;
248 case TRUE :
249 this.v = anomaly;
250 break;
251 default :
252 throw new OrekitInternalError(null);
253 }
254 this.vDot = Double.NaN;
255 }
256
257
258 if (1 + e * FastMath.cos(v) <= 0) {
259 final double vMax = FastMath.acos(-1 / e);
260 throw new OrekitIllegalArgumentException(OrekitMessages.ORBIT_ANOMALY_OUT_OF_HYPERBOLIC_RANGE,
261 v, e, -vMax, vMax);
262 }
263
264 this.partialPV = null;
265
266 }
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282 public KeplerianOrbit(final TimeStampedPVCoordinates pvCoordinates,
283 final Frame frame, final double mu)
284 throws IllegalArgumentException {
285 this(pvCoordinates, frame, mu, hasNonKeplerianAcceleration(pvCoordinates, mu));
286 }
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303 private KeplerianOrbit(final TimeStampedPVCoordinates pvCoordinates,
304 final Frame frame, final double mu,
305 final boolean reliableAcceleration)
306 throws IllegalArgumentException {
307 super(pvCoordinates, frame, mu);
308
309
310 final Vector3D momentum = pvCoordinates.getMomentum();
311 final double m2 = momentum.getNormSq();
312 i = Vector3D.angle(momentum, Vector3D.PLUS_K);
313
314
315 raan = Vector3D.crossProduct(Vector3D.PLUS_K, momentum).getAlpha();
316
317
318 final Vector3D pvP = pvCoordinates.getPosition();
319 final Vector3D pvV = pvCoordinates.getVelocity();
320 final Vector3D pvA = pvCoordinates.getAcceleration();
321 final double r2 = pvP.getNormSq();
322 final double r = FastMath.sqrt(r2);
323 final double V2 = pvV.getNormSq();
324 final double rV2OnMu = r * V2 / mu;
325
326
327 a = r / (2 - rV2OnMu);
328 final double muA = mu * a;
329
330
331 if (a > 0) {
332
333 final double eSE = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(muA);
334 final double eCE = rV2OnMu - 1;
335 e = FastMath.sqrt(eSE * eSE + eCE * eCE);
336 v = ellipticEccentricToTrue(FastMath.atan2(eSE, eCE), e);
337 } else {
338
339 final double eSH = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(-muA);
340 final double eCH = rV2OnMu - 1;
341 e = FastMath.sqrt(1 - m2 / muA);
342 v = hyperbolicEccentricToTrue(FastMath.log((eCH + eSH) / (eCH - eSH)) / 2, e);
343 }
344
345
346 final Vector3D node = new Vector3D(raan, 0.0);
347 final double px = Vector3D.dotProduct(pvP, node);
348 final double py = Vector3D.dotProduct(pvP, Vector3D.crossProduct(momentum, node)) / FastMath.sqrt(m2);
349 pa = FastMath.atan2(py, px) - v;
350
351 partialPV = pvCoordinates;
352
353 if (reliableAcceleration) {
354
355
356 final double[][] jacobian = new double[6][6];
357 getJacobianWrtCartesian(PositionAngle.MEAN, jacobian);
358
359 final Vector3D keplerianAcceleration = new Vector3D(-mu / (r * r2), pvP);
360 final Vector3D nonKeplerianAcceleration = pvA.subtract(keplerianAcceleration);
361 final double aX = nonKeplerianAcceleration.getX();
362 final double aY = nonKeplerianAcceleration.getY();
363 final double aZ = nonKeplerianAcceleration.getZ();
364 aDot = jacobian[0][3] * aX + jacobian[0][4] * aY + jacobian[0][5] * aZ;
365 eDot = jacobian[1][3] * aX + jacobian[1][4] * aY + jacobian[1][5] * aZ;
366 iDot = jacobian[2][3] * aX + jacobian[2][4] * aY + jacobian[2][5] * aZ;
367 paDot = jacobian[3][3] * aX + jacobian[3][4] * aY + jacobian[3][5] * aZ;
368 raanDot = jacobian[4][3] * aX + jacobian[4][4] * aY + jacobian[4][5] * aZ;
369
370
371
372 final double MDot = getKeplerianMeanMotion() +
373 jacobian[5][3] * aX + jacobian[5][4] * aY + jacobian[5][5] * aZ;
374 final DerivativeStructure eDS = FACTORY.build(e, eDot);
375 final DerivativeStructure MDS = FACTORY.build(getMeanAnomaly(), MDot);
376 final DerivativeStructure vDS = (a < 0) ?
377 FieldKeplerianOrbit.hyperbolicEccentricToTrue(FieldKeplerianOrbit.meanToHyperbolicEccentric(MDS, eDS), eDS) :
378 FieldKeplerianOrbit.ellipticEccentricToTrue(FieldKeplerianOrbit.meanToEllipticEccentric(MDS, eDS), eDS);
379 vDot = vDS.getPartialDerivative(1);
380
381 } else {
382
383
384
385 aDot = Double.NaN;
386 eDot = Double.NaN;
387 iDot = Double.NaN;
388 paDot = Double.NaN;
389 raanDot = Double.NaN;
390 vDot = Double.NaN;
391 }
392
393 }
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410 public KeplerianOrbit(final PVCoordinates pvCoordinates,
411 final Frame frame, final AbsoluteDate date, final double mu)
412 throws IllegalArgumentException {
413 this(new TimeStampedPVCoordinates(date, pvCoordinates), frame, mu);
414 }
415
416
417
418
419 public KeplerianOrbit(final Orbit op) {
420 this(op.getPVCoordinates(), op.getFrame(), op.getMu(), op.hasDerivatives());
421 }
422
423
424 public OrbitType getType() {
425 return OrbitType.KEPLERIAN;
426 }
427
428
429 public double getA() {
430 return a;
431 }
432
433
434 public double getADot() {
435 return aDot;
436 }
437
438
439 public double getE() {
440 return e;
441 }
442
443
444 public double getEDot() {
445 return eDot;
446 }
447
448
449 public double getI() {
450 return i;
451 }
452
453
454 public double getIDot() {
455 return iDot;
456 }
457
458
459
460
461 public double getPerigeeArgument() {
462 return pa;
463 }
464
465
466
467
468
469
470
471
472 public double getPerigeeArgumentDot() {
473 return paDot;
474 }
475
476
477
478
479 public double getRightAscensionOfAscendingNode() {
480 return raan;
481 }
482
483
484
485
486
487
488
489
490 public double getRightAscensionOfAscendingNodeDot() {
491 return raanDot;
492 }
493
494
495
496
497 public double getTrueAnomaly() {
498 return v;
499 }
500
501
502
503
504 public double getTrueAnomalyDot() {
505 return vDot;
506 }
507
508
509
510
511 public double getEccentricAnomaly() {
512 return (a < 0) ? trueToHyperbolicEccentric(v, e) : trueToEllipticEccentric(v, e);
513 }
514
515
516
517
518
519 public double getEccentricAnomalyDot() {
520 final DerivativeStructure eDS = FACTORY.build(e, eDot);
521 final DerivativeStructure vDS = FACTORY.build(v, vDot);
522 final DerivativeStructure EDS = (a < 0) ?
523 FieldKeplerianOrbit.trueToHyperbolicEccentric(vDS, eDS) :
524 FieldKeplerianOrbit.trueToEllipticEccentric(vDS, eDS);
525 return EDS.getPartialDerivative(1);
526 }
527
528
529
530
531 public double getMeanAnomaly() {
532 return (a < 0) ?
533 hyperbolicEccentricToMean(trueToHyperbolicEccentric(v, e), e) :
534 ellipticEccentricToMean(trueToEllipticEccentric(v, e), e);
535 }
536
537
538
539
540
541 public double getMeanAnomalyDot() {
542 final DerivativeStructure eDS = FACTORY.build(e, eDot);
543 final DerivativeStructure vDS = FACTORY.build(v, vDot);
544 final DerivativeStructure MDS = (a < 0) ?
545 FieldKeplerianOrbit.hyperbolicEccentricToMean(FieldKeplerianOrbit.trueToHyperbolicEccentric(vDS, eDS), eDS) :
546 FieldKeplerianOrbit.ellipticEccentricToMean(FieldKeplerianOrbit.trueToEllipticEccentric(vDS, eDS), eDS);
547 return MDS.getPartialDerivative(1);
548 }
549
550
551
552
553
554 public double getAnomaly(final PositionAngle type) {
555 return (type == PositionAngle.MEAN) ? getMeanAnomaly() :
556 ((type == PositionAngle.ECCENTRIC) ? getEccentricAnomaly() :
557 getTrueAnomaly());
558 }
559
560
561
562
563
564
565 public double getAnomalyDot(final PositionAngle type) {
566 return (type == PositionAngle.MEAN) ? getMeanAnomalyDot() :
567 ((type == PositionAngle.ECCENTRIC) ? getEccentricAnomalyDot() :
568 getTrueAnomalyDot());
569 }
570
571
572
573
574
575
576
577 public static double ellipticEccentricToTrue(final double E, final double e) {
578 final double beta = e / (1 + FastMath.sqrt((1 - e) * (1 + e)));
579 return E + 2 * FastMath.atan(beta * FastMath.sin(E) / (1 - beta * FastMath.cos(E)));
580 }
581
582
583
584
585
586
587
588 public static double trueToEllipticEccentric(final double v, final double e) {
589 final double beta = e / (1 + FastMath.sqrt(1 - e * e));
590 return v - 2 * FastMath.atan(beta * FastMath.sin(v) / (1 + beta * FastMath.cos(v)));
591 }
592
593
594
595
596
597
598 public static double hyperbolicEccentricToTrue(final double H, final double e) {
599 return 2 * FastMath.atan(FastMath.sqrt((e + 1) / (e - 1)) * FastMath.tanh(H / 2));
600 }
601
602
603
604
605
606
607
608 public static double trueToHyperbolicEccentric(final double v, final double e) {
609 final double sinhH = FastMath.sqrt(e * e - 1) * FastMath.sin(v) / (1 + e * FastMath.cos(v));
610 return FastMath.asinh(sinhH);
611 }
612
613
614
615
616
617
618
619
620
621
622
623 public static double meanToEllipticEccentric(final double M, final double e) {
624
625
626 final double reducedM = MathUtils.normalizeAngle(M, 0.0);
627
628
629 double E;
630 if (FastMath.abs(reducedM) < 1.0 / 6.0) {
631 E = reducedM + e * (FastMath.cbrt(6 * reducedM) - reducedM);
632 } else {
633 if (reducedM < 0) {
634 final double w = FastMath.PI + reducedM;
635 E = reducedM + e * (A * w / (B - w) - FastMath.PI - reducedM);
636 } else {
637 final double w = FastMath.PI - reducedM;
638 E = reducedM + e * (FastMath.PI - A * w / (B - w) - reducedM);
639 }
640 }
641
642 final double e1 = 1 - e;
643 final boolean noCancellationRisk = (e1 + E * E / 6) >= 0.1;
644
645
646 for (int j = 0; j < 2; ++j) {
647 final double f;
648 double fd;
649 final double fdd = e * FastMath.sin(E);
650 final double fddd = e * FastMath.cos(E);
651 if (noCancellationRisk) {
652 f = (E - fdd) - reducedM;
653 fd = 1 - fddd;
654 } else {
655 f = eMeSinE(E, e) - reducedM;
656 final double s = FastMath.sin(0.5 * E);
657 fd = e1 + 2 * e * s * s;
658 }
659 final double dee = f * fd / (0.5 * f * fdd - fd * fd);
660
661
662 final double w = fd + 0.5 * dee * (fdd + dee * fddd / 3);
663 fd += dee * (fdd + 0.5 * dee * fddd);
664 E -= (f - dee * (fd - w)) / fd;
665
666 }
667
668
669 E += M - reducedM;
670
671 return E;
672
673 }
674
675
676
677
678
679
680
681
682
683
684 private static double eMeSinE(final double E, final double e) {
685 double x = (1 - e) * FastMath.sin(E);
686 final double mE2 = -E * E;
687 double term = E;
688 double d = 0;
689
690 for (double x0 = Double.NaN; x != x0;) {
691 d += 2;
692 term *= mE2 / (d * (d + 1));
693 x0 = x;
694 x = x - term;
695 }
696 return x;
697 }
698
699
700
701
702
703
704
705
706
707
708 public static double meanToHyperbolicEccentric(final double M, final double ecc) {
709
710
711
712
713 double H;
714 if (ecc < 1.6) {
715 if ((-FastMath.PI < M && M < 0.) || M > FastMath.PI) {
716 H = M - ecc;
717 } else {
718 H = M + ecc;
719 }
720 } else {
721 if (ecc < 3.6 && FastMath.abs(M) > FastMath.PI) {
722 H = M - FastMath.copySign(ecc, M);
723 } else {
724 H = M / (ecc - 1.);
725 }
726 }
727
728
729 int iter = 0;
730 do {
731 final double f3 = ecc * FastMath.cosh(H);
732 final double f2 = ecc * FastMath.sinh(H);
733 final double f1 = f3 - 1.;
734 final double f0 = f2 - H - M;
735 final double f12 = 2. * f1;
736 final double d = f0 / f12;
737 final double fdf = f1 - d * f2;
738 final double ds = f0 / fdf;
739
740 final double shift = f0 / (fdf + ds * ds * f3 / 6.);
741
742 H -= shift;
743
744 if (FastMath.abs(shift) <= 1.0e-12) {
745 return H;
746 }
747
748 } while (++iter < 50);
749
750 throw new MathIllegalStateException(OrekitMessages.UNABLE_TO_COMPUTE_HYPERBOLIC_ECCENTRIC_ANOMALY,
751 iter);
752 }
753
754
755
756
757
758
759
760 public static double ellipticEccentricToMean(final double E, final double e) {
761 return E - e * FastMath.sin(E);
762 }
763
764
765
766
767
768
769
770 public static double hyperbolicEccentricToMean(final double H, final double e) {
771 return e * FastMath.sinh(H) - H;
772 }
773
774
775 public double getEquinoctialEx() {
776 return e * FastMath.cos(pa + raan);
777 }
778
779
780 public double getEquinoctialExDot() {
781 final double paPraan = pa + raan;
782 final double cos = FastMath.cos(paPraan);
783 final double sin = FastMath.sin(paPraan);
784 return eDot * cos - e * sin * (paDot + raanDot);
785 }
786
787
788 public double getEquinoctialEy() {
789 return e * FastMath.sin(pa + raan);
790 }
791
792
793 public double getEquinoctialEyDot() {
794 final double paPraan = pa + raan;
795 final double cos = FastMath.cos(paPraan);
796 final double sin = FastMath.sin(paPraan);
797 return eDot * sin + e * cos * (paDot + raanDot);
798 }
799
800
801 public double getHx() {
802
803 if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
804 return Double.NaN;
805 }
806 return FastMath.cos(raan) * FastMath.tan(0.5 * i);
807 }
808
809
810 public double getHxDot() {
811
812 if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
813 return Double.NaN;
814 }
815 final double cosRaan = FastMath.cos(raan);
816 final double sinRaan = FastMath.sin(raan);
817 final double tan = FastMath.tan(0.5 * i);
818 return 0.5 * (1 + tan * tan) * cosRaan * iDot - tan * sinRaan * raanDot;
819 }
820
821
822 public double getHy() {
823
824 if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
825 return Double.NaN;
826 }
827 return FastMath.sin(raan) * FastMath.tan(0.5 * i);
828 }
829
830
831 public double getHyDot() {
832
833 if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
834 return Double.NaN;
835 }
836 final double cosRaan = FastMath.cos(raan);
837 final double sinRaan = FastMath.sin(raan);
838 final double tan = FastMath.tan(0.5 * i);
839 return 0.5 * (1 + tan * tan) * sinRaan * iDot + tan * cosRaan * raanDot;
840 }
841
842
843 public double getLv() {
844 return pa + raan + v;
845 }
846
847
848 public double getLvDot() {
849 return paDot + raanDot + vDot;
850 }
851
852
853 public double getLE() {
854 return pa + raan + getEccentricAnomaly();
855 }
856
857
858 public double getLEDot() {
859 return paDot + raanDot + getEccentricAnomalyDot();
860 }
861
862
863 public double getLM() {
864 return pa + raan + getMeanAnomaly();
865 }
866
867
868 public double getLMDot() {
869 return paDot + raanDot + getMeanAnomalyDot();
870 }
871
872
873
874 private void computePVWithoutA() {
875
876 if (partialPV != null) {
877
878 return;
879 }
880
881
882 final double cosRaan = FastMath.cos(raan);
883 final double sinRaan = FastMath.sin(raan);
884 final double cosPa = FastMath.cos(pa);
885 final double sinPa = FastMath.sin(pa);
886 final double cosI = FastMath.cos(i);
887 final double sinI = FastMath.sin(i);
888
889 final double crcp = cosRaan * cosPa;
890 final double crsp = cosRaan * sinPa;
891 final double srcp = sinRaan * cosPa;
892 final double srsp = sinRaan * sinPa;
893
894
895 final Vector3D p = new Vector3D( crcp - cosI * srsp, srcp + cosI * crsp, sinI * sinPa);
896 final Vector3D q = new Vector3D(-crsp - cosI * srcp, -srsp + cosI * crcp, sinI * cosPa);
897
898 if (a > 0) {
899
900
901
902
903 final double uME2 = (1 - e) * (1 + e);
904 final double s1Me2 = FastMath.sqrt(uME2);
905 final double E = getEccentricAnomaly();
906 final double cosE = FastMath.cos(E);
907 final double sinE = FastMath.sin(E);
908
909
910 final double x = a * (cosE - e);
911 final double y = a * sinE * s1Me2;
912 final double factor = FastMath.sqrt(getMu() / a) / (1 - e * cosE);
913 final double xDot = -sinE * factor;
914 final double yDot = cosE * s1Me2 * factor;
915
916 final Vector3D position = new Vector3D(x, p, y, q);
917 final Vector3D velocity = new Vector3D(xDot, p, yDot, q);
918 partialPV = new PVCoordinates(position, velocity);
919
920 } else {
921
922
923
924
925 final double sinV = FastMath.sin(v);
926 final double cosV = FastMath.cos(v);
927 final double f = a * (1 - e * e);
928 final double posFactor = f / (1 + e * cosV);
929 final double velFactor = FastMath.sqrt(getMu() / f);
930
931 final double x = posFactor * cosV;
932 final double y = posFactor * sinV;
933 final double xDot = -velFactor * sinV;
934 final double yDot = velFactor * (e + cosV);
935
936 final Vector3D position = new Vector3D(x, p, y, q);
937 final Vector3D velocity = new Vector3D(xDot, p, yDot, q);
938 partialPV = new PVCoordinates(position, velocity);
939
940 }
941
942 }
943
944
945
946
947
948
949
950 private Vector3D nonKeplerianAcceleration() {
951
952 final double[][] dCdP = new double[6][6];
953 getJacobianWrtParameters(PositionAngle.MEAN, dCdP);
954
955 final double nonKeplerianMeanMotion = getMeanAnomalyDot() - getKeplerianMeanMotion();
956 final double nonKeplerianAx = dCdP[3][0] * aDot + dCdP[3][1] * eDot + dCdP[3][2] * iDot +
957 dCdP[3][3] * paDot + dCdP[3][4] * raanDot + dCdP[3][5] * nonKeplerianMeanMotion;
958 final double nonKeplerianAy = dCdP[4][0] * aDot + dCdP[4][1] * eDot + dCdP[4][2] * iDot +
959 dCdP[4][3] * paDot + dCdP[4][4] * raanDot + dCdP[4][5] * nonKeplerianMeanMotion;
960 final double nonKeplerianAz = dCdP[5][0] * aDot + dCdP[5][1] * eDot + dCdP[5][2] * iDot +
961 dCdP[5][3] * paDot + dCdP[5][4] * raanDot + dCdP[5][5] * nonKeplerianMeanMotion;
962
963 return new Vector3D(nonKeplerianAx, nonKeplerianAy, nonKeplerianAz);
964
965 }
966
967
968 protected TimeStampedPVCoordinates initPVCoordinates() {
969
970
971 computePVWithoutA();
972
973
974 final double r2 = partialPV.getPosition().getNormSq();
975 final Vector3D keplerianAcceleration = new Vector3D(-getMu() / (r2 * FastMath.sqrt(r2)), partialPV.getPosition());
976 final Vector3D acceleration = hasDerivatives() ?
977 keplerianAcceleration.add(nonKeplerianAcceleration()) :
978 keplerianAcceleration;
979
980 return new TimeStampedPVCoordinates(getDate(), partialPV.getPosition(), partialPV.getVelocity(), acceleration);
981
982 }
983
984
985 public KeplerianOrbit shiftedBy(final double dt) {
986
987
988 final KeplerianOrbitbit">KeplerianOrbit keplerianShifted = new KeplerianOrbit(a, e, i, pa, raan,
989 getMeanAnomaly() + getKeplerianMeanMotion() * dt,
990 PositionAngle.MEAN, getFrame(),
991 getDate().shiftedBy(dt), getMu());
992
993 if (hasDerivatives()) {
994
995
996 final Vector3D nonKeplerianAcceleration = nonKeplerianAcceleration();
997
998
999 keplerianShifted.computePVWithoutA();
1000 final Vector3D fixedP = new Vector3D(1, keplerianShifted.partialPV.getPosition(),
1001 0.5 * dt * dt, nonKeplerianAcceleration);
1002 final double fixedR2 = fixedP.getNormSq();
1003 final double fixedR = FastMath.sqrt(fixedR2);
1004 final Vector3D fixedV = new Vector3D(1, keplerianShifted.partialPV.getVelocity(),
1005 dt, nonKeplerianAcceleration);
1006 final Vector3D fixedA = new Vector3D(-getMu() / (fixedR2 * fixedR), keplerianShifted.partialPV.getPosition(),
1007 1, nonKeplerianAcceleration);
1008
1009
1010 return new KeplerianOrbit(new TimeStampedPVCoordinates(keplerianShifted.getDate(),
1011 fixedP, fixedV, fixedA),
1012 keplerianShifted.getFrame(), keplerianShifted.getMu());
1013
1014 } else {
1015
1016 return keplerianShifted;
1017 }
1018
1019 }
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041 public KeplerianOrbit interpolate(final AbsoluteDate date, final Stream<Orbit> sample) {
1042
1043
1044 final List<Orbit> list = sample.collect(Collectors.toList());
1045 boolean useDerivatives = true;
1046 for (final Orbit orbit : list) {
1047 useDerivatives = useDerivatives && orbit.hasDerivatives();
1048 }
1049
1050
1051 final HermiteInterpolator interpolator = new HermiteInterpolator();
1052
1053
1054 AbsoluteDate previousDate = null;
1055 double previousPA = Double.NaN;
1056 double previousRAAN = Double.NaN;
1057 double previousM = Double.NaN;
1058 for (final Orbit orbit : list) {
1059 final KeplerianOrbit/../org/orekit/orbits/KeplerianOrbit.html#KeplerianOrbit">KeplerianOrbit kep = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(orbit);
1060 final double continuousPA;
1061 final double continuousRAAN;
1062 final double continuousM;
1063 if (previousDate == null) {
1064 continuousPA = kep.getPerigeeArgument();
1065 continuousRAAN = kep.getRightAscensionOfAscendingNode();
1066 continuousM = kep.getMeanAnomaly();
1067 } else {
1068 final double dt = kep.getDate().durationFrom(previousDate);
1069 final double keplerM = previousM + kep.getKeplerianMeanMotion() * dt;
1070 continuousPA = MathUtils.normalizeAngle(kep.getPerigeeArgument(), previousPA);
1071 continuousRAAN = MathUtils.normalizeAngle(kep.getRightAscensionOfAscendingNode(), previousRAAN);
1072 continuousM = MathUtils.normalizeAngle(kep.getMeanAnomaly(), keplerM);
1073 }
1074 previousDate = kep.getDate();
1075 previousPA = continuousPA;
1076 previousRAAN = continuousRAAN;
1077 previousM = continuousM;
1078 if (useDerivatives) {
1079 interpolator.addSamplePoint(kep.getDate().durationFrom(date),
1080 new double[] {
1081 kep.getA(),
1082 kep.getE(),
1083 kep.getI(),
1084 continuousPA,
1085 continuousRAAN,
1086 continuousM
1087 }, new double[] {
1088 kep.getADot(),
1089 kep.getEDot(),
1090 kep.getIDot(),
1091 kep.getPerigeeArgumentDot(),
1092 kep.getRightAscensionOfAscendingNodeDot(),
1093 kep.getMeanAnomalyDot()
1094 });
1095 } else {
1096 interpolator.addSamplePoint(kep.getDate().durationFrom(date),
1097 new double[] {
1098 kep.getA(),
1099 kep.getE(),
1100 kep.getI(),
1101 continuousPA,
1102 continuousRAAN,
1103 continuousM
1104 });
1105 }
1106 }
1107
1108
1109 final double[][] interpolated = interpolator.derivatives(0.0, 1);
1110
1111
1112 return new KeplerianOrbit(interpolated[0][0], interpolated[0][1], interpolated[0][2],
1113 interpolated[0][3], interpolated[0][4], interpolated[0][5],
1114 interpolated[1][0], interpolated[1][1], interpolated[1][2],
1115 interpolated[1][3], interpolated[1][4], interpolated[1][5],
1116 PositionAngle.MEAN, getFrame(), date, getMu());
1117
1118 }
1119
1120
1121 protected double[][] computeJacobianMeanWrtCartesian() {
1122 if (a > 0) {
1123 return computeJacobianMeanWrtCartesianElliptical();
1124 } else {
1125 return computeJacobianMeanWrtCartesianHyperbolic();
1126 }
1127 }
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137 private double[][] computeJacobianMeanWrtCartesianElliptical() {
1138
1139 final double[][] jacobian = new double[6][6];
1140
1141
1142 computePVWithoutA();
1143 final Vector3D position = partialPV.getPosition();
1144 final Vector3D velocity = partialPV.getVelocity();
1145 final Vector3D momentum = partialPV.getMomentum();
1146 final double v2 = velocity.getNormSq();
1147 final double r2 = position.getNormSq();
1148 final double r = FastMath.sqrt(r2);
1149 final double r3 = r * r2;
1150
1151 final double px = position.getX();
1152 final double py = position.getY();
1153 final double pz = position.getZ();
1154 final double vx = velocity.getX();
1155 final double vy = velocity.getY();
1156 final double vz = velocity.getZ();
1157 final double mx = momentum.getX();
1158 final double my = momentum.getY();
1159 final double mz = momentum.getZ();
1160
1161 final double mu = getMu();
1162 final double sqrtMuA = FastMath.sqrt(a * mu);
1163 final double sqrtAoMu = FastMath.sqrt(a / mu);
1164 final double a2 = a * a;
1165 final double twoA = 2 * a;
1166 final double rOnA = r / a;
1167
1168 final double oMe2 = 1 - e * e;
1169 final double epsilon = FastMath.sqrt(oMe2);
1170 final double sqrtRec = 1 / epsilon;
1171
1172 final double cosI = FastMath.cos(i);
1173 final double sinI = FastMath.sin(i);
1174 final double cosPA = FastMath.cos(pa);
1175 final double sinPA = FastMath.sin(pa);
1176
1177 final double pv = Vector3D.dotProduct(position, velocity);
1178 final double cosE = (a - r) / (a * e);
1179 final double sinE = pv / (e * sqrtMuA);
1180
1181
1182 final Vector3D vectorAR = new Vector3D(2 * a2 / r3, position);
1183 final Vector3D vectorARDot = velocity.scalarMultiply(2 * a2 / mu);
1184 fillHalfRow(1, vectorAR, jacobian[0], 0);
1185 fillHalfRow(1, vectorARDot, jacobian[0], 3);
1186
1187
1188 final double factorER3 = pv / twoA;
1189 final Vector3D vectorER = new Vector3D(cosE * v2 / (r * mu), position,
1190 sinE / sqrtMuA, velocity,
1191 -factorER3 * sinE / sqrtMuA, vectorAR);
1192 final Vector3D vectorERDot = new Vector3D(sinE / sqrtMuA, position,
1193 cosE * 2 * r / mu, velocity,
1194 -factorER3 * sinE / sqrtMuA, vectorARDot);
1195 fillHalfRow(1, vectorER, jacobian[1], 0);
1196 fillHalfRow(1, vectorERDot, jacobian[1], 3);
1197
1198
1199 final double coefE = cosE / (e * sqrtMuA);
1200 final Vector3D vectorEAnR =
1201 new Vector3D(-sinE * v2 / (e * r * mu), position, coefE, velocity,
1202 -factorER3 * coefE, vectorAR);
1203
1204
1205 final Vector3D vectorEAnRDot =
1206 new Vector3D(-sinE * 2 * r / (e * mu), velocity, coefE, position,
1207 -factorER3 * coefE, vectorARDot);
1208
1209
1210 final double s1 = -sinE * pz / r - cosE * vz * sqrtAoMu;
1211 final double s2 = -cosE * pz / r3;
1212 final double s3 = -sinE * vz / (2 * sqrtMuA);
1213 final double t1 = sqrtRec * (cosE * pz / r - sinE * vz * sqrtAoMu);
1214 final double t2 = sqrtRec * (-sinE * pz / r3);
1215 final double t3 = sqrtRec * (cosE - e) * vz / (2 * sqrtMuA);
1216 final double t4 = sqrtRec * (e * sinI * cosPA * sqrtRec - vz * sqrtAoMu);
1217 final Vector3D s = new Vector3D(cosE / r, Vector3D.PLUS_K,
1218 s1, vectorEAnR,
1219 s2, position,
1220 s3, vectorAR);
1221 final Vector3D sDot = new Vector3D(-sinE * sqrtAoMu, Vector3D.PLUS_K,
1222 s1, vectorEAnRDot,
1223 s3, vectorARDot);
1224 final Vector3D t =
1225 new Vector3D(sqrtRec * sinE / r, Vector3D.PLUS_K).add(new Vector3D(t1, vectorEAnR,
1226 t2, position,
1227 t3, vectorAR,
1228 t4, vectorER));
1229 final Vector3D tDot = new Vector3D(sqrtRec * (cosE - e) * sqrtAoMu, Vector3D.PLUS_K,
1230 t1, vectorEAnRDot,
1231 t3, vectorARDot,
1232 t4, vectorERDot);
1233
1234
1235 final double factorI1 = -sinI * sqrtRec / sqrtMuA;
1236 final double i1 = factorI1;
1237 final double i2 = -factorI1 * mz / twoA;
1238 final double i3 = factorI1 * mz * e / oMe2;
1239 final double i4 = cosI * sinPA;
1240 final double i5 = cosI * cosPA;
1241 fillHalfRow(i1, new Vector3D(vy, -vx, 0), i2, vectorAR, i3, vectorER, i4, s, i5, t,
1242 jacobian[2], 0);
1243 fillHalfRow(i1, new Vector3D(-py, px, 0), i2, vectorARDot, i3, vectorERDot, i4, sDot, i5, tDot,
1244 jacobian[2], 3);
1245
1246
1247 fillHalfRow(cosPA / sinI, s, -sinPA / sinI, t, jacobian[3], 0);
1248 fillHalfRow(cosPA / sinI, sDot, -sinPA / sinI, tDot, jacobian[3], 3);
1249
1250
1251 final double factorRaanR = 1 / (mu * a * oMe2 * sinI * sinI);
1252 fillHalfRow(-factorRaanR * my, new Vector3D( 0, vz, -vy),
1253 factorRaanR * mx, new Vector3D(-vz, 0, vx),
1254 jacobian[4], 0);
1255 fillHalfRow(-factorRaanR * my, new Vector3D( 0, -pz, py),
1256 factorRaanR * mx, new Vector3D(pz, 0, -px),
1257 jacobian[4], 3);
1258
1259
1260 fillHalfRow(rOnA, vectorEAnR, -sinE, vectorER, jacobian[5], 0);
1261 fillHalfRow(rOnA, vectorEAnRDot, -sinE, vectorERDot, jacobian[5], 3);
1262
1263 return jacobian;
1264
1265 }
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275 private double[][] computeJacobianMeanWrtCartesianHyperbolic() {
1276
1277 final double[][] jacobian = new double[6][6];
1278
1279
1280 computePVWithoutA();
1281 final Vector3D position = partialPV.getPosition();
1282 final Vector3D velocity = partialPV.getVelocity();
1283 final Vector3D momentum = partialPV.getMomentum();
1284 final double r2 = position.getNormSq();
1285 final double r = FastMath.sqrt(r2);
1286 final double r3 = r * r2;
1287
1288 final double x = position.getX();
1289 final double y = position.getY();
1290 final double z = position.getZ();
1291 final double vx = velocity.getX();
1292 final double vy = velocity.getY();
1293 final double vz = velocity.getZ();
1294 final double mx = momentum.getX();
1295 final double my = momentum.getY();
1296 final double mz = momentum.getZ();
1297
1298 final double mu = getMu();
1299 final double absA = -a;
1300 final double sqrtMuA = FastMath.sqrt(absA * mu);
1301 final double a2 = a * a;
1302 final double rOa = r / absA;
1303
1304 final double cosI = FastMath.cos(i);
1305 final double sinI = FastMath.sin(i);
1306
1307 final double pv = Vector3D.dotProduct(position, velocity);
1308
1309
1310 final Vector3D vectorAR = new Vector3D(-2 * a2 / r3, position);
1311 final Vector3D vectorARDot = velocity.scalarMultiply(-2 * a2 / mu);
1312 fillHalfRow(-1, vectorAR, jacobian[0], 0);
1313 fillHalfRow(-1, vectorARDot, jacobian[0], 3);
1314
1315
1316 final double m = momentum.getNorm();
1317 final double oOm = 1 / m;
1318 final Vector3D dcXP = new Vector3D( 0, vz, -vy);
1319 final Vector3D dcYP = new Vector3D(-vz, 0, vx);
1320 final Vector3D dcZP = new Vector3D( vy, -vx, 0);
1321 final Vector3D dcXV = new Vector3D( 0, -z, y);
1322 final Vector3D dcYV = new Vector3D( z, 0, -x);
1323 final Vector3D dcZV = new Vector3D( -y, x, 0);
1324 final Vector3D dCP = new Vector3D(mx * oOm, dcXP, my * oOm, dcYP, mz * oOm, dcZP);
1325 final Vector3D dCV = new Vector3D(mx * oOm, dcXV, my * oOm, dcYV, mz * oOm, dcZV);
1326
1327
1328 final double mOMu = m / mu;
1329 final Vector3D dpP = new Vector3D(2 * mOMu, dCP);
1330 final Vector3D dpV = new Vector3D(2 * mOMu, dCV);
1331
1332
1333 final double p = m * mOMu;
1334 final double moO2ae = 1 / (2 * absA * e);
1335 final double m2OaMu = -p / absA;
1336 fillHalfRow(moO2ae, dpP, m2OaMu * moO2ae, vectorAR, jacobian[1], 0);
1337 fillHalfRow(moO2ae, dpV, m2OaMu * moO2ae, vectorARDot, jacobian[1], 3);
1338
1339
1340 final double cI1 = 1 / (m * sinI);
1341 final double cI2 = cosI * cI1;
1342 fillHalfRow(cI2, dCP, -cI1, dcZP, jacobian[2], 0);
1343 fillHalfRow(cI2, dCV, -cI1, dcZV, jacobian[2], 3);
1344
1345
1346 final double cP1 = y * oOm;
1347 final double cP2 = -x * oOm;
1348 final double cP3 = -(mx * cP1 + my * cP2);
1349 final double cP4 = cP3 * oOm;
1350 final double cP5 = -1 / (r2 * sinI * sinI);
1351 final double cP6 = z * cP5;
1352 final double cP7 = cP3 * cP5;
1353 final Vector3D dacP = new Vector3D(cP1, dcXP, cP2, dcYP, cP4, dCP, oOm, new Vector3D(-my, mx, 0));
1354 final Vector3D dacV = new Vector3D(cP1, dcXV, cP2, dcYV, cP4, dCV);
1355 final Vector3D dpoP = new Vector3D(cP6, dacP, cP7, Vector3D.PLUS_K);
1356 final Vector3D dpoV = new Vector3D(cP6, dacV);
1357
1358 final double re2 = r2 * e * e;
1359 final double recOre2 = (p - r) / re2;
1360 final double resOre2 = (pv * mOMu) / re2;
1361 final Vector3D dreP = new Vector3D(mOMu, velocity, pv / mu, dCP);
1362 final Vector3D dreV = new Vector3D(mOMu, position, pv / mu, dCV);
1363 final Vector3D davP = new Vector3D(-resOre2, dpP, recOre2, dreP, resOre2 / r, position);
1364 final Vector3D davV = new Vector3D(-resOre2, dpV, recOre2, dreV);
1365 fillHalfRow(1, dpoP, -1, davP, jacobian[3], 0);
1366 fillHalfRow(1, dpoV, -1, davV, jacobian[3], 3);
1367
1368
1369 final double cO0 = cI1 * cI1;
1370 final double cO1 = mx * cO0;
1371 final double cO2 = -my * cO0;
1372 fillHalfRow(cO1, dcYP, cO2, dcXP, jacobian[4], 0);
1373 fillHalfRow(cO1, dcYV, cO2, dcXV, jacobian[4], 3);
1374
1375
1376 final double s2a = pv / (2 * absA);
1377 final double oObux = 1 / FastMath.sqrt(m * m + mu * absA);
1378 final double scasbu = pv * oObux;
1379 final Vector3D dauP = new Vector3D(1 / sqrtMuA, velocity, -s2a / sqrtMuA, vectorAR);
1380 final Vector3D dauV = new Vector3D(1 / sqrtMuA, position, -s2a / sqrtMuA, vectorARDot);
1381 final Vector3D dbuP = new Vector3D(oObux * mu / 2, vectorAR, m * oObux, dCP);
1382 final Vector3D dbuV = new Vector3D(oObux * mu / 2, vectorARDot, m * oObux, dCV);
1383 final Vector3D dcuP = new Vector3D(oObux, velocity, -scasbu * oObux, dbuP);
1384 final Vector3D dcuV = new Vector3D(oObux, position, -scasbu * oObux, dbuV);
1385 fillHalfRow(1, dauP, -e / (1 + rOa), dcuP, jacobian[5], 0);
1386 fillHalfRow(1, dauV, -e / (1 + rOa), dcuV, jacobian[5], 3);
1387
1388 return jacobian;
1389
1390 }
1391
1392
1393 protected double[][] computeJacobianEccentricWrtCartesian() {
1394 if (a > 0) {
1395 return computeJacobianEccentricWrtCartesianElliptical();
1396 } else {
1397 return computeJacobianEccentricWrtCartesianHyperbolic();
1398 }
1399 }
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409 private double[][] computeJacobianEccentricWrtCartesianElliptical() {
1410
1411
1412 final double[][] jacobian = computeJacobianMeanWrtCartesianElliptical();
1413
1414
1415
1416
1417
1418 final double eccentricAnomaly = getEccentricAnomaly();
1419 final double cosE = FastMath.cos(eccentricAnomaly);
1420 final double sinE = FastMath.sin(eccentricAnomaly);
1421 final double aOr = 1 / (1 - e * cosE);
1422
1423
1424 final double[] eRow = jacobian[1];
1425 final double[] anomalyRow = jacobian[5];
1426 for (int j = 0; j < anomalyRow.length; ++j) {
1427 anomalyRow[j] = aOr * (anomalyRow[j] + sinE * eRow[j]);
1428 }
1429
1430 return jacobian;
1431
1432 }
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442 private double[][] computeJacobianEccentricWrtCartesianHyperbolic() {
1443
1444
1445 final double[][] jacobian = computeJacobianMeanWrtCartesianHyperbolic();
1446
1447
1448
1449
1450
1451 final double H = getEccentricAnomaly();
1452 final double coshH = FastMath.cosh(H);
1453 final double sinhH = FastMath.sinh(H);
1454 final double absaOr = 1 / (e * coshH - 1);
1455
1456
1457 final double[] eRow = jacobian[1];
1458 final double[] anomalyRow = jacobian[5];
1459 for (int j = 0; j < anomalyRow.length; ++j) {
1460 anomalyRow[j] = absaOr * (anomalyRow[j] - sinhH * eRow[j]);
1461 }
1462
1463 return jacobian;
1464
1465 }
1466
1467
1468 protected double[][] computeJacobianTrueWrtCartesian() {
1469 if (a > 0) {
1470 return computeJacobianTrueWrtCartesianElliptical();
1471 } else {
1472 return computeJacobianTrueWrtCartesianHyperbolic();
1473 }
1474 }
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484 private double[][] computeJacobianTrueWrtCartesianElliptical() {
1485
1486
1487 final double[][] jacobian = computeJacobianEccentricWrtCartesianElliptical();
1488
1489
1490
1491
1492
1493
1494 final double e2 = e * e;
1495 final double oMe2 = 1 - e2;
1496 final double epsilon = FastMath.sqrt(oMe2);
1497 final double eccentricAnomaly = getEccentricAnomaly();
1498 final double cosE = FastMath.cos(eccentricAnomaly);
1499 final double sinE = FastMath.sin(eccentricAnomaly);
1500 final double aOr = 1 / (1 - e * cosE);
1501 final double aFactor = epsilon * aOr;
1502 final double eFactor = sinE * aOr / epsilon;
1503
1504
1505 final double[] eRow = jacobian[1];
1506 final double[] anomalyRow = jacobian[5];
1507 for (int j = 0; j < anomalyRow.length; ++j) {
1508 anomalyRow[j] = aFactor * anomalyRow[j] + eFactor * eRow[j];
1509 }
1510
1511 return jacobian;
1512
1513 }
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523 private double[][] computeJacobianTrueWrtCartesianHyperbolic() {
1524
1525
1526 final double[][] jacobian = computeJacobianEccentricWrtCartesianHyperbolic();
1527
1528
1529
1530
1531
1532
1533 final double e2 = e * e;
1534 final double e2Mo = e2 - 1;
1535 final double epsilon = FastMath.sqrt(e2Mo);
1536 final double H = getEccentricAnomaly();
1537 final double coshH = FastMath.cosh(H);
1538 final double sinhH = FastMath.sinh(H);
1539 final double aOr = 1 / (e * coshH - 1);
1540 final double aFactor = epsilon * aOr;
1541 final double eFactor = sinhH * aOr / epsilon;
1542
1543
1544 final double[] eRow = jacobian[1];
1545 final double[] anomalyRow = jacobian[5];
1546 for (int j = 0; j < anomalyRow.length; ++j) {
1547 anomalyRow[j] = aFactor * anomalyRow[j] - eFactor * eRow[j];
1548 }
1549
1550 return jacobian;
1551
1552 }
1553
1554
1555 public void addKeplerContribution(final PositionAngle type, final double gm,
1556 final double[] pDot) {
1557 final double oMe2;
1558 final double ksi;
1559 final double absA = FastMath.abs(a);
1560 final double n = FastMath.sqrt(gm / absA) / absA;
1561 switch (type) {
1562 case MEAN :
1563 pDot[5] += n;
1564 break;
1565 case ECCENTRIC :
1566 oMe2 = FastMath.abs(1 - e * e);
1567 ksi = 1 + e * FastMath.cos(v);
1568 pDot[5] += n * ksi / oMe2;
1569 break;
1570 case TRUE :
1571 oMe2 = FastMath.abs(1 - e * e);
1572 ksi = 1 + e * FastMath.cos(v);
1573 pDot[5] += n * ksi * ksi / (oMe2 * FastMath.sqrt(oMe2));
1574 break;
1575 default :
1576 throw new OrekitInternalError(null);
1577 }
1578 }
1579
1580
1581
1582
1583 public String toString() {
1584 return new StringBuffer().append("Keplerian parameters: ").append('{').
1585 append("a: ").append(a).
1586 append("; e: ").append(e).
1587 append("; i: ").append(FastMath.toDegrees(i)).
1588 append("; pa: ").append(FastMath.toDegrees(pa)).
1589 append("; raan: ").append(FastMath.toDegrees(raan)).
1590 append("; v: ").append(FastMath.toDegrees(v)).
1591 append(";}").toString();
1592 }
1593
1594
1595
1596
1597 private Object writeReplace() {
1598 return new DTO(this);
1599 }
1600
1601
1602 private static class DTO implements Serializable {
1603
1604
1605 private static final long serialVersionUID = 20170414L;
1606
1607
1608 private double[] d;
1609
1610
1611 private final Frame frame;
1612
1613
1614
1615
1616 private DTO(final KeplerianOrbit orbit) {
1617
1618 final TimeStampedPVCoordinates pv = orbit.getPVCoordinates();
1619
1620
1621 final double epoch = FastMath.floor(pv.getDate().durationFrom(AbsoluteDate.J2000_EPOCH));
1622 final double offset = pv.getDate().durationFrom(AbsoluteDate.J2000_EPOCH.shiftedBy(epoch));
1623
1624 if (orbit.hasDerivatives()) {
1625
1626 this.d = new double[] {
1627 epoch, offset, orbit.getMu(),
1628 orbit.a, orbit.e, orbit.i,
1629 orbit.pa, orbit.raan, orbit.v,
1630 orbit.aDot, orbit.eDot, orbit.iDot,
1631 orbit.paDot, orbit.raanDot, orbit.vDot
1632 };
1633 } else {
1634
1635 this.d = new double[] {
1636 epoch, offset, orbit.getMu(),
1637 orbit.a, orbit.e, orbit.i,
1638 orbit.pa, orbit.raan, orbit.v
1639 };
1640 }
1641
1642 this.frame = orbit.getFrame();
1643
1644 }
1645
1646
1647
1648
1649 private Object readResolve() {
1650 if (d.length >= 15) {
1651
1652 return new KeplerianOrbit(d[ 3], d[ 4], d[ 5], d[ 6], d[ 7], d[ 8],
1653 d[ 9], d[10], d[11], d[12], d[13], d[14],
1654 PositionAngle.TRUE,
1655 frame, AbsoluteDate.J2000_EPOCH.shiftedBy(d[0]).shiftedBy(d[1]),
1656 d[2]);
1657 } else {
1658
1659 return new KeplerianOrbit(d[3], d[4], d[5], d[6], d[7], d[8], PositionAngle.TRUE,
1660 frame, AbsoluteDate.J2000_EPOCH.shiftedBy(d[0]).shiftedBy(d[1]),
1661 d[2]);
1662 }
1663 }
1664
1665 }
1666
1667 }