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