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