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