1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.utils;
18
19 import java.io.BufferedReader;
20 import java.io.IOException;
21 import java.io.InputStream;
22 import java.io.InputStreamReader;
23 import java.io.Serializable;
24 import java.util.List;
25
26 import org.hipparchus.RealFieldElement;
27 import org.hipparchus.analysis.interpolation.FieldHermiteInterpolator;
28 import org.hipparchus.analysis.interpolation.HermiteInterpolator;
29 import org.hipparchus.util.FastMath;
30 import org.hipparchus.util.MathArrays;
31 import org.hipparchus.util.MathUtils;
32 import org.orekit.data.BodiesElements;
33 import org.orekit.data.DelaunayArguments;
34 import org.orekit.data.FieldBodiesElements;
35 import org.orekit.data.FieldDelaunayArguments;
36 import org.orekit.data.FundamentalNutationArguments;
37 import org.orekit.data.PoissonSeries;
38 import org.orekit.data.PoissonSeriesParser;
39 import org.orekit.data.PolynomialNutation;
40 import org.orekit.data.PolynomialParser;
41 import org.orekit.data.PolynomialParser.Unit;
42 import org.orekit.data.SimpleTimeStampedTableParser;
43 import org.orekit.errors.OrekitException;
44 import org.orekit.errors.OrekitInternalError;
45 import org.orekit.errors.OrekitMessages;
46 import org.orekit.errors.TimeStampedCacheException;
47 import org.orekit.frames.EOPHistory;
48 import org.orekit.frames.FieldPoleCorrection;
49 import org.orekit.frames.PoleCorrection;
50 import org.orekit.time.AbsoluteDate;
51 import org.orekit.time.DateComponents;
52 import org.orekit.time.FieldAbsoluteDate;
53 import org.orekit.time.TimeComponents;
54 import org.orekit.time.TimeScalarFunction;
55 import org.orekit.time.TimeScale;
56 import org.orekit.time.TimeScalesFactory;
57 import org.orekit.time.TimeStamped;
58 import org.orekit.time.TimeVectorFunction;
59
60
61
62
63
64
65 public enum IERSConventions {
66
67
68 IERS_1996 {
69
70
71 private static final String NUTATION_ARGUMENTS = IERS_BASE + "1996/nutation-arguments.txt";
72
73
74 private static final String X_Y_SERIES = IERS_BASE + "1996/tab5.4.txt";
75
76
77 private static final String PSI_EPSILON_SERIES = IERS_BASE + "1996/tab5.1.txt";
78
79
80 private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "1996/tab8.4.txt";
81
82
83 private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "1996/tab8.3.txt";
84
85
86 private static final String LOVE_NUMBERS = IERS_BASE + "1996/tab6.1.txt";
87
88
89 private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2b.txt";
90
91
92 private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2a.txt";
93
94
95 private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2c.txt";
96
97
98 @Override
99 public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
100 throws OrekitException {
101 return new FundamentalNutationArguments(this, timeScale,
102 getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
103 }
104
105
106 @Override
107 public TimeScalarFunction getMeanObliquityFunction() throws OrekitException {
108
109
110 final PolynomialNutation epsilonA =
111 new PolynomialNutation(84381.448 * Constants.ARC_SECONDS_TO_RADIANS,
112 -46.8150 * Constants.ARC_SECONDS_TO_RADIANS,
113 -0.00059 * Constants.ARC_SECONDS_TO_RADIANS,
114 0.001813 * Constants.ARC_SECONDS_TO_RADIANS);
115
116 return new TimeScalarFunction() {
117
118
119 @Override
120 public double value(final AbsoluteDate date) {
121 return epsilonA.value(evaluateTC(date));
122 }
123
124
125 @Override
126 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
127 return epsilonA.value(evaluateTC(date));
128 }
129
130 };
131
132 }
133
134
135 @Override
136 public TimeVectorFunction getXYSpXY2Function()
137 throws OrekitException {
138
139
140 final FundamentalNutationArguments arguments = getNutationArguments(null);
141
142
143
144
145 final PolynomialNutation xPolynomial =
146 new PolynomialNutation(0,
147 2004.3109 * Constants.ARC_SECONDS_TO_RADIANS,
148 -0.42665 * Constants.ARC_SECONDS_TO_RADIANS,
149 -0.198656 * Constants.ARC_SECONDS_TO_RADIANS,
150 0.0000140 * Constants.ARC_SECONDS_TO_RADIANS);
151
152 final double fXCosOm = 0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
153 final double fXSinOm = 0.00204 * Constants.ARC_SECONDS_TO_RADIANS;
154 final double fXSin2FDOm = 0.00016 * Constants.ARC_SECONDS_TO_RADIANS;
155 final double sinEps0 = FastMath.sin(getMeanObliquityFunction().value(getNutationReferenceEpoch()));
156
157 final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
158 final PoissonSeriesParser baseParser =
159 new PoissonSeriesParser(12).withFirstDelaunay(1);
160
161 final PoissonSeriesParser xParser =
162 baseParser.
163 withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
164 withSinCos(1, 8, deciMilliAS, 9, deciMilliAS);
165 final PoissonSeries xSum = xParser.parse(getStream(X_Y_SERIES), X_Y_SERIES);
166
167
168
169
170 final PolynomialNutation yPolynomial =
171 new PolynomialNutation(-0.00013 * Constants.ARC_SECONDS_TO_RADIANS,
172 0.0,
173 -22.40992 * Constants.ARC_SECONDS_TO_RADIANS,
174 0.001836 * Constants.ARC_SECONDS_TO_RADIANS,
175 0.0011130 * Constants.ARC_SECONDS_TO_RADIANS);
176
177 final double fYCosOm = -0.00231 * Constants.ARC_SECONDS_TO_RADIANS;
178 final double fYCos2FDOm = -0.00014 * Constants.ARC_SECONDS_TO_RADIANS;
179
180 final PoissonSeriesParser yParser =
181 baseParser.
182 withSinCos(0, -1, deciMilliAS, 10, deciMilliAS).
183 withSinCos(1, 12, deciMilliAS, 11, deciMilliAS);
184 final PoissonSeries ySum = yParser.parse(getStream(X_Y_SERIES), X_Y_SERIES);
185
186 final PoissonSeries.CompiledSeries xySum =
187 PoissonSeries.compile(xSum, ySum);
188
189
190
191 final double fST = 0.00385 * Constants.ARC_SECONDS_TO_RADIANS;
192 final double fST3 = -0.07259 * Constants.ARC_SECONDS_TO_RADIANS;
193 final double fSSinOm = -0.00264 * Constants.ARC_SECONDS_TO_RADIANS;
194 final double fSSin2Om = -0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
195 final double fST2SinOm = 0.00074 * Constants.ARC_SECONDS_TO_RADIANS;
196 final double fST2Sin2FDOm = 0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
197
198 return new TimeVectorFunction() {
199
200
201 @Override
202 public double[] value(final AbsoluteDate date) {
203
204 final BodiesElements elements = arguments.evaluateAll(date);
205 final double[] xy = xySum.value(elements);
206
207 final double omega = elements.getOmega();
208 final double f = elements.getF();
209 final double d = elements.getD();
210 final double t = elements.getTC();
211
212 final double cosOmega = FastMath.cos(omega);
213 final double sinOmega = FastMath.sin(omega);
214 final double sin2Omega = FastMath.sin(2 * omega);
215 final double cos2FDOm = FastMath.cos(2 * (f - d + omega));
216 final double sin2FDOm = FastMath.sin(2 * (f - d + omega));
217
218 final double x = xPolynomial.value(t) + sinEps0 * xy[0] +
219 t * t * (fXCosOm * cosOmega + fXSinOm * sinOmega + fXSin2FDOm * cos2FDOm);
220 final double y = yPolynomial.value(t) + xy[1] +
221 t * t * (fYCosOm * cosOmega + fYCos2FDOm * cos2FDOm);
222 final double sPxy2 = fSSinOm * sinOmega + fSSin2Om * sin2Omega +
223 t * (fST + t * (fST2SinOm * sinOmega + fST2Sin2FDOm * sin2FDOm + t * fST3));
224
225 return new double[] {
226 x, y, sPxy2
227 };
228
229 }
230
231
232 @Override
233 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
234
235 final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
236 final T[] xy = xySum.value(elements);
237
238 final T omega = elements.getOmega();
239 final T f = elements.getF();
240 final T d = elements.getD();
241 final T t = elements.getTC();
242 final T t2 = t.multiply(t);
243
244 final T cosOmega = omega.cos();
245 final T sinOmega = omega.sin();
246 final T sin2Omega = omega.multiply(2).sin();
247 final T fMDPO2 = f.subtract(d).add(omega).multiply(2);
248 final T cos2FDOm = fMDPO2.cos();
249 final T sin2FDOm = fMDPO2.sin();
250
251 final T x = xPolynomial.value(t).
252 add(xy[0].multiply(sinEps0)).
253 add(t2.multiply(cosOmega.multiply(fXCosOm).add(sinOmega.multiply(fXSinOm)).add(cos2FDOm.multiply(fXSin2FDOm))));
254 final T y = yPolynomial.value(t).
255 add(xy[1]).
256 add(t2.multiply(cosOmega.multiply(fYCosOm).add(cos2FDOm.multiply(fYCos2FDOm))));
257 final T sPxy2 = sinOmega.multiply(fSSinOm).
258 add(sin2Omega.multiply(fSSin2Om)).
259 add(t.multiply(fST3).add(sinOmega.multiply(fST2SinOm)).add(sin2FDOm.multiply(fST2Sin2FDOm)).multiply(t).add(fST).multiply(t));
260
261 final T[] a = MathArrays.buildArray(date.getField(), 3);
262 a[0] = x;
263 a[1] = y;
264 a[2] = sPxy2;
265 return a;
266
267 }
268
269 };
270
271 }
272
273
274 @Override
275 public TimeVectorFunction getPrecessionFunction() throws OrekitException {
276
277
278
279
280
281
282 final PolynomialNutation psiA =
283 new PolynomialNutation( 0.0,
284 5038.7784 * Constants.ARC_SECONDS_TO_RADIANS,
285 -1.07259 * Constants.ARC_SECONDS_TO_RADIANS,
286 -0.001147 * Constants.ARC_SECONDS_TO_RADIANS);
287 final PolynomialNutation omegaA =
288 new PolynomialNutation(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
289 0.0,
290 0.05127 * Constants.ARC_SECONDS_TO_RADIANS,
291 -0.007726 * Constants.ARC_SECONDS_TO_RADIANS);
292 final PolynomialNutation chiA =
293 new PolynomialNutation( 0.0,
294 10.5526 * Constants.ARC_SECONDS_TO_RADIANS,
295 -2.38064 * Constants.ARC_SECONDS_TO_RADIANS,
296 -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);
297
298 return new PrecessionFunction(psiA, omegaA, chiA);
299
300 }
301
302
303 @Override
304 public TimeVectorFunction getNutationFunction()
305 throws OrekitException {
306
307
308 final FundamentalNutationArguments arguments = getNutationArguments(null);
309
310
311 final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
312 final PoissonSeriesParser baseParser =
313 new PoissonSeriesParser(10).withFirstDelaunay(1);
314
315 final PoissonSeriesParser psiParser =
316 baseParser.
317 withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
318 withSinCos(1, 8, deciMilliAS, -1, deciMilliAS);
319 final PoissonSeries psiSeries = psiParser.parse(getStream(PSI_EPSILON_SERIES), PSI_EPSILON_SERIES);
320
321 final PoissonSeriesParser epsilonParser =
322 baseParser.
323 withSinCos(0, -1, deciMilliAS, 9, deciMilliAS).
324 withSinCos(1, -1, deciMilliAS, 10, deciMilliAS);
325 final PoissonSeries epsilonSeries = epsilonParser.parse(getStream(PSI_EPSILON_SERIES), PSI_EPSILON_SERIES);
326
327 final PoissonSeries.CompiledSeries psiEpsilonSeries =
328 PoissonSeries.compile(psiSeries, epsilonSeries);
329
330 return new TimeVectorFunction() {
331
332
333 @Override
334 public double[] value(final AbsoluteDate date) {
335 final BodiesElements elements = arguments.evaluateAll(date);
336 final double[] psiEpsilon = psiEpsilonSeries.value(elements);
337 return new double[] {
338 psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements)
339 };
340 }
341
342
343 @Override
344 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
345 final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
346 final T[] psiEpsilon = psiEpsilonSeries.value(elements);
347 final T[] result = MathArrays.buildArray(date.getField(), 3);
348 result[0] = psiEpsilon[0];
349 result[1] = psiEpsilon[1];
350 result[2] = IAU1994ResolutionC7.value(elements);
351 return result;
352 }
353
354 };
355
356 }
357
358
359 @Override
360 public TimeScalarFunction getGMSTFunction(final TimeScale ut1)
361 throws OrekitException {
362
363
364 final double radiansPerSecond = MathUtils.TWO_PI / Constants.JULIAN_DAY;
365
366
367
368 final AbsoluteDate gmstReference =
369 new AbsoluteDate(DateComponents.J2000_EPOCH, TimeComponents.H12, TimeScalesFactory.getTAI());
370 final double gmst0 = 24110.54841;
371 final double gmst1 = 8640184.812866;
372 final double gmst2 = 0.093104;
373 final double gmst3 = -6.2e-6;
374
375 return new TimeScalarFunction() {
376
377
378 @Override
379 public double value(final AbsoluteDate date) {
380
381
382 final double dtai = date.durationFrom(gmstReference);
383 final double tut1 = dtai + ut1.offsetFromTAI(date);
384 final double tt = tut1 / Constants.JULIAN_CENTURY;
385
386
387
388 final double sd = FastMath.IEEEremainder(tut1 + Constants.JULIAN_DAY / 2, Constants.JULIAN_DAY);
389
390
391 return ((((((tt * gmst3 + gmst2) * tt) + gmst1) * tt) + gmst0) + sd) * radiansPerSecond;
392
393 }
394
395
396 @Override
397 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
398
399
400 final T dtai = date.durationFrom(gmstReference);
401 final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
402 final T tt = tut1.divide(Constants.JULIAN_CENTURY);
403
404
405
406 final T sd = tut1.add(Constants.JULIAN_DAY / 2).remainder(Constants.JULIAN_DAY);
407
408
409 return tt.multiply(gmst3).add(gmst2).multiply(tt).add(gmst1).multiply(tt).add(gmst0).add(sd).multiply(radiansPerSecond);
410
411 }
412
413 };
414
415 }
416
417
418 @Override
419 public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1)
420 throws OrekitException {
421
422
423 final double radiansPerSecond = MathUtils.TWO_PI / Constants.JULIAN_DAY;
424
425
426
427 final AbsoluteDate gmstReference =
428 new AbsoluteDate(DateComponents.J2000_EPOCH, TimeComponents.H12, TimeScalesFactory.getTAI());
429 final double gmst1 = 8640184.812866;
430 final double gmst2 = 0.093104;
431 final double gmst3 = -6.2e-6;
432
433 return new TimeScalarFunction() {
434
435
436 @Override
437 public double value(final AbsoluteDate date) {
438
439
440 final double dtai = date.durationFrom(gmstReference);
441 final double tut1 = dtai + ut1.offsetFromTAI(date);
442 final double tt = tut1 / Constants.JULIAN_CENTURY;
443
444
445 return ((((tt * 3 * gmst3 + 2 * gmst2) * tt) + gmst1) / Constants.JULIAN_CENTURY + 1) * radiansPerSecond;
446
447 }
448
449
450 @Override
451 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
452
453
454 final T dtai = date.durationFrom(gmstReference);
455 final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
456 final T tt = tut1.divide(Constants.JULIAN_CENTURY);
457
458
459 return tt.multiply(3 * gmst3).add(2 * gmst2).multiply(tt).add(gmst1).divide(Constants.JULIAN_CENTURY).add(1).multiply(radiansPerSecond);
460
461 }
462
463 };
464
465 }
466
467
468 @Override
469 public TimeScalarFunction getGASTFunction(final TimeScale ut1, final EOPHistory eopHistory)
470 throws OrekitException {
471
472
473 final TimeScalarFunction epsilonA = getMeanObliquityFunction();
474
475
476 final TimeScalarFunction gmst = getGMSTFunction(ut1);
477
478
479 final TimeVectorFunction nutation = getNutationFunction();
480
481 return new TimeScalarFunction() {
482
483
484 @Override
485 public double value(final AbsoluteDate date) {
486
487
488 final double[] angles = nutation.value(date);
489 double deltaPsi = angles[0];
490 if (eopHistory != null) {
491 deltaPsi += eopHistory.getEquinoxNutationCorrection(date)[0];
492 }
493 final double eqe = deltaPsi * FastMath.cos(epsilonA.value(date)) + angles[2];
494
495
496 return gmst.value(date) + eqe;
497
498 }
499
500
501 @Override
502 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
503
504
505 final T[] angles = nutation.value(date);
506 T deltaPsi = angles[0];
507 if (eopHistory != null) {
508 deltaPsi = deltaPsi.add(eopHistory.getEquinoxNutationCorrection(date)[0]);
509 }
510 final T eqe = deltaPsi.multiply(epsilonA.value(date).cos()).add(angles[2]);
511
512
513 return gmst.value(date).add(eqe);
514
515 }
516
517 };
518
519 }
520
521
522 @Override
523 public TimeVectorFunction getEOPTidalCorrection()
524 throws OrekitException {
525
526
527
528
529
530
531 final FundamentalNutationArguments arguments = getNutationArguments(TimeScalesFactory.getTT());
532
533
534 final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
535 final PoissonSeriesParser xyParser = new PoissonSeriesParser(17).
536 withOptionalColumn(1).
537 withGamma(7).
538 withFirstDelaunay(2);
539 final PoissonSeries xSeries =
540 xyParser.
541 withSinCos(0, 14, milliAS, 15, milliAS).
542 parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
543 final PoissonSeries ySeries =
544 xyParser.
545 withSinCos(0, 16, milliAS, 17, milliAS).
546 parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES),
547 TIDAL_CORRECTION_XP_YP_SERIES);
548
549 final double deciMilliS = 1.0e-4;
550 final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(17).
551 withOptionalColumn(1).
552 withGamma(7).
553 withFirstDelaunay(2).
554 withSinCos(0, 16, deciMilliS, 17, deciMilliS);
555 final PoissonSeries ut1Series =
556 ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);
557
558 return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
559
560 }
561
562
563 public LoveNumbers getLoveNumbers() throws OrekitException {
564 return loadLoveNumbers(LOVE_NUMBERS);
565 }
566
567
568 public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1)
569 throws OrekitException {
570
571
572 final FundamentalNutationArguments arguments = getNutationArguments(ut1);
573
574
575 final PoissonSeriesParser k20Parser =
576 new PoissonSeriesParser(18).
577 withOptionalColumn(1).
578 withDoodson(4, 2).
579 withFirstDelaunay(10);
580 final PoissonSeriesParser k21Parser =
581 new PoissonSeriesParser(18).
582 withOptionalColumn(1).
583 withDoodson(4, 3).
584 withFirstDelaunay(10);
585 final PoissonSeriesParser k22Parser =
586 new PoissonSeriesParser(16).
587 withOptionalColumn(1).
588 withDoodson(4, 2).
589 withFirstDelaunay(10);
590
591 final double pico = 1.0e-12;
592 final PoissonSeries c20Series =
593 k20Parser.
594 withSinCos(0, 18, -pico, 16, pico).
595 parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
596 final PoissonSeries c21Series =
597 k21Parser.
598 withSinCos(0, 17, pico, 18, pico).
599 parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
600 final PoissonSeries s21Series =
601 k21Parser.
602 withSinCos(0, 18, -pico, 17, pico).
603 parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
604 final PoissonSeries c22Series =
605 k22Parser.
606 withSinCos(0, -1, pico, 16, pico).
607 parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
608 final PoissonSeries s22Series =
609 k22Parser.
610 withSinCos(0, 16, -pico, -1, pico).
611 parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
612
613 return new TideFrequencyDependenceFunction(arguments,
614 c20Series,
615 c21Series, s21Series,
616 c22Series, s22Series);
617
618 }
619
620
621 @Override
622 public double getPermanentTide() throws OrekitException {
623 return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
624 }
625
626
627 @Override
628 public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {
629
630
631 final double globalFactor = -1.348e-9 / Constants.ARC_SECONDS_TO_RADIANS;
632 final double coupling = 0.00112;
633
634 return new TimeVectorFunction() {
635
636
637 @Override
638 public double[] value(final AbsoluteDate date) {
639 final PoleCorrection pole = eopHistory.getPoleCorrection(date);
640 return new double[] {
641 globalFactor * (pole.getXp() + coupling * pole.getYp()),
642 globalFactor * (coupling * pole.getXp() - pole.getYp()),
643 };
644 }
645
646
647 @Override
648 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
649 final FieldPoleCorrection<T> pole = eopHistory.getPoleCorrection(date);
650 final T[] a = MathArrays.buildArray(date.getField(), 2);
651 a[0] = pole.getXp().add(pole.getYp().multiply(coupling)).multiply(globalFactor);
652 a[1] = pole.getXp().multiply(coupling).subtract(pole.getYp()).multiply(globalFactor);
653 return a;
654 }
655
656 };
657 }
658
659
660 @Override
661 public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory)
662 throws OrekitException {
663
664 return new TimeVectorFunction() {
665
666
667 @Override
668 public double[] value(final AbsoluteDate date) {
669
670 return new double[] {
671 0.0, 0.0
672 };
673 }
674
675
676 @Override
677 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
678
679 return MathArrays.buildArray(date.getField(), 2);
680 }
681
682 };
683 }
684
685 },
686
687
688 IERS_2003 {
689
690
691 private static final String NUTATION_ARGUMENTS = IERS_BASE + "2003/nutation-arguments.txt";
692
693
694 private static final String X_SERIES = IERS_BASE + "2003/tab5.2a.txt";
695
696
697 private static final String Y_SERIES = IERS_BASE + "2003/tab5.2b.txt";
698
699
700 private static final String S_SERIES = IERS_BASE + "2003/tab5.2c.txt";
701
702
703 private static final String LUNI_SOLAR_SERIES = IERS_BASE + "2003/tab5.3a-first-table.txt";
704
705
706 private static final String PLANETARY_SERIES = IERS_BASE + "2003/tab5.3b.txt";
707
708
709 private static final String GST_SERIES = IERS_BASE + "2003/tab5.4.txt";
710
711
712 private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "2003/tab8.2ab.txt";
713
714
715 private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "2003/tab8.3ab.txt";
716
717
718 private static final String LOVE_NUMBERS = IERS_BASE + "2003/tab6.1.txt";
719
720
721 private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3b.txt";
722
723
724 private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3a.txt";
725
726
727 private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3c.txt";
728
729
730 private static final String ANNUAL_POLE = IERS_BASE + "2003/annual.pole";
731
732
733 public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
734 throws OrekitException {
735 return new FundamentalNutationArguments(this, timeScale,
736 getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
737 }
738
739
740 @Override
741 public TimeScalarFunction getMeanObliquityFunction() throws OrekitException {
742
743
744 final PolynomialNutation epsilonA =
745 new PolynomialNutation(84381.448 * Constants.ARC_SECONDS_TO_RADIANS,
746 -46.84024 * Constants.ARC_SECONDS_TO_RADIANS,
747 -0.00059 * Constants.ARC_SECONDS_TO_RADIANS,
748 0.001813 * Constants.ARC_SECONDS_TO_RADIANS);
749
750 return new TimeScalarFunction() {
751
752
753 @Override
754 public double value(final AbsoluteDate date) {
755 return epsilonA.value(evaluateTC(date));
756 }
757
758
759 @Override
760 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
761 return epsilonA.value(evaluateTC(date));
762 }
763
764 };
765
766 }
767
768
769 @Override
770 public TimeVectorFunction getXYSpXY2Function()
771 throws OrekitException {
772
773
774 final FundamentalNutationArguments arguments = getNutationArguments(null);
775
776
777 final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
778 final PoissonSeriesParser parser =
779 new PoissonSeriesParser(17).
780 withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
781 withFirstDelaunay(4).
782 withFirstPlanetary(9).
783 withSinCos(0, 2, microAS, 3, microAS);
784
785 final PoissonSeries xSeries = parser.parse(getStream(X_SERIES), X_SERIES);
786 final PoissonSeries ySeries = parser.parse(getStream(Y_SERIES), Y_SERIES);
787 final PoissonSeries sSeries = parser.parse(getStream(S_SERIES), S_SERIES);
788 final PoissonSeries.CompiledSeries xys = PoissonSeries.compile(xSeries, ySeries, sSeries);
789
790
791 return new TimeVectorFunction() {
792
793
794 @Override
795 public double[] value(final AbsoluteDate date) {
796 return xys.value(arguments.evaluateAll(date));
797 }
798
799
800 @Override
801 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
802 return xys.value(arguments.evaluateAll(date));
803 }
804
805 };
806
807 }
808
809
810
811 @Override
812 public TimeVectorFunction getPrecessionFunction() throws OrekitException {
813
814
815
816 final PolynomialNutation psiA =
817 new PolynomialNutation( 0.0,
818 5038.47875 * Constants.ARC_SECONDS_TO_RADIANS,
819 -1.07259 * Constants.ARC_SECONDS_TO_RADIANS,
820 -0.001147 * Constants.ARC_SECONDS_TO_RADIANS);
821 final PolynomialNutation omegaA =
822 new PolynomialNutation(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
823 -0.02524 * Constants.ARC_SECONDS_TO_RADIANS,
824 0.05127 * Constants.ARC_SECONDS_TO_RADIANS,
825 -0.007726 * Constants.ARC_SECONDS_TO_RADIANS);
826 final PolynomialNutation chiA =
827 new PolynomialNutation( 0.0,
828 10.5526 * Constants.ARC_SECONDS_TO_RADIANS,
829 -2.38064 * Constants.ARC_SECONDS_TO_RADIANS,
830 -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);
831
832 return new PrecessionFunction(psiA, omegaA, chiA);
833
834 }
835
836
837 @Override
838 public TimeVectorFunction getNutationFunction()
839 throws OrekitException {
840
841
842 final FundamentalNutationArguments arguments = getNutationArguments(null);
843
844
845 final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
846 final PoissonSeriesParser luniSolarParser =
847 new PoissonSeriesParser(14).withFirstDelaunay(1);
848 final PoissonSeriesParser luniSolarPsiParser =
849 luniSolarParser.
850 withSinCos(0, 7, milliAS, 11, milliAS).
851 withSinCos(1, 8, milliAS, 12, milliAS);
852 final PoissonSeries psiLuniSolarSeries =
853 luniSolarPsiParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);
854 final PoissonSeriesParser luniSolarEpsilonParser =
855 luniSolarParser.
856 withSinCos(0, 13, milliAS, 9, milliAS).
857 withSinCos(1, 14, milliAS, 10, milliAS);
858 final PoissonSeries epsilonLuniSolarSeries =
859 luniSolarEpsilonParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);
860
861 final PoissonSeriesParser planetaryParser =
862 new PoissonSeriesParser(21).
863 withFirstDelaunay(2).
864 withFirstPlanetary(7);
865 final PoissonSeriesParser planetaryPsiParser =
866 planetaryParser.withSinCos(0, 17, milliAS, 18, milliAS);
867 final PoissonSeries psiPlanetarySeries =
868 planetaryPsiParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);
869 final PoissonSeriesParser planetaryEpsilonParser =
870 planetaryParser.withSinCos(0, 19, milliAS, 20, milliAS);
871 final PoissonSeries epsilonPlanetarySeries =
872 planetaryEpsilonParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);
873
874 final PoissonSeries.CompiledSeries luniSolarSeries =
875 PoissonSeries.compile(psiLuniSolarSeries, epsilonLuniSolarSeries);
876 final PoissonSeries.CompiledSeries planetarySeries =
877 PoissonSeries.compile(psiPlanetarySeries, epsilonPlanetarySeries);
878
879 return new TimeVectorFunction() {
880
881
882 @Override
883 public double[] value(final AbsoluteDate date) {
884 final BodiesElements elements = arguments.evaluateAll(date);
885 final double[] luniSolar = luniSolarSeries.value(elements);
886 final double[] planetary = planetarySeries.value(elements);
887 return new double[] {
888 luniSolar[0] + planetary[0], luniSolar[1] + planetary[1],
889 IAU1994ResolutionC7.value(elements)
890 };
891 }
892
893
894 @Override
895 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
896 final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
897 final T[] luniSolar = luniSolarSeries.value(elements);
898 final T[] planetary = planetarySeries.value(elements);
899 final T[] result = MathArrays.buildArray(date.getField(), 3);
900 result[0] = luniSolar[0].add(planetary[0]);
901 result[1] = luniSolar[1].add(planetary[1]);
902 result[2] = IAU1994ResolutionC7.value(elements);
903 return result;
904 }
905
906 };
907
908 }
909
910
911 @Override
912 public TimeScalarFunction getGMSTFunction(final TimeScale ut1)
913 throws OrekitException {
914
915
916 final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);
917
918
919
920 final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
921 final PoissonSeriesParser parser =
922 new PoissonSeriesParser(17).
923 withFirstDelaunay(4).
924 withFirstPlanetary(9).
925 withSinCos(0, 2, microAS, 3, microAS).
926 withPolynomialPart('t', Unit.ARC_SECONDS);
927 final PolynomialNutation minusEO =
928 parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();
929
930
931 return new TimeScalarFunction() {
932
933
934 @Override
935 public double value(final AbsoluteDate date) {
936 return era.value(date) + minusEO.value(evaluateTC(date));
937 }
938
939
940 @Override
941 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
942 return era.value(date).add(minusEO.value(evaluateTC(date)));
943 }
944
945 };
946
947 }
948
949
950 @Override
951 public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1)
952 throws OrekitException {
953
954
955 final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);
956
957
958
959 final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
960 final PoissonSeriesParser parser =
961 new PoissonSeriesParser(17).
962 withFirstDelaunay(4).
963 withFirstPlanetary(9).
964 withSinCos(0, 2, microAS, 3, microAS).
965 withPolynomialPart('t', Unit.ARC_SECONDS);
966 final PolynomialNutation minusEO =
967 parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();
968
969
970 return new TimeScalarFunction() {
971
972
973 @Override
974 public double value(final AbsoluteDate date) {
975 return era.getRate() + minusEO.derivative(evaluateTC(date));
976 }
977
978
979 @Override
980 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
981 return minusEO.derivative(evaluateTC(date)).add(era.getRate());
982 }
983
984 };
985
986 }
987
988
989 @Override
990 public TimeScalarFunction getGASTFunction(final TimeScale ut1, final EOPHistory eopHistory)
991 throws OrekitException {
992
993
994 final FundamentalNutationArguments arguments = getNutationArguments(null);
995
996
997 final TimeScalarFunction epsilon = getMeanObliquityFunction();
998
999
1000 final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
1001 final PoissonSeriesParser luniSolarPsiParser =
1002 new PoissonSeriesParser(14).
1003 withFirstDelaunay(1).
1004 withSinCos(0, 7, milliAS, 11, milliAS).
1005 withSinCos(1, 8, milliAS, 12, milliAS);
1006 final PoissonSeries psiLuniSolarSeries =
1007 luniSolarPsiParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);
1008
1009 final PoissonSeriesParser planetaryPsiParser =
1010 new PoissonSeriesParser(21).
1011 withFirstDelaunay(2).
1012 withFirstPlanetary(7).
1013 withSinCos(0, 17, milliAS, 18, milliAS);
1014 final PoissonSeries psiPlanetarySeries =
1015 planetaryPsiParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);
1016
1017 final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1018 final PoissonSeriesParser gstParser =
1019 new PoissonSeriesParser(17).
1020 withFirstDelaunay(4).
1021 withFirstPlanetary(9).
1022 withSinCos(0, 2, microAS, 3, microAS).
1023 withPolynomialPart('t', Unit.ARC_SECONDS);
1024 final PoissonSeries gstSeries = gstParser.parse(getStream(GST_SERIES), GST_SERIES);
1025 final PoissonSeries.CompiledSeries psiGstSeries =
1026 PoissonSeries.compile(psiLuniSolarSeries, psiPlanetarySeries, gstSeries);
1027
1028
1029 final TimeScalarFunction era = getEarthOrientationAngleFunction(ut1);
1030
1031 return new TimeScalarFunction() {
1032
1033
1034 @Override
1035 public double value(final AbsoluteDate date) {
1036
1037
1038 final BodiesElements elements = arguments.evaluateAll(date);
1039 final double[] angles = psiGstSeries.value(elements);
1040 final double ddPsi = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
1041 final double deltaPsi = angles[0] + angles[1] + ddPsi;
1042 final double epsilonA = epsilon.value(date);
1043
1044
1045
1046 return era.value(date) + deltaPsi * FastMath.cos(epsilonA) + angles[2];
1047
1048 }
1049
1050
1051 @Override
1052 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
1053
1054
1055 final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
1056 final T[] angles = psiGstSeries.value(elements);
1057 final T ddPsi = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
1058 final T deltaPsi = angles[0].add(angles[1]).add(ddPsi);
1059 final T epsilonA = epsilon.value(date);
1060
1061
1062
1063 return era.value(date).add(deltaPsi.multiply(epsilonA.cos())).add(angles[2]);
1064
1065 }
1066
1067 };
1068
1069 }
1070
1071
1072 @Override
1073 public TimeVectorFunction getEOPTidalCorrection()
1074 throws OrekitException {
1075
1076
1077
1078
1079
1080
1081 final FundamentalNutationArguments arguments = getNutationArguments(TimeScalesFactory.getTT());
1082
1083
1084 final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1085 final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
1086 withOptionalColumn(1).
1087 withGamma(2).
1088 withFirstDelaunay(3);
1089 final PoissonSeries xSeries =
1090 xyParser.
1091 withSinCos(0, 10, microAS, 11, microAS).
1092 parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
1093 final PoissonSeries ySeries =
1094 xyParser.
1095 withSinCos(0, 12, microAS, 13, microAS).
1096 parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
1097
1098 final double microS = 1.0e-6;
1099 final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
1100 withOptionalColumn(1).
1101 withGamma(2).
1102 withFirstDelaunay(3).
1103 withSinCos(0, 10, microS, 11, microS);
1104 final PoissonSeries ut1Series =
1105 ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);
1106
1107 return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
1108
1109 }
1110
1111
1112 public LoveNumbers getLoveNumbers() throws OrekitException {
1113 return loadLoveNumbers(LOVE_NUMBERS);
1114 }
1115
1116
1117 public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1)
1118 throws OrekitException {
1119
1120
1121 final FundamentalNutationArguments arguments = getNutationArguments(ut1);
1122
1123
1124 final PoissonSeriesParser k20Parser =
1125 new PoissonSeriesParser(18).
1126 withOptionalColumn(1).
1127 withDoodson(4, 2).
1128 withFirstDelaunay(10);
1129 final PoissonSeriesParser k21Parser =
1130 new PoissonSeriesParser(18).
1131 withOptionalColumn(1).
1132 withDoodson(4, 3).
1133 withFirstDelaunay(10);
1134 final PoissonSeriesParser k22Parser =
1135 new PoissonSeriesParser(16).
1136 withOptionalColumn(1).
1137 withDoodson(4, 2).
1138 withFirstDelaunay(10);
1139
1140 final double pico = 1.0e-12;
1141 final PoissonSeries c20Series =
1142 k20Parser.
1143 withSinCos(0, 18, -pico, 16, pico).
1144 parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
1145 final PoissonSeries c21Series =
1146 k21Parser.
1147 withSinCos(0, 17, pico, 18, pico).
1148 parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
1149 final PoissonSeries s21Series =
1150 k21Parser.
1151 withSinCos(0, 18, -pico, 17, pico).
1152 parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
1153 final PoissonSeries c22Series =
1154 k22Parser.
1155 withSinCos(0, -1, pico, 16, pico).
1156 parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
1157 final PoissonSeries s22Series =
1158 k22Parser.
1159 withSinCos(0, 16, -pico, -1, pico).
1160 parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
1161
1162 return new TideFrequencyDependenceFunction(arguments,
1163 c20Series,
1164 c21Series, s21Series,
1165 c22Series, s22Series);
1166
1167 }
1168
1169
1170 @Override
1171 public double getPermanentTide() throws OrekitException {
1172 return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
1173 }
1174
1175
1176 @Override
1177 public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory)
1178 throws OrekitException {
1179
1180
1181 final TimeScale utc = TimeScalesFactory.getUTC();
1182 final SimpleTimeStampedTableParser.RowConverter<MeanPole> converter =
1183 new SimpleTimeStampedTableParser.RowConverter<MeanPole>() {
1184
1185 @Override
1186 public MeanPole convert(final double[] rawFields) throws OrekitException {
1187 return new MeanPole(new AbsoluteDate((int) rawFields[0], 1, 1, utc),
1188 rawFields[1] * Constants.ARC_SECONDS_TO_RADIANS,
1189 rawFields[2] * Constants.ARC_SECONDS_TO_RADIANS);
1190 }
1191 };
1192 final SimpleTimeStampedTableParser<MeanPole> parser =
1193 new SimpleTimeStampedTableParser<MeanPole>(3, converter);
1194 final List<MeanPole> annualPoleList = parser.parse(getStream(ANNUAL_POLE), ANNUAL_POLE);
1195 final AbsoluteDate firstAnnualPoleDate = annualPoleList.get(0).getDate();
1196 final AbsoluteDate lastAnnualPoleDate = annualPoleList.get(annualPoleList.size() - 1).getDate();
1197 final ImmutableTimeStampedCache<MeanPole> annualCache =
1198 new ImmutableTimeStampedCache<MeanPole>(2, annualPoleList);
1199
1200
1201 final double xp0 = 0.054 * Constants.ARC_SECONDS_TO_RADIANS;
1202 final double xp0Dot = 0.00083 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
1203 final double yp0 = 0.357 * Constants.ARC_SECONDS_TO_RADIANS;
1204 final double yp0Dot = 0.00395 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
1205
1206
1207 final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
1208 final double ratio = 0.00115;
1209
1210 return new TimeVectorFunction() {
1211
1212
1213 @Override
1214 public double[] value(final AbsoluteDate date) {
1215
1216
1217 if (date.compareTo(firstAnnualPoleDate) <= 0) {
1218 return new double[] {
1219 0.0, 0.0
1220 };
1221 }
1222
1223
1224 double meanPoleX = 0;
1225 double meanPoleY = 0;
1226 if (date.compareTo(lastAnnualPoleDate) <= 0) {
1227
1228
1229 try {
1230 final HermiteInterpolator interpolator = new HermiteInterpolator();
1231 annualCache.getNeighbors(date).forEach(neighbor ->
1232 interpolator.addSamplePoint(neighbor.getDate().durationFrom(date),
1233 new double[] {
1234 neighbor.getX(), neighbor.getY()
1235 }));
1236 final double[] interpolated = interpolator.value(0);
1237 meanPoleX = interpolated[0];
1238 meanPoleY = interpolated[1];
1239 } catch (TimeStampedCacheException tsce) {
1240
1241 throw new OrekitInternalError(tsce);
1242 }
1243 } else {
1244
1245
1246
1247 final double t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
1248 meanPoleX = xp0 + t * xp0Dot;
1249 meanPoleY = yp0 + t * yp0Dot;
1250
1251 }
1252
1253
1254 final PoleCorrection correction = eopHistory.getPoleCorrection(date);
1255 final double m1 = correction.getXp() - meanPoleX;
1256 final double m2 = meanPoleY - correction.getYp();
1257
1258 return new double[] {
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285 globalFactor * (m1 - ratio * m2),
1286 globalFactor * (m2 + ratio * m1),
1287 };
1288
1289 }
1290
1291
1292 @Override
1293 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
1294
1295 final AbsoluteDate aDate = date.toAbsoluteDate();
1296
1297
1298 if (aDate.compareTo(firstAnnualPoleDate) <= 0) {
1299 return MathArrays.buildArray(date.getField(), 2);
1300 }
1301
1302
1303 T meanPoleX = date.getField().getZero();
1304 T meanPoleY = date.getField().getZero();
1305 if (aDate.compareTo(lastAnnualPoleDate) <= 0) {
1306
1307
1308 try {
1309 final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
1310 final T[] y = MathArrays.buildArray(date.getField(), 2);
1311 final T zero = date.getField().getZero();
1312 final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero);
1313
1314
1315 annualCache.getNeighbors(aDate).forEach(neighbor -> {
1316 y[0] = zero.add(neighbor.getX());
1317 y[1] = zero.add(neighbor.getY());
1318 interpolator.addSamplePoint(central.durationFrom(neighbor.getDate()).negate(), y);
1319 });
1320 final T[] interpolated = interpolator.value(date.durationFrom(central));
1321 meanPoleX = interpolated[0];
1322 meanPoleY = interpolated[1];
1323 } catch (TimeStampedCacheException tsce) {
1324
1325 throw new OrekitInternalError(tsce);
1326 }
1327 } else {
1328
1329
1330
1331 final T t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
1332 meanPoleX = t.multiply(xp0Dot).add(xp0);
1333 meanPoleY = t.multiply(yp0Dot).add(yp0);
1334
1335 }
1336
1337
1338 final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
1339 final T m1 = correction.getXp().subtract(meanPoleX);
1340 final T m2 = meanPoleY.subtract(correction.getYp());
1341
1342 final T[] a = MathArrays.buildArray(date.getField(), 2);
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370 a[0] = m1.add(m2.multiply(-ratio)).multiply(globalFactor);
1371 a[1] = m2.add(m1.multiply( ratio)).multiply(globalFactor);
1372
1373 return a;
1374
1375 }
1376
1377 };
1378
1379 }
1380
1381
1382 @Override
1383 public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory)
1384 throws OrekitException {
1385
1386 return new TimeVectorFunction() {
1387
1388
1389 @Override
1390 public double[] value(final AbsoluteDate date) {
1391
1392 return new double[] {
1393 0.0, 0.0
1394 };
1395 }
1396
1397
1398 @Override
1399 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
1400
1401 return MathArrays.buildArray(date.getField(), 2);
1402 }
1403
1404 };
1405 }
1406
1407 },
1408
1409
1410 IERS_2010 {
1411
1412
1413 private static final String NUTATION_ARGUMENTS = IERS_BASE + "2010/nutation-arguments.txt";
1414
1415
1416 private static final String X_SERIES = IERS_BASE + "2010/tab5.2a.txt";
1417
1418
1419 private static final String Y_SERIES = IERS_BASE + "2010/tab5.2b.txt";
1420
1421
1422 private static final String S_SERIES = IERS_BASE + "2010/tab5.2d.txt";
1423
1424
1425 private static final String PSI_SERIES = IERS_BASE + "2010/tab5.3a.txt";
1426
1427
1428 private static final String EPSILON_SERIES = IERS_BASE + "2010/tab5.3b.txt";
1429
1430
1431 private static final String GST_SERIES = IERS_BASE + "2010/tab5.2e.txt";
1432
1433
1434 private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "2010/tab8.2ab.txt";
1435
1436
1437 private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "2010/tab8.3ab.txt";
1438
1439
1440 private static final String LOVE_NUMBERS = IERS_BASE + "2010/tab6.3.txt";
1441
1442
1443 private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5b.txt";
1444
1445
1446 private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5a.txt";
1447
1448
1449 private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5c.txt";
1450
1451
1452 public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
1453 throws OrekitException {
1454 return new FundamentalNutationArguments(this, timeScale,
1455 getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
1456 }
1457
1458
1459 @Override
1460 public TimeScalarFunction getMeanObliquityFunction() throws OrekitException {
1461
1462
1463 final PolynomialNutation epsilonA =
1464 new PolynomialNutation(84381.406 * Constants.ARC_SECONDS_TO_RADIANS,
1465 -46.836769 * Constants.ARC_SECONDS_TO_RADIANS,
1466 -0.0001831 * Constants.ARC_SECONDS_TO_RADIANS,
1467 0.00200340 * Constants.ARC_SECONDS_TO_RADIANS,
1468 -0.000000576 * Constants.ARC_SECONDS_TO_RADIANS,
1469 -0.0000000434 * Constants.ARC_SECONDS_TO_RADIANS);
1470
1471 return new TimeScalarFunction() {
1472
1473
1474 @Override
1475 public double value(final AbsoluteDate date) {
1476 return epsilonA.value(evaluateTC(date));
1477 }
1478
1479
1480 @Override
1481 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
1482 return epsilonA.value(evaluateTC(date));
1483 }
1484
1485 };
1486
1487 }
1488
1489
1490 @Override
1491 public TimeVectorFunction getXYSpXY2Function() throws OrekitException {
1492
1493
1494 final FundamentalNutationArguments arguments = getNutationArguments(null);
1495
1496
1497 final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1498 final PoissonSeriesParser parser =
1499 new PoissonSeriesParser(17).
1500 withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
1501 withFirstDelaunay(4).
1502 withFirstPlanetary(9).
1503 withSinCos(0, 2, microAS, 3, microAS);
1504 final PoissonSeries xSeries = parser.parse(getStream(X_SERIES), X_SERIES);
1505 final PoissonSeries ySeries = parser.parse(getStream(Y_SERIES), Y_SERIES);
1506 final PoissonSeries sSeries = parser.parse(getStream(S_SERIES), S_SERIES);
1507 final PoissonSeries.CompiledSeries xys = PoissonSeries.compile(xSeries, ySeries, sSeries);
1508
1509
1510 return new TimeVectorFunction() {
1511
1512
1513 @Override
1514 public double[] value(final AbsoluteDate date) {
1515 return xys.value(arguments.evaluateAll(date));
1516 }
1517
1518
1519 @Override
1520 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
1521 return xys.value(arguments.evaluateAll(date));
1522 }
1523
1524 };
1525
1526 }
1527
1528
1529 public LoveNumbers getLoveNumbers() throws OrekitException {
1530 return loadLoveNumbers(LOVE_NUMBERS);
1531 }
1532
1533
1534 public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1)
1535 throws OrekitException {
1536
1537
1538 final FundamentalNutationArguments arguments = getNutationArguments(ut1);
1539
1540
1541 final PoissonSeriesParser k20Parser =
1542 new PoissonSeriesParser(18).
1543 withOptionalColumn(1).
1544 withDoodson(4, 2).
1545 withFirstDelaunay(10);
1546 final PoissonSeriesParser k21Parser =
1547 new PoissonSeriesParser(18).
1548 withOptionalColumn(1).
1549 withDoodson(4, 3).
1550 withFirstDelaunay(10);
1551 final PoissonSeriesParser k22Parser =
1552 new PoissonSeriesParser(16).
1553 withOptionalColumn(1).
1554 withDoodson(4, 2).
1555 withFirstDelaunay(10);
1556
1557 final double pico = 1.0e-12;
1558 final PoissonSeries c20Series =
1559 k20Parser.
1560 withSinCos(0, 18, -pico, 16, pico).
1561 parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
1562 final PoissonSeries c21Series =
1563 k21Parser.
1564 withSinCos(0, 17, pico, 18, pico).
1565 parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
1566 final PoissonSeries s21Series =
1567 k21Parser.
1568 withSinCos(0, 18, -pico, 17, pico).
1569 parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
1570 final PoissonSeries c22Series =
1571 k22Parser.
1572 withSinCos(0, -1, pico, 16, pico).
1573 parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
1574 final PoissonSeries s22Series =
1575 k22Parser.
1576 withSinCos(0, 16, -pico, -1, pico).
1577 parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
1578
1579 return new TideFrequencyDependenceFunction(arguments,
1580 c20Series,
1581 c21Series, s21Series,
1582 c22Series, s22Series);
1583
1584 }
1585
1586
1587 @Override
1588 public double getPermanentTide() throws OrekitException {
1589 return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
1590 }
1591
1592
1593
1594
1595
1596
1597 private double[] computePoleWobble(final AbsoluteDate date, final EOPHistory eopHistory) {
1598
1599
1600 final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
1601 final double f1 = f0 / Constants.JULIAN_YEAR;
1602 final double f2 = f1 / Constants.JULIAN_YEAR;
1603 final double f3 = f2 / Constants.JULIAN_YEAR;
1604 final AbsoluteDate changeDate = new AbsoluteDate(2010, 1, 1, TimeScalesFactory.getTT());
1605
1606
1607 final double[] xPolynomial;
1608 final double[] yPolynomial;
1609 if (date.compareTo(changeDate) <= 0) {
1610 xPolynomial = new double[] {
1611 55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
1612 };
1613 yPolynomial = new double[] {
1614 346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
1615 };
1616 } else {
1617 xPolynomial = new double[] {
1618 23.513 * f0, 7.6141 * f1
1619 };
1620 yPolynomial = new double[] {
1621 358.891 * f0, -0.6287 * f1
1622 };
1623 }
1624 double meanPoleX = 0;
1625 double meanPoleY = 0;
1626 final double t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
1627 for (int i = xPolynomial.length - 1; i >= 0; --i) {
1628 meanPoleX = meanPoleX * t + xPolynomial[i];
1629 }
1630 for (int i = yPolynomial.length - 1; i >= 0; --i) {
1631 meanPoleY = meanPoleY * t + yPolynomial[i];
1632 }
1633
1634
1635 final PoleCorrection correction = eopHistory.getPoleCorrection(date);
1636 final double m1 = correction.getXp() - meanPoleX;
1637 final double m2 = meanPoleY - correction.getYp();
1638
1639 return new double[] {
1640 m1, m2
1641 };
1642
1643 }
1644
1645
1646
1647
1648
1649
1650
1651 private <T extends RealFieldElement<T>> T[] computePoleWobble(final FieldAbsoluteDate<T> date, final EOPHistory eopHistory) {
1652
1653
1654 final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
1655 final double f1 = f0 / Constants.JULIAN_YEAR;
1656 final double f2 = f1 / Constants.JULIAN_YEAR;
1657 final double f3 = f2 / Constants.JULIAN_YEAR;
1658 final AbsoluteDate changeDate = new AbsoluteDate(2010, 1, 1, TimeScalesFactory.getTT());
1659
1660
1661 final double[] xPolynomial;
1662 final double[] yPolynomial;
1663 if (date.toAbsoluteDate().compareTo(changeDate) <= 0) {
1664 xPolynomial = new double[] {
1665 55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
1666 };
1667 yPolynomial = new double[] {
1668 346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
1669 };
1670 } else {
1671 xPolynomial = new double[] {
1672 23.513 * f0, 7.6141 * f1
1673 };
1674 yPolynomial = new double[] {
1675 358.891 * f0, -0.6287 * f1
1676 };
1677 }
1678 T meanPoleX = date.getField().getZero();
1679 T meanPoleY = date.getField().getZero();
1680 final T t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
1681 for (int i = xPolynomial.length - 1; i >= 0; --i) {
1682 meanPoleX = meanPoleX.multiply(t).add(xPolynomial[i]);
1683 }
1684 for (int i = yPolynomial.length - 1; i >= 0; --i) {
1685 meanPoleY = meanPoleY.multiply(t).add(yPolynomial[i]);
1686 }
1687
1688
1689 final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
1690 final T[] m = MathArrays.buildArray(date.getField(), 2);
1691 m[0] = correction.getXp().subtract(meanPoleX);
1692 m[1] = meanPoleY.subtract(correction.getYp());
1693
1694 return m;
1695
1696 }
1697
1698
1699 @Override
1700 public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory)
1701 throws OrekitException {
1702
1703
1704 final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
1705 final double ratio = 0.00115;
1706
1707 return new TimeVectorFunction() {
1708
1709
1710 @Override
1711 public double[] value(final AbsoluteDate date) {
1712
1713
1714 final double[] wobbleM = computePoleWobble(date, eopHistory);
1715
1716 return new double[] {
1717
1718
1719
1720
1721
1722
1723
1724
1725 globalFactor * (wobbleM[0] + ratio * wobbleM[1]),
1726 globalFactor * (wobbleM[1] - ratio * wobbleM[0])
1727 };
1728
1729 }
1730
1731
1732 @Override
1733 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
1734
1735
1736 final T[] wobbleM = computePoleWobble(date, eopHistory);
1737
1738 final T[] a = MathArrays.buildArray(date.getField(), 2);
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748 a[0] = wobbleM[0].add(wobbleM[1].multiply( ratio)).multiply(globalFactor);
1749 a[1] = wobbleM[1].add(wobbleM[0].multiply(-ratio)).multiply(globalFactor);
1750
1751 return a;
1752
1753 }
1754
1755 };
1756
1757 }
1758
1759
1760 @Override
1761 public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory)
1762 throws OrekitException {
1763
1764 return new TimeVectorFunction() {
1765
1766
1767 @Override
1768 public double[] value(final AbsoluteDate date) {
1769
1770
1771 final double[] wobbleM = computePoleWobble(date, eopHistory);
1772
1773 return new double[] {
1774
1775
1776
1777
1778 -2.1778e-10 * (wobbleM[0] - 0.01724 * wobbleM[1]) / Constants.ARC_SECONDS_TO_RADIANS,
1779 -1.7232e-10 * (wobbleM[1] - 0.03365 * wobbleM[0]) / Constants.ARC_SECONDS_TO_RADIANS
1780 };
1781
1782 }
1783
1784
1785 @Override
1786 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
1787
1788 final T[] v = MathArrays.buildArray(date.getField(), 2);
1789
1790
1791 final T[] wobbleM = computePoleWobble(date, eopHistory);
1792
1793
1794
1795
1796
1797 v[0] = wobbleM[0].subtract(wobbleM[1].multiply(0.01724)).multiply(-2.1778e-10 / Constants.ARC_SECONDS_TO_RADIANS);
1798 v[1] = wobbleM[1].subtract(wobbleM[0].multiply(0.03365)).multiply(-1.7232e-10 / Constants.ARC_SECONDS_TO_RADIANS);
1799
1800 return v;
1801
1802 }
1803
1804 };
1805
1806 }
1807
1808
1809 @Override
1810 public TimeVectorFunction getPrecessionFunction() throws OrekitException {
1811
1812
1813
1814 final PolynomialNutation psiA =
1815 new PolynomialNutation( 0.0,
1816 5038.481507 * Constants.ARC_SECONDS_TO_RADIANS,
1817 -1.0790069 * Constants.ARC_SECONDS_TO_RADIANS,
1818 -0.00114045 * Constants.ARC_SECONDS_TO_RADIANS,
1819 0.000132851 * Constants.ARC_SECONDS_TO_RADIANS,
1820 -0.0000000951 * Constants.ARC_SECONDS_TO_RADIANS);
1821 final PolynomialNutation omegaA =
1822 new PolynomialNutation(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
1823 -0.025754 * Constants.ARC_SECONDS_TO_RADIANS,
1824 0.0512623 * Constants.ARC_SECONDS_TO_RADIANS,
1825 -0.00772503 * Constants.ARC_SECONDS_TO_RADIANS,
1826 -0.000000467 * Constants.ARC_SECONDS_TO_RADIANS,
1827 0.0000003337 * Constants.ARC_SECONDS_TO_RADIANS);
1828 final PolynomialNutation chiA =
1829 new PolynomialNutation( 0.0,
1830 10.556403 * Constants.ARC_SECONDS_TO_RADIANS,
1831 -2.3814292 * Constants.ARC_SECONDS_TO_RADIANS,
1832 -0.00121197 * Constants.ARC_SECONDS_TO_RADIANS,
1833 0.000170663 * Constants.ARC_SECONDS_TO_RADIANS,
1834 -0.0000000560 * Constants.ARC_SECONDS_TO_RADIANS);
1835
1836 return new PrecessionFunction(psiA, omegaA, chiA);
1837
1838 }
1839
1840
1841 @Override
1842 public TimeVectorFunction getNutationFunction()
1843 throws OrekitException {
1844
1845
1846 final FundamentalNutationArguments arguments = getNutationArguments(null);
1847
1848
1849 final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1850 final PoissonSeriesParser parser =
1851 new PoissonSeriesParser(17).
1852 withFirstDelaunay(4).
1853 withFirstPlanetary(9).
1854 withSinCos(0, 2, microAS, 3, microAS);
1855 final PoissonSeries psiSeries = parser.parse(getStream(PSI_SERIES), PSI_SERIES);
1856 final PoissonSeries epsilonSeries = parser.parse(getStream(EPSILON_SERIES), EPSILON_SERIES);
1857 final PoissonSeries.CompiledSeries psiEpsilonSeries =
1858 PoissonSeries.compile(psiSeries, epsilonSeries);
1859
1860 return new TimeVectorFunction() {
1861
1862
1863 @Override
1864 public double[] value(final AbsoluteDate date) {
1865 final BodiesElements elements = arguments.evaluateAll(date);
1866 final double[] psiEpsilon = psiEpsilonSeries.value(elements);
1867 return new double[] {
1868 psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements)
1869 };
1870 }
1871
1872
1873 @Override
1874 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
1875 final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
1876 final T[] psiEpsilon = psiEpsilonSeries.value(elements);
1877 final T[] result = MathArrays.buildArray(date.getField(), 3);
1878 result[0] = psiEpsilon[0];
1879 result[1] = psiEpsilon[1];
1880 result[2] = IAU1994ResolutionC7.value(elements);
1881 return result;
1882 }
1883
1884 };
1885
1886 }
1887
1888
1889 @Override
1890 public TimeScalarFunction getGMSTFunction(final TimeScale ut1) throws OrekitException {
1891
1892
1893 final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);
1894
1895
1896
1897 final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1898 final PoissonSeriesParser parser =
1899 new PoissonSeriesParser(17).
1900 withFirstDelaunay(4).
1901 withFirstPlanetary(9).
1902 withSinCos(0, 2, microAS, 3, microAS).
1903 withPolynomialPart('t', Unit.ARC_SECONDS);
1904 final PolynomialNutation minusEO =
1905 parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();
1906
1907
1908 return new TimeScalarFunction() {
1909
1910
1911 @Override
1912 public double value(final AbsoluteDate date) {
1913 return era.value(date) + minusEO.value(evaluateTC(date));
1914 }
1915
1916
1917 @Override
1918 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
1919 return era.value(date).add(minusEO.value(evaluateTC(date)));
1920 }
1921
1922 };
1923
1924 }
1925
1926
1927 @Override
1928 public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1) throws OrekitException {
1929
1930
1931 final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);
1932
1933
1934
1935 final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1936 final PoissonSeriesParser parser =
1937 new PoissonSeriesParser(17).
1938 withFirstDelaunay(4).
1939 withFirstPlanetary(9).
1940 withSinCos(0, 2, microAS, 3, microAS).
1941 withPolynomialPart('t', Unit.ARC_SECONDS);
1942 final PolynomialNutation minusEO =
1943 parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();
1944
1945
1946 return new TimeScalarFunction() {
1947
1948
1949 @Override
1950 public double value(final AbsoluteDate date) {
1951 return era.getRate() + minusEO.derivative(evaluateTC(date));
1952 }
1953
1954
1955 @Override
1956 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
1957 return minusEO.derivative(evaluateTC(date)).add(era.getRate());
1958 }
1959
1960 };
1961
1962 }
1963
1964
1965 @Override
1966 public TimeScalarFunction getGASTFunction(final TimeScale ut1, final EOPHistory eopHistory)
1967 throws OrekitException {
1968
1969
1970 final FundamentalNutationArguments arguments = getNutationArguments(null);
1971
1972
1973 final TimeScalarFunction epsilon = getMeanObliquityFunction();
1974
1975
1976 final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1977 final PoissonSeriesParser baseParser =
1978 new PoissonSeriesParser(17).
1979 withFirstDelaunay(4).
1980 withFirstPlanetary(9).
1981 withSinCos(0, 2, microAS, 3, microAS);
1982 final PoissonSeriesParser gstParser = baseParser.withPolynomialPart('t', Unit.ARC_SECONDS);
1983 final PoissonSeries psiSeries = baseParser.parse(getStream(PSI_SERIES), PSI_SERIES);
1984 final PoissonSeries gstSeries = gstParser.parse(getStream(GST_SERIES), GST_SERIES);
1985 final PoissonSeries.CompiledSeries psiGstSeries =
1986 PoissonSeries.compile(psiSeries, gstSeries);
1987
1988
1989 final TimeScalarFunction era = getEarthOrientationAngleFunction(ut1);
1990
1991 return new TimeScalarFunction() {
1992
1993
1994 @Override
1995 public double value(final AbsoluteDate date) {
1996
1997
1998 final BodiesElements elements = arguments.evaluateAll(date);
1999 final double[] angles = psiGstSeries.value(elements);
2000 final double ddPsi = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
2001 final double deltaPsi = angles[0] + ddPsi;
2002 final double epsilonA = epsilon.value(date);
2003
2004
2005
2006 return era.value(date) + deltaPsi * FastMath.cos(epsilonA) + angles[1];
2007
2008 }
2009
2010
2011 @Override
2012 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
2013
2014
2015 final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
2016 final T[] angles = psiGstSeries.value(elements);
2017 final T ddPsi = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
2018 final T deltaPsi = angles[0].add(ddPsi);
2019 final T epsilonA = epsilon.value(date);
2020
2021
2022
2023 return era.value(date).add(deltaPsi.multiply(epsilonA.cos())).add(angles[1]);
2024
2025 }
2026
2027 };
2028
2029 }
2030
2031
2032 @Override
2033 public TimeVectorFunction getEOPTidalCorrection()
2034 throws OrekitException {
2035
2036
2037
2038
2039
2040
2041 final FundamentalNutationArguments arguments = getNutationArguments(TimeScalesFactory.getTT());
2042
2043
2044 final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
2045 final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
2046 withOptionalColumn(1).
2047 withGamma(2).
2048 withFirstDelaunay(3);
2049 final PoissonSeries xSeries =
2050 xyParser.
2051 withSinCos(0, 10, microAS, 11, microAS).
2052 parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
2053 final PoissonSeries ySeries =
2054 xyParser.
2055 withSinCos(0, 12, microAS, 13, microAS).
2056 parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
2057
2058 final double microS = 1.0e-6;
2059 final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
2060 withOptionalColumn(1).
2061 withGamma(2).
2062 withFirstDelaunay(3).
2063 withSinCos(0, 10, microS, 11, microS);
2064 final PoissonSeries ut1Series =
2065 ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);
2066
2067 return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
2068
2069 }
2070
2071 };
2072
2073
2074 private static final String IERS_BASE = "/assets/org/orekit/IERS-conventions/";
2075
2076
2077
2078
2079
2080 public AbsoluteDate getNutationReferenceEpoch() {
2081
2082 return AbsoluteDate.J2000_EPOCH;
2083 }
2084
2085
2086
2087
2088
2089
2090 public double evaluateTC(final AbsoluteDate date) {
2091 return date.durationFrom(getNutationReferenceEpoch()) / Constants.JULIAN_CENTURY;
2092 }
2093
2094
2095
2096
2097
2098
2099
2100 public <T extends RealFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date) {
2101 return date.durationFrom(getNutationReferenceEpoch()).divide(Constants.JULIAN_CENTURY);
2102 }
2103
2104
2105
2106
2107
2108
2109
2110
2111 public abstract FundamentalNutationArguments getNutationArguments(TimeScale timeScale)
2112 throws OrekitException;
2113
2114
2115
2116
2117
2118
2119 public abstract TimeScalarFunction getMeanObliquityFunction() throws OrekitException;
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129 public abstract TimeVectorFunction getXYSpXY2Function()
2130 throws OrekitException;
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143 public TimeScalarFunction getEarthOrientationAngleFunction(final TimeScale ut1) {
2144 return new StellarAngleCapitaine(ut1);
2145 }
2146
2147
2148
2149
2150
2151
2152
2153
2154
2155
2156
2157
2158
2159
2160
2161 public abstract TimeVectorFunction getPrecessionFunction() throws OrekitException;
2162
2163
2164
2165
2166
2167
2168
2169
2170
2171
2172
2173
2174 public abstract TimeVectorFunction getNutationFunction()
2175 throws OrekitException;
2176
2177
2178
2179
2180
2181
2182
2183 public abstract TimeScalarFunction getGMSTFunction(TimeScale ut1)
2184 throws OrekitException;
2185
2186
2187
2188
2189
2190
2191
2192 public abstract TimeScalarFunction getGMSTRateFunction(TimeScale ut1)
2193 throws OrekitException;
2194
2195
2196
2197
2198
2199
2200
2201
2202 public abstract TimeScalarFunction getGASTFunction(TimeScale ut1, EOPHistory eopHistory)
2203 throws OrekitException;
2204
2205
2206
2207
2208
2209
2210
2211 public abstract TimeVectorFunction getEOPTidalCorrection()
2212 throws OrekitException;
2213
2214
2215
2216
2217
2218
2219 public abstract LoveNumbers getLoveNumbers()
2220 throws OrekitException;
2221
2222
2223
2224
2225
2226
2227
2228 public abstract TimeVectorFunction getTideFrequencyDependenceFunction(TimeScale ut1)
2229 throws OrekitException;
2230
2231
2232
2233
2234
2235 public abstract double getPermanentTide() throws OrekitException;
2236
2237
2238
2239
2240
2241
2242
2243 public abstract TimeVectorFunction getSolidPoleTide(EOPHistory eopHistory)
2244 throws OrekitException;
2245
2246
2247
2248
2249
2250
2251
2252 public abstract TimeVectorFunction getOceanPoleTide(EOPHistory eopHistory)
2253 throws OrekitException;
2254
2255
2256
2257
2258
2259
2260
2261
2262
2263 public interface NutationCorrectionConverter {
2264
2265
2266
2267
2268
2269
2270
2271
2272 double[] toNonRotating(AbsoluteDate date, double ddPsi, double ddEpsilon)
2273 throws OrekitException;
2274
2275
2276
2277
2278
2279
2280
2281
2282 double[] toEquinox(AbsoluteDate date, double dX, double dY)
2283 throws OrekitException;
2284
2285 }
2286
2287
2288
2289
2290
2291
2292
2293
2294
2295
2296
2297 public NutationCorrectionConverter getNutationCorrectionConverter()
2298 throws OrekitException {
2299
2300
2301 final TimeVectorFunction precessionFunction = getPrecessionFunction();
2302 final TimeScalarFunction epsilonAFunction = getMeanObliquityFunction();
2303 final AbsoluteDate date0 = getNutationReferenceEpoch();
2304 final double cosE0 = FastMath.cos(epsilonAFunction.value(date0));
2305
2306 return new NutationCorrectionConverter() {
2307
2308
2309 @Override
2310 public double[] toNonRotating(final AbsoluteDate date,
2311 final double ddPsi, final double ddEpsilon)
2312 throws OrekitException {
2313
2314 final double[] angles = precessionFunction.value(date);
2315
2316
2317 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
2318 final double c = angles[0] * cosE0 - angles[2];
2319
2320
2321 return new double[] {
2322 sinEA * ddPsi + c * ddEpsilon,
2323 ddEpsilon - c * sinEA * ddPsi
2324 };
2325
2326 }
2327
2328
2329 @Override
2330 public double[] toEquinox(final AbsoluteDate date,
2331 final double dX, final double dY)
2332 throws OrekitException {
2333
2334 final double[] angles = precessionFunction.value(date);
2335
2336
2337 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
2338 final double c = angles[0] * cosE0 - angles[2];
2339 final double opC2 = 1 + c * c;
2340
2341
2342 return new double[] {
2343 (dX - c * dY) / (sinEA * opC2),
2344 (dY + c * dX) / opC2
2345 };
2346 }
2347
2348 };
2349
2350 }
2351
2352
2353
2354
2355
2356
2357
2358 protected LoveNumbers loadLoveNumbers(final String nameLove) throws OrekitException {
2359 try {
2360
2361
2362 final double[][] real = new double[4][];
2363 final double[][] imaginary = new double[4][];
2364 final double[][] plus = new double[4][];
2365 for (int i = 0; i < real.length; ++i) {
2366 real[i] = new double[i + 1];
2367 imaginary[i] = new double[i + 1];
2368 plus[i] = new double[i + 1];
2369 }
2370
2371 try (InputStream stream = IERSConventions.class.getResourceAsStream(nameLove)) {
2372
2373 if (stream == null) {
2374
2375 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, nameLove);
2376 }
2377
2378
2379 try (BufferedReader reader = new BufferedReader(new InputStreamReader(stream, "UTF-8"))) {
2380
2381 String line = reader.readLine();
2382 int lineNumber = 1;
2383
2384
2385 while (line != null) {
2386
2387 line = line.trim();
2388 if (!(line.isEmpty() || line.startsWith("#"))) {
2389 final String[] fields = line.split("\\p{Space}+");
2390 if (fields.length != 5) {
2391
2392 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
2393 lineNumber, nameLove, line);
2394 }
2395 final int n = Integer.parseInt(fields[0]);
2396 final int m = Integer.parseInt(fields[1]);
2397 if ((n < 2) || (n > 3) || (m < 0) || (m > n)) {
2398
2399 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
2400 lineNumber, nameLove, line);
2401
2402 }
2403 real[n][m] = Double.parseDouble(fields[2]);
2404 imaginary[n][m] = Double.parseDouble(fields[3]);
2405 plus[n][m] = Double.parseDouble(fields[4]);
2406 }
2407
2408
2409 lineNumber++;
2410 line = reader.readLine();
2411
2412 }
2413 }
2414 }
2415
2416 return new LoveNumbers(real, imaginary, plus);
2417
2418 } catch (IOException ioe) {
2419
2420 throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, nameLove);
2421 }
2422 }
2423
2424
2425
2426
2427
2428 private static InputStream getStream(final String name) {
2429 return IERSConventions.class.getResourceAsStream(name);
2430 }
2431
2432
2433
2434
2435
2436
2437 private static class IAU1994ResolutionC7 {
2438
2439
2440 private static final double EQE1 = 0.00264 * Constants.ARC_SECONDS_TO_RADIANS;
2441
2442
2443 private static final double EQE2 = 0.000063 * Constants.ARC_SECONDS_TO_RADIANS;
2444
2445
2446
2447
2448 private static final AbsoluteDate MODEL_START =
2449 new AbsoluteDate(1997, 2, 27, 0, 0, 30, TimeScalesFactory.getTAI());
2450
2451
2452
2453
2454
2455 public static double value(final DelaunayArguments arguments) {
2456 if (arguments.getDate().compareTo(MODEL_START) >= 0) {
2457
2458
2459
2460
2461
2462 final double om = arguments.getOmega();
2463
2464
2465 return EQE1 * FastMath.sin(om) + EQE2 * FastMath.sin(om + om);
2466
2467 } else {
2468 return 0.0;
2469 }
2470 }
2471
2472
2473
2474
2475
2476
2477 public static <T extends RealFieldElement<T>> T value(final FieldDelaunayArguments<T> arguments) {
2478 if (arguments.getDate().toAbsoluteDate().compareTo(MODEL_START) >= 0) {
2479
2480
2481
2482
2483
2484 final T om = arguments.getOmega();
2485
2486
2487 return om.sin().multiply(EQE1).add(om.add(om).sin().multiply(EQE2));
2488
2489 } else {
2490 return arguments.getDate().getField().getZero();
2491 }
2492 }
2493
2494 };
2495
2496
2497
2498
2499
2500
2501
2502
2503
2504
2505
2506
2507
2508
2509
2510
2511 private static class StellarAngleCapitaine implements TimeScalarFunction {
2512
2513
2514 private static final AbsoluteDate REFERENCE_DATE = new AbsoluteDate(DateComponents.J2000_EPOCH,
2515 TimeComponents.H12,
2516 TimeScalesFactory.getTAI());
2517
2518
2519 private static final double ERA_0 = MathUtils.TWO_PI * 0.7790572732640;
2520
2521
2522
2523 private static final double ERA_1A = MathUtils.TWO_PI / Constants.JULIAN_DAY;
2524
2525
2526
2527 private static final double ERA_1B = ERA_1A * 0.00273781191135448;
2528
2529
2530 private final TimeScale ut1;
2531
2532
2533
2534
2535 StellarAngleCapitaine(final TimeScale ut1) {
2536 this.ut1 = ut1;
2537 }
2538
2539
2540
2541
2542 public double getRate() {
2543 return ERA_1A + ERA_1B;
2544 }
2545
2546
2547 @Override
2548 public double value(final AbsoluteDate date) {
2549
2550
2551 final int secondsInDay = 86400;
2552 final double dt = date.durationFrom(REFERENCE_DATE);
2553 final long days = ((long) dt) / secondsInDay;
2554 final double dtA = secondsInDay * days;
2555 final double dtB = (dt - dtA) + ut1.offsetFromTAI(date);
2556
2557 return ERA_0 + ERA_1A * dtB + ERA_1B * (dtA + dtB);
2558
2559 }
2560
2561
2562 @Override
2563 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
2564
2565
2566 final int secondsInDay = 86400;
2567 final T dt = date.durationFrom(REFERENCE_DATE);
2568 final long days = ((long) dt.getReal()) / secondsInDay;
2569 final double dtA = secondsInDay * days;
2570 final T dtB = dt.subtract(dtA).add(ut1.offsetFromTAI(date.toAbsoluteDate()));
2571
2572 return dtB.add(dtA).multiply(ERA_1B).add(dtB.multiply(ERA_1A)).add(ERA_0);
2573
2574 }
2575
2576 }
2577
2578
2579 private static class MeanPole implements TimeStamped, Serializable {
2580
2581
2582 private static final long serialVersionUID = 20131028l;
2583
2584
2585 private final AbsoluteDate date;
2586
2587
2588 private double x;
2589
2590
2591 private double y;
2592
2593
2594
2595
2596
2597
2598 MeanPole(final AbsoluteDate date, final double x, final double y) {
2599 this.date = date;
2600 this.x = x;
2601 this.y = y;
2602 }
2603
2604
2605 @Override
2606 public AbsoluteDate getDate() {
2607 return date;
2608 }
2609
2610
2611
2612
2613 public double getX() {
2614 return x;
2615 }
2616
2617
2618
2619
2620 public double getY() {
2621 return y;
2622 }
2623
2624 }
2625
2626
2627 private class PrecessionFunction implements TimeVectorFunction {
2628
2629
2630 private final PolynomialNutation psiA;
2631
2632
2633 private final PolynomialNutation omegaA;
2634
2635
2636 private final PolynomialNutation chiA;
2637
2638
2639
2640
2641
2642
2643 PrecessionFunction(final PolynomialNutation psiA,
2644 final PolynomialNutation omegaA,
2645 final PolynomialNutation chiA) {
2646 this.psiA = psiA;
2647 this.omegaA = omegaA;
2648 this.chiA = chiA;
2649 }
2650
2651
2652
2653 @Override
2654 public double[] value(final AbsoluteDate date) {
2655 final double tc = evaluateTC(date);
2656 return new double[] {
2657 psiA.value(tc), omegaA.value(tc), chiA.value(tc)
2658 };
2659 }
2660
2661
2662 @Override
2663 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
2664 final T[] a = MathArrays.buildArray(date.getField(), 3);
2665 final T tc = evaluateTC(date);
2666 a[0] = psiA.value(tc);
2667 a[1] = omegaA.value(tc);
2668 a[2] = chiA.value(tc);
2669 return a;
2670 }
2671
2672 }
2673
2674
2675 private static class TideFrequencyDependenceFunction implements TimeVectorFunction {
2676
2677
2678 private final FundamentalNutationArguments arguments;
2679
2680
2681 private final PoissonSeries.CompiledSeries kSeries;
2682
2683
2684
2685
2686
2687
2688
2689
2690
2691 TideFrequencyDependenceFunction(final FundamentalNutationArguments arguments,
2692 final PoissonSeries c20Series,
2693 final PoissonSeries c21Series,
2694 final PoissonSeries s21Series,
2695 final PoissonSeries c22Series,
2696 final PoissonSeries s22Series) {
2697 this.arguments = arguments;
2698 this.kSeries = PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);
2699 }
2700
2701
2702
2703 @Override
2704 public double[] value(final AbsoluteDate date) {
2705 return kSeries.value(arguments.evaluateAll(date));
2706 }
2707
2708
2709 @Override
2710 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
2711 return kSeries.value(arguments.evaluateAll(date));
2712 }
2713
2714 }
2715
2716
2717 private static class EOPTidalCorrection implements TimeVectorFunction {
2718
2719
2720 private final FundamentalNutationArguments arguments;
2721
2722
2723 private final PoissonSeries.CompiledSeries correctionSeries;
2724
2725
2726
2727
2728
2729
2730
2731 EOPTidalCorrection(final FundamentalNutationArguments arguments,
2732 final PoissonSeries xSeries,
2733 final PoissonSeries ySeries,
2734 final PoissonSeries ut1Series) {
2735 this.arguments = arguments;
2736 this.correctionSeries = PoissonSeries.compile(xSeries, ySeries, ut1Series);
2737 }
2738
2739
2740 @Override
2741 public double[] value(final AbsoluteDate date) {
2742 final BodiesElements elements = arguments.evaluateAll(date);
2743 final double[] correction = correctionSeries.value(elements);
2744 final double[] correctionDot = correctionSeries.derivative(elements);
2745 return new double[] {
2746 correction[0],
2747 correction[1],
2748 correction[2],
2749 correctionDot[2] * (-Constants.JULIAN_DAY)
2750 };
2751 }
2752
2753
2754 @Override
2755 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
2756
2757 final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
2758 final T[] correction = correctionSeries.value(elements);
2759 final T[] correctionDot = correctionSeries.derivative(elements);
2760 final T[] a = MathArrays.buildArray(date.getField(), 4);
2761 a[0] = correction[0];
2762 a[1] = correction[1];
2763 a[2] = correction[2];
2764 a[3] = correctionDot[2].multiply(-Constants.JULIAN_DAY);
2765 return a;
2766 }
2767
2768 }
2769
2770 }