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