1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.analytical.tle;
18
19 import java.util.ArrayList;
20 import java.util.Collections;
21 import java.util.HashMap;
22 import java.util.List;
23 import java.util.Map;
24
25 import org.hipparchus.geometry.euclidean.threed.Vector3D;
26 import org.hipparchus.linear.RealMatrix;
27 import org.hipparchus.util.FastMath;
28 import org.hipparchus.util.MathUtils;
29 import org.hipparchus.util.SinCos;
30 import org.orekit.annotation.DefaultDataContext;
31 import org.orekit.attitudes.Attitude;
32 import org.orekit.attitudes.AttitudeProvider;
33 import org.orekit.attitudes.FrameAlignedProvider;
34 import org.orekit.data.DataContext;
35 import org.orekit.errors.OrekitException;
36 import org.orekit.errors.OrekitMessages;
37 import org.orekit.frames.Frame;
38 import org.orekit.orbits.CartesianOrbit;
39 import org.orekit.orbits.Orbit;
40 import org.orekit.propagation.AbstractMatricesHarvester;
41 import org.orekit.propagation.MatricesHarvester;
42 import org.orekit.propagation.SpacecraftState;
43 import org.orekit.propagation.analytical.AbstractAnalyticalPropagator;
44 import org.orekit.propagation.analytical.tle.generation.FixedPointTleGenerationAlgorithm;
45 import org.orekit.propagation.analytical.tle.generation.TleGenerationAlgorithm;
46 import org.orekit.time.AbsoluteDate;
47 import org.orekit.time.TimeScale;
48 import org.orekit.utils.DoubleArrayDictionary;
49 import org.orekit.utils.PVCoordinates;
50 import org.orekit.utils.ParameterDriver;
51 import org.orekit.utils.TimeSpanMap;
52 import org.orekit.utils.TimeSpanMap.Span;
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 abstract class TLEPropagator extends AbstractAnalyticalPropagator {
80
81
82
83
84 protected TLE tle;
85
86
87 protected final TimeScale utc;
88
89
90 protected double xnode;
91
92
93 protected double a;
94
95
96 protected double e;
97
98
99 protected double i;
100
101
102 protected double omega;
103
104
105 protected double xl;
106
107
108 protected double a0dp;
109
110
111 protected double xn0dp;
112
113
114 protected double cosi0;
115
116
117 protected double theta2;
118
119
120 protected double sini0;
121
122
123 protected double xmdot;
124
125
126 protected double omgdot;
127
128
129 protected double xnodot;
130
131
132 protected double e0sq;
133
134 protected double beta02;
135
136
137 protected double beta0;
138
139
140 protected double perige;
141
142
143 protected double etasq;
144
145
146 protected double eeta;
147
148
149 protected double s4;
150
151
152 protected double tsi;
153
154
155 protected double eta;
156
157
158 protected double coef;
159
160
161 protected double coef1;
162
163
164 protected double c1;
165
166
167 protected double c2;
168
169
170 protected double c4;
171
172
173 protected double xnodcf;
174
175
176 protected double t2cof;
177
178
179
180
181 private final Frame teme;
182
183
184 private Map<TLE, Double> masses;
185
186
187 private TimeSpanMap<TLE> tles;
188
189
190
191
192
193
194
195
196
197
198 @DefaultDataContext
199 protected TLEPropagator(final TLE initialTLE, final AttitudeProvider attitudeProvider,
200 final double mass) {
201 this(initialTLE, attitudeProvider, mass,
202 DataContext.getDefault().getFrames().getTEME());
203 }
204
205
206
207
208
209
210
211
212 protected TLEPropagator(final TLE initialTLE,
213 final AttitudeProvider attitudeProvider,
214 final double mass,
215 final Frame teme) {
216 super(attitudeProvider);
217 setStartDate(initialTLE.getDate());
218 this.utc = initialTLE.getUtc();
219 initializeTle(initialTLE);
220 this.teme = teme;
221 this.tles = new TimeSpanMap<>(tle);
222 this.masses = new HashMap<>();
223 this.masses.put(tle, mass);
224
225
226 final Orbit orbit = propagateOrbit(initialTLE.getDate());
227 final Attitude attitude = attitudeProvider.getAttitude(orbit, orbit.getDate(), orbit.getFrame());
228 super.resetInitialState(new SpacecraftState(orbit, attitude, mass));
229 }
230
231
232
233
234
235
236
237
238
239 @DefaultDataContext
240 public static TLEPropagator selectExtrapolator(final TLE tle) {
241 return selectExtrapolator(tle, DataContext.getDefault().getFrames().getTEME());
242 }
243
244
245
246
247
248
249
250
251 public static TLEPropagator selectExtrapolator(final TLE tle, final Frame teme) {
252 return selectExtrapolator(tle, teme, FrameAlignedProvider.of(teme));
253 }
254
255
256
257
258
259
260
261
262 public static TLEPropagator selectExtrapolator(final TLE tle, final Frame teme, final AttitudeProvider attitudeProvider) {
263 return selectExtrapolator(tle, attitudeProvider, DEFAULT_MASS, teme);
264 }
265
266
267
268
269
270
271
272
273
274
275
276 @DefaultDataContext
277 public static TLEPropagator selectExtrapolator(final TLE tle, final AttitudeProvider attitudeProvider,
278 final double mass) {
279 return selectExtrapolator(tle, attitudeProvider, mass,
280 DataContext.getDefault().getFrames().getTEME());
281 }
282
283
284
285
286
287
288
289
290
291 public static TLEPropagator selectExtrapolator(final TLE tle,
292 final AttitudeProvider attitudeProvider,
293 final double mass,
294 final Frame teme) {
295
296 final double a1 = FastMath.pow( TLEConstants.XKE / (tle.getMeanMotion() * 60.0), TLEConstants.TWO_THIRD);
297 final double cosi0 = FastMath.cos(tle.getI());
298 final double temp = TLEConstants.CK2 * 1.5 * (3 * cosi0 * cosi0 - 1.0) *
299 FastMath.pow(1.0 - tle.getE() * tle.getE(), -1.5);
300 final double delta1 = temp / (a1 * a1);
301 final double a0 = a1 * (1.0 - delta1 * (TLEConstants.ONE_THIRD + delta1 * (delta1 * 134.0 / 81.0 + 1.0)));
302 final double delta0 = temp / (a0 * a0);
303
304
305 final double xn0dp = tle.getMeanMotion() * 60.0 / (delta0 + 1.0);
306
307
308 if (MathUtils.TWO_PI / (xn0dp * TLEConstants.MINUTES_PER_DAY) >= (1.0 / 6.4)) {
309 return new DeepSDP4(tle, attitudeProvider, mass, teme);
310 } else {
311 return new SGP4(tle, attitudeProvider, mass, teme);
312 }
313 }
314
315
316
317
318 public static double getMU() {
319 return TLEConstants.MU;
320 }
321
322
323
324
325
326 public PVCoordinates getPVCoordinates(final AbsoluteDate date) {
327
328 sxpPropagate(date.durationFrom(tle.getDate()) / 60.0);
329
330
331 return computePVCoordinates();
332 }
333
334
335
336 private void initializeCommons() {
337
338
339 final SinCos scI0 = FastMath.sinCos(tle.getI());
340
341 final double a1 = FastMath.pow(TLEConstants.XKE / (tle.getMeanMotion() * 60.0), TLEConstants.TWO_THIRD);
342 cosi0 = scI0.cos();
343 theta2 = cosi0 * cosi0;
344 final double x3thm1 = 3.0 * theta2 - 1.0;
345 e0sq = tle.getE() * tle.getE();
346 beta02 = 1.0 - e0sq;
347 beta0 = FastMath.sqrt(beta02);
348 final double tval = TLEConstants.CK2 * 1.5 * x3thm1 / (beta0 * beta02);
349 final double delta1 = tval / (a1 * a1);
350 final double a0 = a1 * (1.0 - delta1 * (TLEConstants.ONE_THIRD + delta1 * (1.0 + 134.0 / 81.0 * delta1)));
351 final double delta0 = tval / (a0 * a0);
352
353
354 xn0dp = tle.getMeanMotion() * 60.0 / (delta0 + 1.0);
355 a0dp = a0 / (1.0 - delta0);
356
357
358 s4 = TLEConstants.S;
359 double q0ms24 = TLEConstants.QOMS2T;
360
361 perige = (a0dp * (1 - tle.getE()) - TLEConstants.NORMALIZED_EQUATORIAL_RADIUS) * TLEConstants.EARTH_RADIUS;
362
363
364 if (perige < 156.0) {
365 if (perige <= 98.0) {
366 s4 = 20.0;
367 } else {
368 s4 = perige - 78.0;
369 }
370 final double temp_val = (120.0 - s4) * TLEConstants.NORMALIZED_EQUATORIAL_RADIUS / TLEConstants.EARTH_RADIUS;
371 final double temp_val_squared = temp_val * temp_val;
372 q0ms24 = temp_val_squared * temp_val_squared;
373 s4 = s4 / TLEConstants.EARTH_RADIUS + TLEConstants.NORMALIZED_EQUATORIAL_RADIUS;
374 }
375
376 final double pinv = 1.0 / (a0dp * beta02);
377 final double pinvsq = pinv * pinv;
378 tsi = 1.0 / (a0dp - s4);
379 eta = a0dp * tle.getE() * tsi;
380 etasq = eta * eta;
381 eeta = tle.getE() * eta;
382
383 final double psisq = FastMath.abs(1.0 - etasq);
384 final double tsi_squared = tsi * tsi;
385 coef = q0ms24 * tsi_squared * tsi_squared;
386 coef1 = coef / FastMath.pow(psisq, 3.5);
387
388
389 c2 = coef1 * xn0dp * (a0dp * (1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) +
390 0.75 * TLEConstants.CK2 * tsi / psisq * x3thm1 * (8.0 + 3.0 * etasq * (8.0 + etasq)));
391 c1 = tle.getBStar() * c2;
392 sini0 = scI0.sin();
393
394 final double x1mth2 = 1.0 - theta2;
395
396
397 c4 = 2.0 * xn0dp * coef1 * a0dp * beta02 * (eta * (2.0 + 0.5 * etasq) +
398 tle.getE() * (0.5 + 2.0 * etasq) -
399 2 * TLEConstants.CK2 * tsi / (a0dp * psisq) *
400 (-3.0 * x3thm1 * (1.0 - 2.0 * eeta + etasq * (1.5 - 0.5 * eeta)) +
401 0.75 * x1mth2 * (2.0 * etasq - eeta * (1.0 + etasq)) * FastMath.cos(2.0 * tle.getPerigeeArgument())));
402
403 final double theta4 = theta2 * theta2;
404 final double temp1 = 3 * TLEConstants.CK2 * pinvsq * xn0dp;
405 final double temp2 = temp1 * TLEConstants.CK2 * pinvsq;
406 final double temp3 = 1.25 * TLEConstants.CK4 * pinvsq * pinvsq * xn0dp;
407
408
409 xmdot = xn0dp +
410 0.5 * temp1 * beta0 * x3thm1 +
411 0.0625 * temp2 * beta0 * (13.0 - 78.0 * theta2 + 137.0 * theta4);
412
413 final double x1m5th = 1.0 - 5.0 * theta2;
414
415 omgdot = -0.5 * temp1 * x1m5th +
416 0.0625 * temp2 * (7.0 - 114.0 * theta2 + 395.0 * theta4) +
417 temp3 * (3.0 - 36.0 * theta2 + 49.0 * theta4);
418
419 final double xhdot1 = -temp1 * cosi0;
420
421 xnodot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * theta2) + 2.0 * temp3 * (3.0 - 7.0 * theta2)) * cosi0;
422 xnodcf = 3.5 * beta02 * xhdot1 * c1;
423 t2cof = 1.5 * c1;
424
425 }
426
427
428
429
430 private PVCoordinates computePVCoordinates() {
431
432
433 final SinCos scOmega = FastMath.sinCos(omega);
434
435
436 final double axn = e * scOmega.cos();
437 double temp = 1.0 / (a * (1.0 - e * e));
438 final double xlcof = 0.125 * TLEConstants.A3OVK2 * sini0 * (3.0 + 5.0 * cosi0) / (1.0 + cosi0);
439 final double aycof = 0.25 * TLEConstants.A3OVK2 * sini0;
440 final double xll = temp * xlcof * axn;
441 final double aynl = temp * aycof;
442 final double xlt = xl + xll;
443 final double ayn = e * scOmega.sin() + aynl;
444 final double elsq = axn * axn + ayn * ayn;
445 final double capu = MathUtils.normalizeAngle(xlt - xnode, FastMath.PI);
446 double epw = capu;
447 double ecosE = 0;
448 double esinE = 0;
449 double sinEPW = 0;
450 double cosEPW = 0;
451
452
453 final double cosi0Sq = cosi0 * cosi0;
454 final double x3thm1 = 3.0 * cosi0Sq - 1.0;
455 final double x1mth2 = 1.0 - cosi0Sq;
456 final double x7thm1 = 7.0 * cosi0Sq - 1.0;
457
458 if (e > (1 - 1e-6)) {
459 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL, e);
460 }
461
462
463 final double newtonRaphsonEpsilon = 1e-12;
464 for (int j = 0; j < 10; j++) {
465
466 boolean doSecondOrderNewtonRaphson = true;
467
468 final SinCos scEPW = FastMath.sinCos(epw);
469 sinEPW = scEPW.sin();
470 cosEPW = scEPW.cos();
471 ecosE = axn * cosEPW + ayn * sinEPW;
472 esinE = axn * sinEPW - ayn * cosEPW;
473 final double f = capu - epw + esinE;
474 if (FastMath.abs(f) < newtonRaphsonEpsilon) {
475 break;
476 }
477 final double fdot = 1.0 - ecosE;
478 double delta_epw = f / fdot;
479 if (j == 0) {
480 final double maxNewtonRaphson = 1.25 * FastMath.abs(e);
481 doSecondOrderNewtonRaphson = false;
482 if (delta_epw > maxNewtonRaphson) {
483 delta_epw = maxNewtonRaphson;
484 } else if (delta_epw < -maxNewtonRaphson) {
485 delta_epw = -maxNewtonRaphson;
486 } else {
487 doSecondOrderNewtonRaphson = true;
488 }
489 }
490 if (doSecondOrderNewtonRaphson) {
491 delta_epw = f / (fdot + 0.5 * esinE * delta_epw);
492 }
493 epw += delta_epw;
494 }
495
496
497 temp = 1.0 - elsq;
498 final double pl = a * temp;
499 final double r = a * (1.0 - ecosE);
500 double temp2 = a / r;
501 final double betal = FastMath.sqrt(temp);
502 temp = esinE / (1.0 + betal);
503 final double cosu = temp2 * (cosEPW - axn + ayn * temp);
504 final double sinu = temp2 * (sinEPW - ayn - axn * temp);
505 final double u = FastMath.atan2(sinu, cosu);
506 final double sin2u = 2.0 * sinu * cosu;
507 final double cos2u = 2.0 * cosu * cosu - 1.0;
508 final double temp1 = TLEConstants.CK2 / pl;
509 temp2 = temp1 / pl;
510
511
512 final double rk = r * (1.0 - 1.5 * temp2 * betal * x3thm1) + 0.5 * temp1 * x1mth2 * cos2u;
513 final double uk = u - 0.25 * temp2 * x7thm1 * sin2u;
514 final double xnodek = xnode + 1.5 * temp2 * cosi0 * sin2u;
515 final double xinck = i + 1.5 * temp2 * cosi0 * sini0 * cos2u;
516
517
518 final SinCos scuk = FastMath.sinCos(uk);
519 final SinCos scik = FastMath.sinCos(xinck);
520 final SinCos scnok = FastMath.sinCos(xnodek);
521 final double sinuk = scuk.sin();
522 final double cosuk = scuk.cos();
523 final double sinik = scik.sin();
524 final double cosik = scik.cos();
525 final double sinnok = scnok.sin();
526 final double cosnok = scnok.cos();
527 final double xmx = -sinnok * cosik;
528 final double xmy = cosnok * cosik;
529 final double ux = xmx * sinuk + cosnok * cosuk;
530 final double uy = xmy * sinuk + sinnok * cosuk;
531 final double uz = sinik * sinuk;
532
533
534 final double cr = 1000 * rk * TLEConstants.EARTH_RADIUS;
535 final Vector3D pos = new Vector3D(cr * ux, cr * uy, cr * uz);
536
537 final double rdot = TLEConstants.XKE * FastMath.sqrt(a) * esinE / r;
538 final double rfdot = TLEConstants.XKE * FastMath.sqrt(pl) / r;
539 final double xn = TLEConstants.XKE / (a * FastMath.sqrt(a));
540 final double rdotk = rdot - xn * temp1 * x1mth2 * sin2u;
541 final double rfdotk = rfdot + xn * temp1 * (x1mth2 * cos2u + 1.5 * x3thm1);
542 final double vx = xmx * cosuk - cosnok * sinuk;
543 final double vy = xmy * cosuk - sinnok * sinuk;
544 final double vz = sinik * cosuk;
545
546 final double cv = 1000.0 * TLEConstants.EARTH_RADIUS / 60.0;
547 final Vector3D vel = new Vector3D(cv * (rdotk * ux + rfdotk * vx),
548 cv * (rdotk * uy + rfdotk * vy),
549 cv * (rdotk * uz + rfdotk * vz));
550
551 return new PVCoordinates(pos, vel);
552
553 }
554
555
556
557 protected abstract void sxpInitialize();
558
559
560
561
562 protected abstract void sxpPropagate(double t);
563
564
565
566
567
568
569
570
571
572 public void resetInitialState(final SpacecraftState state) {
573 super.resetInitialState(state);
574 resetTle(state);
575 masses = new HashMap<>();
576 masses.put(tle, state.getMass());
577 tles = new TimeSpanMap<>(tle);
578 }
579
580
581 protected void resetIntermediateState(final SpacecraftState state, final boolean forward) {
582 resetTle(state);
583 if (forward) {
584 tles.addValidAfter(tle, state.getDate(), false);
585 } else {
586 tles.addValidBefore(tle, state.getDate(), false);
587 }
588 stateChanged(state);
589 masses.put(tle, state.getMass());
590 }
591
592
593
594
595 private void resetTle(final SpacecraftState state) {
596 final TleGenerationAlgorithm algorithm = getDefaultTleGenerationAlgorithm(utc, teme);
597 final TLE newTle = algorithm.generate(state, tle);
598 initializeTle(newTle);
599 }
600
601
602
603
604 private void initializeTle(final TLE newTle) {
605 tle = newTle;
606 initializeCommons();
607 sxpInitialize();
608 }
609
610
611 protected double getMass(final AbsoluteDate date) {
612 return masses.get(tles.get(date));
613 }
614
615
616 protected Orbit propagateOrbit(final AbsoluteDate date) {
617 final TLE closestTle = tles.get(date);
618 if (!tle.equals(closestTle)) {
619 initializeTle(closestTle);
620 }
621 return new CartesianOrbit(getPVCoordinates(date), teme, date, TLEConstants.MU);
622 }
623
624
625
626
627
628
629 public TLE getTLE() {
630 return tle;
631 }
632
633
634 public Frame getFrame() {
635 return teme;
636 }
637
638
639 @Override
640 protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
641 final DoubleArrayDictionary initialJacobianColumns) {
642
643 final TLEHarvester harvester = new TLEHarvester(this, stmName, initialStm, initialJacobianColumns);
644
645 addAdditionalStateProvider(harvester);
646
647 return harvester;
648 }
649
650
651
652
653
654 protected List<String> getJacobiansColumnsNames() {
655 final List<String> columnsNames = new ArrayList<>();
656 for (final ParameterDriver driver : tle.getParametersDrivers()) {
657
658 if (driver.isSelected() && !columnsNames.contains(driver.getNamesSpanMap().getFirstSpan().getData())) {
659
660
661 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
662 columnsNames.add(span.getData());
663 }
664 }
665 }
666 Collections.sort(columnsNames);
667 return columnsNames;
668 }
669
670
671
672
673
674
675
676
677 public static TleGenerationAlgorithm getDefaultTleGenerationAlgorithm(final TimeScale utc, final Frame teme) {
678 return new FixedPointTleGenerationAlgorithm(FixedPointTleGenerationAlgorithm.EPSILON_DEFAULT,
679 FixedPointTleGenerationAlgorithm.MAX_ITERATIONS_DEFAULT,
680 FixedPointTleGenerationAlgorithm.SCALE_DEFAULT, utc, teme);
681 }
682
683 }