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