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