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.Collection;
21
22 import org.apache.commons.math3.analysis.interpolation.HermiteInterpolator;
23 import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
24 import org.apache.commons.math3.util.FastMath;
25 import org.apache.commons.math3.util.MathUtils;
26 import org.orekit.errors.OrekitException;
27 import org.orekit.errors.OrekitMessages;
28 import org.orekit.frames.Frame;
29 import org.orekit.time.AbsoluteDate;
30 import org.orekit.utils.PVCoordinates;
31 import org.orekit.utils.TimeStampedPVCoordinates;
32
33
34
35
36
37
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 public class CircularOrbit
74 extends Orbit {
75
76
77 private static final long serialVersionUID = 20141228L;
78
79
80 private final double a;
81
82
83 private final double ex;
84
85
86 private final double ey;
87
88
89 private final double i;
90
91
92 private final double raan;
93
94
95 private final double alphaV;
96
97
98 private final boolean serializePV;
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115 public CircularOrbit(final double a, final double ex, final double ey,
116 final double i, final double raan,
117 final double alpha, final PositionAngle type,
118 final Frame frame, final AbsoluteDate date, final double mu)
119 throws IllegalArgumentException {
120 super(frame, date, mu);
121 if (ex * ex + ey * ey >= 1.0) {
122 throw OrekitException.createIllegalArgumentException(
123 OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS, getClass().getName());
124 }
125 this.a = a;
126 this.ex = ex;
127 this.ey = ey;
128 this.i = i;
129 this.raan = raan;
130
131 switch (type) {
132 case MEAN :
133 this.alphaV = eccentricToTrue(meanToEccentric(alpha));
134 break;
135 case ECCENTRIC :
136 this.alphaV = eccentricToTrue(alpha);
137 break;
138 case TRUE :
139 this.alphaV = alpha;
140 break;
141 default :
142 throw OrekitException.createInternalError(null);
143 }
144
145 serializePV = false;
146
147 }
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164 public CircularOrbit(final double a, final double ex, final double ey,
165 final double i, final double raan,
166 final double alpha, final PositionAngle type,
167 final TimeStampedPVCoordinates pvCoordinates, final Frame frame,
168 final double mu)
169 throws IllegalArgumentException {
170 super(pvCoordinates, frame, mu);
171 if (ex * ex + ey * ey >= 1.0) {
172 throw OrekitException.createIllegalArgumentException(
173 OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS, getClass().getName());
174 }
175 this.a = a;
176 this.ex = ex;
177 this.ey = ey;
178 this.i = i;
179 this.raan = raan;
180
181 switch (type) {
182 case MEAN :
183 this.alphaV = eccentricToTrue(meanToEccentric(alpha));
184 break;
185 case ECCENTRIC :
186 this.alphaV = eccentricToTrue(alpha);
187 break;
188 case TRUE :
189 this.alphaV = alpha;
190 break;
191 default :
192 throw OrekitException.createInternalError(null);
193 }
194
195 serializePV = true;
196
197 }
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213 public CircularOrbit(final TimeStampedPVCoordinates pvCoordinates, final Frame frame, final double mu)
214 throws IllegalArgumentException {
215 super(pvCoordinates, frame, mu);
216
217
218 final Vector3D pvP = pvCoordinates.getPosition();
219 final Vector3D pvV = pvCoordinates.getVelocity();
220 final double r = pvP.getNorm();
221 final double V2 = pvV.getNormSq();
222 final double rV2OnMu = r * V2 / mu;
223
224 if (rV2OnMu > 2) {
225 throw OrekitException.createIllegalArgumentException(
226 OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS, getClass().getName());
227 }
228
229 a = r / (2 - rV2OnMu);
230
231
232 final Vector3D momentum = pvCoordinates.getMomentum();
233 i = Vector3D.angle(momentum, Vector3D.PLUS_K);
234
235
236 final Vector3D node = Vector3D.crossProduct(Vector3D.PLUS_K, momentum);
237 raan = FastMath.atan2(node.getY(), node.getX());
238
239
240 final double cosRaan = FastMath.cos(raan);
241 final double sinRaan = FastMath.sin(raan);
242 final double cosI = FastMath.cos(i);
243 final double sinI = FastMath.sin(i);
244 final double xP = pvP.getX();
245 final double yP = pvP.getY();
246 final double zP = pvP.getZ();
247 final double x2 = (xP * cosRaan + yP * sinRaan) / a;
248 final double y2 = ((yP * cosRaan - xP * sinRaan) * cosI + zP * sinI) / a;
249
250
251 final double eSE = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(mu * a);
252 final double eCE = rV2OnMu - 1;
253 final double e2 = eCE * eCE + eSE * eSE;
254 final double f = eCE - e2;
255 final double g = FastMath.sqrt(1 - e2) * eSE;
256 final double aOnR = a / r;
257 final double a2OnR2 = aOnR * aOnR;
258 ex = a2OnR2 * (f * x2 + g * y2);
259 ey = a2OnR2 * (f * y2 - g * x2);
260
261
262 final double beta = 1 / (1 + FastMath.sqrt(1 - ex * ex - ey * ey));
263 alphaV = eccentricToTrue(FastMath.atan2(y2 + ey + eSE * beta * ex, x2 + ex - eSE * beta * ey));
264
265 serializePV = true;
266
267 }
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284 public CircularOrbit(final PVCoordinates pvCoordinates, final Frame frame,
285 final AbsoluteDate date, final double mu)
286 throws IllegalArgumentException {
287 this(new TimeStampedPVCoordinates(date, pvCoordinates), frame, mu);
288 }
289
290
291
292
293 public CircularOrbit(final Orbit op) {
294 super(op.getFrame(), op.getDate(), op.getMu());
295 a = op.getA();
296 i = op.getI();
297 raan = FastMath.atan2(op.getHy(), op.getHx());
298 final double cosRaan = FastMath.cos(raan);
299 final double sinRaan = FastMath.sin(raan);
300 final double equiEx = op.getEquinoctialEx();
301 final double equiEy = op.getEquinoctialEy();
302 ex = equiEx * cosRaan + equiEy * sinRaan;
303 ey = equiEy * cosRaan - equiEx * sinRaan;
304 this.alphaV = op.getLv() - raan;
305 serializePV = false;
306 }
307
308
309 public OrbitType getType() {
310 return OrbitType.CIRCULAR;
311 }
312
313
314 public double getA() {
315 return a;
316 }
317
318
319 public double getEquinoctialEx() {
320 return ex * FastMath.cos(raan) - ey * FastMath.sin(raan);
321 }
322
323
324 public double getEquinoctialEy() {
325 return ey * FastMath.cos(raan) + ex * FastMath.sin(raan);
326 }
327
328
329
330
331 public double getCircularEx() {
332 return ex;
333 }
334
335
336
337
338 public double getCircularEy() {
339 return ey;
340 }
341
342
343 public double getHx() {
344 return FastMath.cos(raan) * FastMath.tan(i / 2);
345 }
346
347
348 public double getHy() {
349 return FastMath.sin(raan) * FastMath.tan(i / 2);
350 }
351
352
353
354
355 public double getAlphaV() {
356 return alphaV;
357 }
358
359
360
361
362
363 public double getAlpha(final PositionAngle type) {
364 return (type == PositionAngle.MEAN) ? getAlphaM() :
365 ((type == PositionAngle.ECCENTRIC) ? getAlphaE() :
366 getAlphaV());
367 }
368
369
370
371
372 public double getAlphaE() {
373 final double epsilon = FastMath.sqrt(1 - ex * ex - ey * ey);
374 final double cosAlphaV = FastMath.cos(alphaV);
375 final double sinAlphaV = FastMath.sin(alphaV);
376 return alphaV + 2 * FastMath.atan((ey * cosAlphaV - ex * sinAlphaV) /
377 (epsilon + 1 + ex * cosAlphaV + ey * sinAlphaV));
378 }
379
380
381
382
383
384 private double eccentricToTrue(final double alphaE) {
385 final double epsilon = FastMath.sqrt(1 - ex * ex - ey * ey);
386 final double cosAlphaE = FastMath.cos(alphaE);
387 final double sinAlphaE = FastMath.sin(alphaE);
388 return alphaE + 2 * FastMath.atan((ex * sinAlphaE - ey * cosAlphaE) /
389 (epsilon + 1 - ex * cosAlphaE - ey * sinAlphaE));
390 }
391
392
393
394
395 public double getAlphaM() {
396 final double alphaE = getAlphaE();
397 return alphaE - ex * FastMath.sin(alphaE) + ey * FastMath.cos(alphaE);
398 }
399
400
401
402
403
404 private double meanToEccentric(final double alphaM) {
405
406
407
408 double alphaE = alphaM;
409 double shift = 0.0;
410 double alphaEMalphaM = 0.0;
411 double cosAlphaE = FastMath.cos(alphaE);
412 double sinAlphaE = FastMath.sin(alphaE);
413 int iter = 0;
414 do {
415 final double f2 = ex * sinAlphaE - ey * cosAlphaE;
416 final double f1 = 1.0 - ex * cosAlphaE - ey * sinAlphaE;
417 final double f0 = alphaEMalphaM - f2;
418
419 final double f12 = 2.0 * f1;
420 shift = f0 * f12 / (f1 * f12 - f0 * f2);
421
422 alphaEMalphaM -= shift;
423 alphaE = alphaM + alphaEMalphaM;
424 cosAlphaE = FastMath.cos(alphaE);
425 sinAlphaE = FastMath.sin(alphaE);
426
427 } while ((++iter < 50) && (FastMath.abs(shift) > 1.0e-12));
428
429 return alphaE;
430
431 }
432
433
434 public double getE() {
435 return FastMath.sqrt(ex * ex + ey * ey);
436 }
437
438
439 public double getI() {
440 return i;
441 }
442
443
444
445
446 public double getRightAscensionOfAscendingNode() {
447 return raan;
448 }
449
450
451 public double getLv() {
452 return alphaV + raan;
453 }
454
455
456 public double getLE() {
457 return getAlphaE() + raan;
458 }
459
460
461 public double getLM() {
462 return getAlphaM() + raan;
463 }
464
465
466 protected TimeStampedPVCoordinates initPVCoordinates() {
467
468
469 final double equEx = getEquinoctialEx();
470 final double equEy = getEquinoctialEy();
471 final double hx = getHx();
472 final double hy = getHy();
473 final double lE = getLE();
474
475
476 final double hx2 = hx * hx;
477 final double hy2 = hy * hy;
478 final double factH = 1. / (1 + hx2 + hy2);
479
480
481 final double ux = (1 + hx2 - hy2) * factH;
482 final double uy = 2 * hx * hy * factH;
483 final double uz = -2 * hy * factH;
484
485 final double vx = uy;
486 final double vy = (1 - hx2 + hy2) * factH;
487 final double vz = 2 * hx * factH;
488
489
490 final double exey = equEx * equEy;
491 final double ex2 = equEx * equEx;
492 final double ey2 = equEy * equEy;
493 final double e2 = ex2 + ey2;
494 final double eta = 1 + FastMath.sqrt(1 - e2);
495 final double beta = 1. / eta;
496
497
498 final double cLe = FastMath.cos(lE);
499 final double sLe = FastMath.sin(lE);
500 final double exCeyS = equEx * cLe + equEy * sLe;
501
502
503 final double x = a * ((1 - beta * ey2) * cLe + beta * exey * sLe - equEx);
504 final double y = a * ((1 - beta * ex2) * sLe + beta * exey * cLe - equEy);
505
506 final double factor = FastMath.sqrt(getMu() / a) / (1 - exCeyS);
507 final double xdot = factor * (-sLe + beta * equEy * exCeyS);
508 final double ydot = factor * ( cLe - beta * equEx * exCeyS);
509
510 final Vector3D position =
511 new Vector3D(x * ux + y * vx, x * uy + y * vy, x * uz + y * vz);
512 final double r2 = position.getNormSq();
513 final Vector3D velocity =
514 new Vector3D(xdot * ux + ydot * vx, xdot * uy + ydot * vy, xdot * uz + ydot * vz);
515 final Vector3D acceleration = new Vector3D(-getMu() / (r2 * FastMath.sqrt(r2)), position);
516 return new TimeStampedPVCoordinates(getDate(), position, velocity, acceleration);
517
518 }
519
520
521 public CircularOrbit shiftedBy(final double dt) {
522 return new CircularOrbit(a, ex, ey, i, raan,
523 getAlphaM() + getKeplerianMeanMotion() * dt,
524 PositionAngle.MEAN, getFrame(),
525 getDate().shiftedBy(dt), getMu());
526 }
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548 public CircularOrbit interpolate(final AbsoluteDate date, final Collection<Orbit> sample) {
549
550
551 final HermiteInterpolator interpolator = new HermiteInterpolator();
552
553
554 AbsoluteDate previousDate = null;
555 double previousRAAN = Double.NaN;
556 double previousAlphaM = Double.NaN;
557 for (final Orbit orbit : sample) {
558 final CircularOrbit circ = (CircularOrbit) OrbitType.CIRCULAR.convertType(orbit);
559 final double continuousRAAN;
560 final double continuousAlphaM;
561 if (previousDate == null) {
562 continuousRAAN = circ.getRightAscensionOfAscendingNode();
563 continuousAlphaM = circ.getAlphaM();
564 } else {
565 final double dt = circ.getDate().durationFrom(previousDate);
566 final double keplerAM = previousAlphaM + circ.getKeplerianMeanMotion() * dt;
567 continuousRAAN = MathUtils.normalizeAngle(circ.getRightAscensionOfAscendingNode(), previousRAAN);
568 continuousAlphaM = MathUtils.normalizeAngle(circ.getAlphaM(), keplerAM);
569 }
570 previousDate = circ.getDate();
571 previousRAAN = continuousRAAN;
572 previousAlphaM = continuousAlphaM;
573 interpolator.addSamplePoint(circ.getDate().durationFrom(date),
574 new double[] {
575 circ.getA(),
576 circ.getCircularEx(),
577 circ.getCircularEy(),
578 circ.getI(),
579 continuousRAAN,
580 continuousAlphaM
581 });
582 }
583
584
585 final double[] interpolated = interpolator.value(0);
586
587
588 return new CircularOrbit(interpolated[0], interpolated[1], interpolated[2],
589 interpolated[3], interpolated[4], interpolated[5],
590 PositionAngle.MEAN, getFrame(), date, getMu());
591
592 }
593
594
595 protected double[][] computeJacobianMeanWrtCartesian() {
596
597 final double[][] jacobian = new double[6][6];
598
599
600 final PVCoordinates pvc = getPVCoordinates();
601 final Vector3D position = pvc.getPosition();
602 final Vector3D velocity = pvc.getVelocity();
603
604 final double x = position.getX();
605 final double y = position.getY();
606 final double z = position.getZ();
607 final double vx = velocity.getX();
608 final double vy = velocity.getY();
609 final double vz = velocity.getZ();
610 final double pv = Vector3D.dotProduct(position, velocity);
611 final double r2 = position.getNormSq();
612 final double r = FastMath.sqrt(r2);
613 final double v2 = velocity.getNormSq();
614
615 final double mu = getMu();
616 final double oOsqrtMuA = 1 / FastMath.sqrt(mu * a);
617 final double rOa = r / a;
618 final double aOr = a / r;
619 final double aOr2 = a / r2;
620 final double a2 = a * a;
621
622 final double ex2 = ex * ex;
623 final double ey2 = ey * ey;
624 final double e2 = ex2 + ey2;
625 final double epsilon = FastMath.sqrt(1 - e2);
626 final double beta = 1 / (1 + epsilon);
627
628 final double eCosE = 1 - rOa;
629 final double eSinE = pv * oOsqrtMuA;
630
631 final double cosI = FastMath.cos(i);
632 final double sinI = FastMath.sin(i);
633 final double cosRaan = FastMath.cos(raan);
634 final double sinRaan = FastMath.sin(raan);
635
636
637 fillHalfRow(2 * aOr * aOr2, position, jacobian[0], 0);
638 fillHalfRow(2 * a2 / mu, velocity, jacobian[0], 3);
639
640
641 final Vector3D danP = new Vector3D(v2, position, -pv, velocity);
642 final Vector3D danV = new Vector3D(r2, velocity, -pv, position);
643 final double recip = 1 / pvc.getMomentum().getNorm();
644 final double recip2 = recip * recip;
645 final Vector3D dwXP = new Vector3D(recip, new Vector3D( 0, vz, -vy), -recip2 * sinRaan * sinI, danP);
646 final Vector3D dwYP = new Vector3D(recip, new Vector3D(-vz, 0, vx), recip2 * cosRaan * sinI, danP);
647 final Vector3D dwZP = new Vector3D(recip, new Vector3D( vy, -vx, 0), -recip2 * cosI, danP);
648 final Vector3D dwXV = new Vector3D(recip, new Vector3D( 0, -z, y), -recip2 * sinRaan * sinI, danV);
649 final Vector3D dwYV = new Vector3D(recip, new Vector3D( z, 0, -x), recip2 * cosRaan * sinI, danV);
650 final Vector3D dwZV = new Vector3D(recip, new Vector3D( -y, x, 0), -recip2 * cosI, danV);
651
652
653 fillHalfRow(sinRaan * cosI, dwXP, -cosRaan * cosI, dwYP, -sinI, dwZP, jacobian[3], 0);
654 fillHalfRow(sinRaan * cosI, dwXV, -cosRaan * cosI, dwYV, -sinI, dwZV, jacobian[3], 3);
655
656
657 fillHalfRow(sinRaan / sinI, dwYP, cosRaan / sinI, dwXP, jacobian[4], 0);
658 fillHalfRow(sinRaan / sinI, dwYV, cosRaan / sinI, dwXV, jacobian[4], 3);
659
660
661
662 final double u = x * cosRaan + y * sinRaan;
663 final double cv = -x * sinRaan + y * cosRaan;
664 final double v = cv * cosI + z * sinI;
665
666
667 final Vector3D duP = new Vector3D(cv * cosRaan / sinI, dwXP,
668 cv * sinRaan / sinI, dwYP,
669 1, new Vector3D(cosRaan, sinRaan, 0));
670 final Vector3D duV = new Vector3D(cv * cosRaan / sinI, dwXV,
671 cv * sinRaan / sinI, dwYV);
672
673
674 final Vector3D dvP = new Vector3D(-u * cosRaan * cosI / sinI + sinRaan * z, dwXP,
675 -u * sinRaan * cosI / sinI - cosRaan * z, dwYP,
676 cv, dwZP,
677 1, new Vector3D(-sinRaan * cosI, cosRaan * cosI, sinI));
678 final Vector3D dvV = new Vector3D(-u * cosRaan * cosI / sinI + sinRaan * z, dwXV,
679 -u * sinRaan * cosI / sinI - cosRaan * z, dwYV,
680 cv, dwZV);
681
682 final Vector3D dc1P = new Vector3D(aOr2 * (2 * eSinE * eSinE + 1 - eCosE) / r2, position,
683 -2 * aOr2 * eSinE * oOsqrtMuA, velocity);
684 final Vector3D dc1V = new Vector3D(-2 * aOr2 * eSinE * oOsqrtMuA, position,
685 2 / mu, velocity);
686 final Vector3D dc2P = new Vector3D(aOr2 * eSinE * (eSinE * eSinE - (1 - e2)) / (r2 * epsilon), position,
687 aOr2 * (1 - e2 - eSinE * eSinE) * oOsqrtMuA / epsilon, velocity);
688 final Vector3D dc2V = new Vector3D(aOr2 * (1 - e2 - eSinE * eSinE) * oOsqrtMuA / epsilon, position,
689 eSinE / (mu * epsilon), velocity);
690
691 final double cof1 = aOr2 * (eCosE - e2);
692 final double cof2 = aOr2 * epsilon * eSinE;
693 final Vector3D dexP = new Vector3D(u, dc1P, v, dc2P, cof1, duP, cof2, dvP);
694 final Vector3D dexV = new Vector3D(u, dc1V, v, dc2V, cof1, duV, cof2, dvV);
695 final Vector3D deyP = new Vector3D(v, dc1P, -u, dc2P, cof1, dvP, -cof2, duP);
696 final Vector3D deyV = new Vector3D(v, dc1V, -u, dc2V, cof1, dvV, -cof2, duV);
697 fillHalfRow(1, dexP, jacobian[1], 0);
698 fillHalfRow(1, dexV, jacobian[1], 3);
699 fillHalfRow(1, deyP, jacobian[2], 0);
700 fillHalfRow(1, deyV, jacobian[2], 3);
701
702 final double cle = u / a + ex - eSinE * beta * ey;
703 final double sle = v / a + ey + eSinE * beta * ex;
704 final double m1 = beta * eCosE;
705 final double m2 = 1 - m1 * eCosE;
706 final double m3 = (u * ey - v * ex) + eSinE * beta * (u * ex + v * ey);
707 final double m4 = -sle + cle * eSinE * beta;
708 final double m5 = cle + sle * eSinE * beta;
709 fillHalfRow((2 * m3 / r + aOr * eSinE + m1 * eSinE * (1 + m1 - (1 + aOr) * m2) / epsilon) / r2, position,
710 (m1 * m2 / epsilon - 1) * oOsqrtMuA, velocity,
711 m4, dexP, m5, deyP, -sle / a, duP, cle / a, dvP,
712 jacobian[5], 0);
713 fillHalfRow((m1 * m2 / epsilon - 1) * oOsqrtMuA, position,
714 (2 * m3 + eSinE * a + m1 * eSinE * r * (eCosE * beta * 2 - aOr * m2) / epsilon) / mu, velocity,
715 m4, dexV, m5, deyV, -sle / a, duV, cle / a, dvV,
716 jacobian[5], 3);
717
718 return jacobian;
719
720 }
721
722
723 protected double[][] computeJacobianEccentricWrtCartesian() {
724
725
726 final double[][] jacobian = computeJacobianMeanWrtCartesian();
727
728
729
730
731
732 final double alphaE = getAlphaE();
733 final double cosAe = FastMath.cos(alphaE);
734 final double sinAe = FastMath.sin(alphaE);
735 final double aOr = 1 / (1 - ex * cosAe - ey * sinAe);
736
737
738 final double[] rowEx = jacobian[1];
739 final double[] rowEy = jacobian[2];
740 final double[] rowL = jacobian[5];
741 for (int j = 0; j < 6; ++j) {
742 rowL[j] = aOr * (rowL[j] + sinAe * rowEx[j] - cosAe * rowEy[j]);
743 }
744
745 return jacobian;
746
747 }
748
749
750 protected double[][] computeJacobianTrueWrtCartesian() {
751
752
753 final double[][] jacobian = computeJacobianEccentricWrtCartesian();
754
755
756
757
758
759
760
761
762
763
764
765
766
767 final double alphaE = getAlphaE();
768 final double cosAe = FastMath.cos(alphaE);
769 final double sinAe = FastMath.sin(alphaE);
770 final double eSinE = ex * sinAe - ey * cosAe;
771 final double ecosE = ex * cosAe + ey * sinAe;
772 final double e2 = ex * ex + ey * ey;
773 final double epsilon = FastMath.sqrt(1 - e2);
774 final double onePeps = 1 + epsilon;
775 final double d = onePeps - ecosE;
776 final double cT = (d * d + eSinE * eSinE) / 2;
777 final double cE = ecosE * onePeps - e2;
778 final double cX = ex * eSinE / epsilon - ey + sinAe * onePeps;
779 final double cY = ey * eSinE / epsilon + ex - cosAe * onePeps;
780 final double factorLe = (cT + cE) / cT;
781 final double factorEx = cX / cT;
782 final double factorEy = cY / cT;
783
784
785 final double[] rowEx = jacobian[1];
786 final double[] rowEy = jacobian[2];
787 final double[] rowA = jacobian[5];
788 for (int j = 0; j < 6; ++j) {
789 rowA[j] = factorLe * rowA[j] + factorEx * rowEx[j] + factorEy * rowEy[j];
790 }
791
792 return jacobian;
793
794 }
795
796
797 public void addKeplerContribution(final PositionAngle type, final double gm,
798 final double[] pDot) {
799 final double oMe2;
800 final double ksi;
801 final double n = FastMath.sqrt(gm / a) / a;
802 switch (type) {
803 case MEAN :
804 pDot[5] += n;
805 break;
806 case ECCENTRIC :
807 oMe2 = 1 - ex * ex - ey * ey;
808 ksi = 1 + ex * FastMath.cos(alphaV) + ey * FastMath.sin(alphaV);
809 pDot[5] += n * ksi / oMe2;
810 break;
811 case TRUE :
812 oMe2 = 1 - ex * ex - ey * ey;
813 ksi = 1 + ex * FastMath.cos(alphaV) + ey * FastMath.sin(alphaV);
814 pDot[5] += n * ksi * ksi / (oMe2 * FastMath.sqrt(oMe2));
815 break;
816 default :
817 throw OrekitException.createInternalError(null);
818 }
819 }
820
821
822
823
824 public String toString() {
825 return new StringBuffer().append("circular parameters: ").append('{').
826 append("a: ").append(a).
827 append(", ex: ").append(ex).append(", ey: ").append(ey).
828 append(", i: ").append(FastMath.toDegrees(i)).
829 append(", raan: ").append(FastMath.toDegrees(raan)).
830 append(", alphaV: ").append(FastMath.toDegrees(alphaV)).
831 append(";}").toString();
832 }
833
834
835
836
837 private Object writeReplace() {
838 return new DTO(this);
839 }
840
841
842 private static class DTO implements Serializable {
843
844
845 private static final long serialVersionUID = 20140617L;
846
847
848 private double[] d;
849
850
851 private final Frame frame;
852
853
854
855
856 private DTO(final CircularOrbit orbit) {
857
858 final AbsoluteDate date = orbit.getDate();
859
860
861 final double epoch = FastMath.floor(date.durationFrom(AbsoluteDate.J2000_EPOCH));
862 final double offset = date.durationFrom(AbsoluteDate.J2000_EPOCH.shiftedBy(epoch));
863
864 if (orbit.serializePV) {
865 final TimeStampedPVCoordinates pv = orbit.getPVCoordinates();
866 this.d = new double[] {
867 epoch, offset, orbit.getMu(),
868 orbit.a, orbit.ex, orbit.ey,
869 orbit.i, orbit.raan, orbit.alphaV,
870 pv.getPosition().getX(), pv.getPosition().getY(), pv.getPosition().getZ(),
871 pv.getVelocity().getX(), pv.getVelocity().getY(), pv.getVelocity().getZ(),
872 pv.getAcceleration().getX(), pv.getAcceleration().getY(), pv.getAcceleration().getZ(),
873 };
874 } else {
875 this.d = new double[] {
876 epoch, offset, orbit.getMu(),
877 orbit.a, orbit.ex, orbit.ey,
878 orbit.i, orbit.raan, orbit.alphaV
879 };
880 }
881
882 this.frame = orbit.getFrame();
883
884 }
885
886
887
888
889 private Object readResolve() {
890 if (d.length > 10) {
891 return new CircularOrbit(d[3], d[4], d[5], d[6], d[7], d[8], PositionAngle.TRUE,
892 new TimeStampedPVCoordinates(AbsoluteDate.J2000_EPOCH.shiftedBy(d[0]).shiftedBy(d[1]),
893 new Vector3D(d[9], d[10], d[11]),
894 new Vector3D(d[12], d[13], d[14]),
895 new Vector3D(d[15], d[16], d[17])),
896 frame,
897 d[2]);
898 } else {
899 return new CircularOrbit(d[3], d[4], d[5], d[6], d[7], d[8], PositionAngle.TRUE,
900 frame, AbsoluteDate.J2000_EPOCH.shiftedBy(d[0]).shiftedBy(d[1]),
901 d[2]);
902
903 }
904 }
905
906 }
907
908 }