1   /* Copyright 2002-2015 CS Systèmes d'Information
2    * Licensed to CS Systèmes d'Information (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.util.List;
25  
26  import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
27  import org.apache.commons.math3.analysis.interpolation.HermiteInterpolator;
28  import org.apache.commons.math3.util.FastMath;
29  import org.apache.commons.math3.util.MathUtils;
30  import org.orekit.data.BodiesElements;
31  import org.orekit.data.DelaunayArguments;
32  import org.orekit.data.FieldBodiesElements;
33  import org.orekit.data.FundamentalNutationArguments;
34  import org.orekit.data.PoissonSeries;
35  import org.orekit.data.PoissonSeriesParser;
36  import org.orekit.data.PolynomialNutation;
37  import org.orekit.data.PolynomialParser;
38  import org.orekit.data.PolynomialParser.Unit;
39  import org.orekit.data.SimpleTimeStampedTableParser;
40  import org.orekit.errors.OrekitException;
41  import org.orekit.errors.OrekitMessages;
42  import org.orekit.errors.TimeStampedCacheException;
43  import org.orekit.frames.EOPHistory;
44  import org.orekit.frames.PoleCorrection;
45  import org.orekit.time.AbsoluteDate;
46  import org.orekit.time.DateComponents;
47  import org.orekit.time.TimeComponents;
48  import org.orekit.time.TimeFunction;
49  import org.orekit.time.TimeScale;
50  import org.orekit.time.TimeScalesFactory;
51  import org.orekit.time.TimeStamped;
52  
53  
54  /** Supported IERS conventions.
55   * @since 6.0
56   * @author Luc Maisonobe
57   */
58  public enum IERSConventions {
59  
60      /** Constant for IERS 1996 conventions. */
61      IERS_1996 {
62  
63          /** Nutation arguments resources. */
64          private static final String NUTATION_ARGUMENTS = IERS_BASE + "1996/nutation-arguments.txt";
65  
66          /** X series resources. */
67          private static final String X_Y_SERIES         = IERS_BASE + "1996/tab5.4.txt";
68  
69          /** Psi series resources. */
70          private static final String PSI_EPSILON_SERIES = IERS_BASE + "1996/tab5.1.txt";
71  
72          /** Tidal correction for xp, yp series resources. */
73          private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "1996/tab8.4.txt";
74  
75          /** Tidal correction for UT1 resources. */
76          private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "1996/tab8.3.txt";
77  
78          /** Love numbers resources. */
79          private static final String LOVE_NUMBERS = IERS_BASE + "1996/tab6.1.txt";
80  
81          /** Frequency dependence model for k₂₀. */
82          private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2b.txt";
83  
84          /** Frequency dependence model for k₂₁. */
85          private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2a.txt";
86  
87          /** Frequency dependence model for k₂₂. */
88          private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2c.txt";
89  
90          /** {@inheritDoc} */
91          @Override
92          public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
93              throws OrekitException {
94              return new FundamentalNutationArguments(this, timeScale,
95                                                      getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
96          }
97  
98          /** {@inheritDoc} */
99          @Override
100         public TimeFunction<Double> getMeanObliquityFunction() throws OrekitException {
101 
102             // value from chapter 5, page 22
103             final PolynomialNutation<DerivativeStructure> epsilonA =
104                     new PolynomialNutation<DerivativeStructure>(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
105                                                                   -46.8150   * Constants.ARC_SECONDS_TO_RADIANS,
106                                                                    -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
107                                                                     0.001813 * Constants.ARC_SECONDS_TO_RADIANS);
108 
109             return new TimeFunction<Double>() {
110 
111                 /** {@inheritDoc} */
112                 @Override
113                 public Double value(final AbsoluteDate date) {
114                     return epsilonA.value(evaluateTC(date));
115                 }
116 
117             };
118 
119         }
120 
121         /** {@inheritDoc} */
122         @Override
123         public TimeFunction<double[]> getXYSpXY2Function()
124             throws OrekitException {
125 
126             // set up nutation arguments
127             final FundamentalNutationArguments arguments = getNutationArguments(null);
128 
129             // X = 2004.3109″t - 0.42665″t² - 0.198656″t³ + 0.0000140″t⁴
130             //     + 0.00006″t² cos Ω + sin ε0 { Σ [(Ai + Ai' t) sin(ARGUMENT) + Ai'' t cos(ARGUMENT)]}
131             //     + 0.00204″t² sin Ω + 0.00016″t² sin 2(F - D + Ω),
132             final PolynomialNutation<DerivativeStructure> xPolynomial =
133                     new PolynomialNutation<DerivativeStructure>(0,
134                                                                 2004.3109 * Constants.ARC_SECONDS_TO_RADIANS,
135                                                                 -0.42665  * Constants.ARC_SECONDS_TO_RADIANS,
136                                                                 -0.198656 * Constants.ARC_SECONDS_TO_RADIANS,
137                                                                 0.0000140 * Constants.ARC_SECONDS_TO_RADIANS);
138 
139             final double fXCosOm    = 0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
140             final double fXSinOm    = 0.00204 * Constants.ARC_SECONDS_TO_RADIANS;
141             final double fXSin2FDOm = 0.00016 * Constants.ARC_SECONDS_TO_RADIANS;
142             final double sinEps0   = FastMath.sin(getMeanObliquityFunction().value(getNutationReferenceEpoch()));
143 
144             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
145             final PoissonSeriesParser<DerivativeStructure> baseParser =
146                     new PoissonSeriesParser<DerivativeStructure>(12).withFirstDelaunay(1);
147 
148             final PoissonSeriesParser<DerivativeStructure> xParser =
149                     baseParser.
150                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
151                     withSinCos(1, 8, deciMilliAS,  9, deciMilliAS);
152             final PoissonSeries<DerivativeStructure> xSum = xParser.parse(getStream(X_Y_SERIES), X_Y_SERIES);
153 
154             // Y = -0.00013″ - 22.40992″t² + 0.001836″t³ + 0.0011130″t⁴
155             //     + Σ [(Bi + Bi' t) cos(ARGUMENT) + Bi'' t sin(ARGUMENT)]
156             //    - 0.00231″t² cos Ω − 0.00014″t² cos 2(F - D + Ω)
157             final PolynomialNutation<DerivativeStructure> yPolynomial =
158                     new PolynomialNutation<DerivativeStructure>(-0.00013  * Constants.ARC_SECONDS_TO_RADIANS,
159                                                                 0.0,
160                                                                 -22.40992 * Constants.ARC_SECONDS_TO_RADIANS,
161                                                                 0.001836  * Constants.ARC_SECONDS_TO_RADIANS,
162                                                                 0.0011130 * Constants.ARC_SECONDS_TO_RADIANS);
163 
164             final double fYCosOm    = -0.00231 * Constants.ARC_SECONDS_TO_RADIANS;
165             final double fYCos2FDOm = -0.00014 * Constants.ARC_SECONDS_TO_RADIANS;
166 
167             final PoissonSeriesParser<DerivativeStructure> yParser =
168                     baseParser.
169                     withSinCos(0, -1, deciMilliAS, 10, deciMilliAS).
170                     withSinCos(1, 12, deciMilliAS, 11, deciMilliAS);
171             final PoissonSeries<DerivativeStructure> ySum = yParser.parse(getStream(X_Y_SERIES), X_Y_SERIES);
172 
173             @SuppressWarnings("unchecked")
174             final PoissonSeries.CompiledSeries<DerivativeStructure> xySum =
175                     PoissonSeries.compile(xSum, ySum);
176 
177             // s = -XY/2 + 0.00385″t - 0.07259″t³ - 0.00264″ sin Ω - 0.00006″ sin 2Ω
178             //     + 0.00074″t² sin Ω + 0.00006″t² sin 2(F - D + Ω)
179             final double fST          =  0.00385 * Constants.ARC_SECONDS_TO_RADIANS;
180             final double fST3         = -0.07259 * Constants.ARC_SECONDS_TO_RADIANS;
181             final double fSSinOm      = -0.00264 * Constants.ARC_SECONDS_TO_RADIANS;
182             final double fSSin2Om     = -0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
183             final double fST2SinOm    =  0.00074 * Constants.ARC_SECONDS_TO_RADIANS;
184             final double fST2Sin2FDOm =  0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
185 
186             return new TimeFunction<double[]>() {
187 
188                 /** {@inheritDoc} */
189                 @Override
190                 public double[] value(final AbsoluteDate date) {
191 
192                     final BodiesElements elements = arguments.evaluateAll(date);
193                     final double[] xy             = xySum.value(elements);
194 
195                     final double omega     = elements.getOmega();
196                     final double f         = elements.getF();
197                     final double d         = elements.getD();
198                     final double t         = elements.getTC();
199 
200                     final double cosOmega  = FastMath.cos(omega);
201                     final double sinOmega  = FastMath.sin(omega);
202                     final double sin2Omega = FastMath.sin(2 * omega);
203                     final double cos2FDOm  = FastMath.cos(2 * (f - d + omega));
204                     final double sin2FDOm  = FastMath.sin(2 * (f - d + omega));
205 
206                     final double x = xPolynomial.value(t) + sinEps0 * xy[0] +
207                             t * t * (fXCosOm * cosOmega + fXSinOm * sinOmega + fXSin2FDOm * cos2FDOm);
208                     final double y = yPolynomial.value(t) + xy[1] +
209                             t * t * (fYCosOm * cosOmega + fYCos2FDOm * cos2FDOm);
210                     final double sPxy2 = fSSinOm * sinOmega + fSSin2Om * sin2Omega +
211                             t * (fST + t * (fST2SinOm * sinOmega + fST2Sin2FDOm * sin2FDOm + t * fST3));
212 
213                     return new double[] {
214                         x, y, sPxy2
215                     };
216 
217                 }
218 
219             };
220 
221         }
222 
223         /** {@inheritDoc} */
224         @Override
225         public TimeFunction<double[]> getPrecessionFunction() throws OrekitException {
226 
227             // set up the conventional polynomials
228             // the following values are from Lieske et al. paper:
229             // Expressions for the precession quantities based upon the IAU(1976) system of astronomical constants
230             // http://articles.adsabs.harvard.edu/full/1977A%26A....58....1L
231             // also available as equation 30 in IERS 2003 conventions
232             final PolynomialNutation<DerivativeStructure> psiA =
233                     new PolynomialNutation<DerivativeStructure>(   0.0,
234                                                                 5038.7784   * Constants.ARC_SECONDS_TO_RADIANS,
235                                                                   -1.07259  * Constants.ARC_SECONDS_TO_RADIANS,
236                                                                   -0.001147 * Constants.ARC_SECONDS_TO_RADIANS);
237             final PolynomialNutation<DerivativeStructure> omegaA =
238                     new PolynomialNutation<DerivativeStructure>(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
239                                                                  0.0,
240                                                                  0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
241                                                                 -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
242             final PolynomialNutation<DerivativeStructure> chiA =
243                     new PolynomialNutation<DerivativeStructure>( 0.0,
244                                                                 10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
245                                                                 -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
246                                                                 -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);
247 
248             return new TimeFunction<double[]>() {
249                 /** {@inheritDoc} */
250                 @Override
251                 public double[] value(final AbsoluteDate date) {
252                     final double tc = evaluateTC(date);
253                     return new double[] {
254                         psiA.value(tc), omegaA.value(tc), chiA.value(tc)
255                     };
256                 }
257             };
258 
259         }
260 
261         /** {@inheritDoc} */
262         @Override
263         public TimeFunction<double[]> getNutationFunction()
264             throws OrekitException {
265 
266             // set up nutation arguments
267             final FundamentalNutationArguments arguments = getNutationArguments(null);
268 
269             // set up Poisson series
270             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
271             final PoissonSeriesParser<DerivativeStructure> baseParser =
272                     new PoissonSeriesParser<DerivativeStructure>(10).withFirstDelaunay(1);
273 
274             final PoissonSeriesParser<DerivativeStructure> psiParser =
275                     baseParser.
276                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
277                     withSinCos(1, 8, deciMilliAS, -1, deciMilliAS);
278             final PoissonSeries<DerivativeStructure> psiSeries = psiParser.parse(getStream(PSI_EPSILON_SERIES), PSI_EPSILON_SERIES);
279 
280             final PoissonSeriesParser<DerivativeStructure> epsilonParser =
281                     baseParser.
282                     withSinCos(0, -1, deciMilliAS, 9, deciMilliAS).
283                     withSinCos(1, -1, deciMilliAS, 10, deciMilliAS);
284             final PoissonSeries<DerivativeStructure> epsilonSeries = epsilonParser.parse(getStream(PSI_EPSILON_SERIES), PSI_EPSILON_SERIES);
285 
286             @SuppressWarnings("unchecked")
287             final PoissonSeries.CompiledSeries<DerivativeStructure> psiEpsilonSeries =
288                     PoissonSeries.compile(psiSeries, epsilonSeries);
289 
290             return new TimeFunction<double[]>() {
291                 /** {@inheritDoc} */
292                 @Override
293                 public double[] value(final AbsoluteDate date) {
294                     final BodiesElements elements = arguments.evaluateAll(date);
295                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
296                     return new double[] {
297                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements)
298                     };
299                 }
300             };
301 
302         }
303 
304         /** {@inheritDoc} */
305         @Override
306         public TimeFunction<DerivativeStructure> getGMSTFunction(final TimeScale ut1)
307             throws OrekitException {
308 
309             // Radians per second of time
310             final double radiansPerSecond = MathUtils.TWO_PI / Constants.JULIAN_DAY;
311 
312             // constants from IERS 1996 page 21
313             // the underlying model is IAU 1982 GMST-UT1
314             final AbsoluteDate gmstReference =
315                 new AbsoluteDate(DateComponents.J2000_EPOCH, TimeComponents.H12, TimeScalesFactory.getTAI());
316             final double gmst0 = 24110.54841;
317             final double gmst1 = 8640184.812866;
318             final double gmst2 = 0.093104;
319             final double gmst3 = -6.2e-6;
320 
321             return new TimeFunction<DerivativeStructure>() {
322 
323                 /** {@inheritDoc} */
324                 @Override
325                 public DerivativeStructure value(final AbsoluteDate date) {
326 
327                     // offset in Julian centuries from J2000 epoch (UT1 scale)
328                     final double dtai = date.durationFrom(gmstReference);
329                     final DerivativeStructure tut1 =
330                             new DerivativeStructure(1, 1, dtai + ut1.offsetFromTAI(date), 1.0);
331                     final DerivativeStructure tt = tut1.divide(Constants.JULIAN_CENTURY);
332 
333                     // Seconds in the day, adjusted by 12 hours because the
334                     // UT1 is supplied as a Julian date beginning at noon.
335                     final DerivativeStructure sd = tut1.add(Constants.JULIAN_DAY / 2).remainder(Constants.JULIAN_DAY);
336 
337                     // compute Greenwich mean sidereal time, in radians
338                     return tt.multiply(gmst3).add(gmst2).multiply(tt).add(gmst1).multiply(tt).add(gmst0).add(sd).multiply(radiansPerSecond);
339 
340                 }
341 
342             };
343 
344         }
345 
346         /** {@inheritDoc} */
347         @Override
348         public TimeFunction<DerivativeStructure> getGASTFunction(final TimeScale ut1,
349                                                                  final EOPHistory eopHistory)
350             throws OrekitException {
351 
352             // obliquity
353             final TimeFunction<Double> epsilonA = getMeanObliquityFunction();
354 
355             // GMST function
356             final TimeFunction<DerivativeStructure> gmst = getGMSTFunction(ut1);
357 
358             // nutation function
359             final TimeFunction<double[]> nutation = getNutationFunction();
360 
361             return new TimeFunction<DerivativeStructure>() {
362 
363                 /** {@inheritDoc} */
364                 @Override
365                 public DerivativeStructure value(final AbsoluteDate date) {
366 
367                     // compute equation of equinoxes
368                     final double[] angles = nutation.value(date);
369                     double deltaPsi = angles[0];
370                     if (eopHistory != null) {
371                         deltaPsi += eopHistory.getEquinoxNutationCorrection(date)[0];
372                     }
373                     final double eqe = deltaPsi  * FastMath.cos(epsilonA.value(date)) + angles[2];
374 
375                     // add mean sidereal time and equation of equinoxes
376                     return gmst.value(date).add(eqe);
377 
378                 }
379 
380             };
381 
382         }
383 
384         /** {@inheritDoc} */
385         @Override
386         public TimeFunction<double[]> getEOPTidalCorrection()
387             throws OrekitException {
388 
389             // set up nutation arguments
390             // BEWARE! Using TT as the time scale here and not UT1 is intentional!
391             // as this correction is used to compute UT1 itself, it is not surprising we cannot use UT1 yet,
392             // however, using the close UTC as would seem logical make the comparison with interp.f from IERS fail
393             // looking in the interp.f code, the same TT scale is used for both Delaunay and gamma argument
394             final FundamentalNutationArguments arguments = getNutationArguments(TimeScalesFactory.getTT());
395 
396             // set up Poisson series
397             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
398             final PoissonSeriesParser<DerivativeStructure> xyParser = new PoissonSeriesParser<DerivativeStructure>(17).
399                     withOptionalColumn(1).
400                     withGamma(7).
401                     withFirstDelaunay(2);
402             final PoissonSeries<DerivativeStructure> xSeries =
403                     xyParser.
404                     withSinCos(0, 14, milliAS, 15, milliAS).
405                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
406             final PoissonSeries<DerivativeStructure> ySeries =
407                     xyParser.
408                     withSinCos(0, 16, milliAS, 17, milliAS).
409                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES),
410                                                          TIDAL_CORRECTION_XP_YP_SERIES);
411 
412             final double deciMilliS = 1.0e-4;
413             final PoissonSeriesParser<DerivativeStructure> ut1Parser = new PoissonSeriesParser<DerivativeStructure>(17).
414                     withOptionalColumn(1).
415                     withGamma(7).
416                     withFirstDelaunay(2).
417                     withSinCos(0, 16, deciMilliS, 17, deciMilliS);
418             final PoissonSeries<DerivativeStructure> ut1Series =
419                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);
420 
421             @SuppressWarnings("unchecked")
422             final PoissonSeries.CompiledSeries<DerivativeStructure> correctionSeries =
423                 PoissonSeries.compile(xSeries, ySeries, ut1Series);
424 
425             return new TimeFunction<double[]>() {
426                 /** {@inheritDoc} */
427                 @Override
428                 public double[] value(final AbsoluteDate date) {
429                     final FieldBodiesElements<DerivativeStructure> elements =
430                             arguments.evaluateDerivative(date);
431                     final DerivativeStructure[] correction = correctionSeries.value(elements);
432                     return new double[] {
433                         correction[0].getValue(),
434                         correction[1].getValue(),
435                         correction[2].getValue(),
436                         -correction[2].getPartialDerivative(1) * Constants.JULIAN_DAY
437                     };
438                 }
439             };
440 
441         }
442 
443         /** {@inheritDoc} */
444         public LoveNumbers getLoveNumbers() throws OrekitException {
445             return loadLoveNumbers(LOVE_NUMBERS);
446         }
447 
448         /** {@inheritDoc} */
449         public TimeFunction<double[]> getTideFrequencyDependenceFunction(final TimeScale ut1)
450             throws OrekitException {
451 
452             // set up nutation arguments
453             final FundamentalNutationArguments arguments = getNutationArguments(ut1);
454 
455             // set up Poisson series
456             final PoissonSeriesParser<DerivativeStructure> k20Parser =
457                     new PoissonSeriesParser<DerivativeStructure>(18).
458                         withOptionalColumn(1).
459                         withDoodson(4, 2).
460                         withFirstDelaunay(10);
461             final PoissonSeriesParser<DerivativeStructure> k21Parser =
462                     new PoissonSeriesParser<DerivativeStructure>(18).
463                         withOptionalColumn(1).
464                         withDoodson(4, 3).
465                         withFirstDelaunay(10);
466             final PoissonSeriesParser<DerivativeStructure> k22Parser =
467                     new PoissonSeriesParser<DerivativeStructure>(16).
468                         withOptionalColumn(1).
469                         withDoodson(4, 2).
470                         withFirstDelaunay(10);
471 
472             final double pico = 1.0e-12;
473             final PoissonSeries<DerivativeStructure> c20Series =
474                     k20Parser.
475                     // TODO: check sign!
476 //                  withSinCos(0, 18, pico, 16, pico).
477                   withSinCos(0, 18, -pico, 16, pico).
478                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
479             final PoissonSeries<DerivativeStructure> c21Series =
480                     k21Parser.
481                     withSinCos(0, 17, pico, 18, pico).
482                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
483             final PoissonSeries<DerivativeStructure> s21Series =
484                     k21Parser.
485                     withSinCos(0, 18, -pico, 17, pico).
486                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
487             final PoissonSeries<DerivativeStructure> c22Series =
488                     k22Parser.
489                     withSinCos(0, -1, pico, 16, pico).
490                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
491             final PoissonSeries<DerivativeStructure> s22Series =
492                     k22Parser.
493                     withSinCos(0, 16, -pico, -1, pico).
494                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
495 
496             @SuppressWarnings("unchecked")
497             final PoissonSeries.CompiledSeries<DerivativeStructure> kSeries =
498                 PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);
499 
500             return new TimeFunction<double[]>() {
501                 /** {@inheritDoc} */
502                 @Override
503                 public double[] value(final AbsoluteDate date) {
504                     return kSeries.value(arguments.evaluateAll(date));
505                 }
506             };
507 
508         }
509 
510         /** {@inheritDoc} */
511         @Override
512         public double getPermanentTide() throws OrekitException {
513             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
514         }
515 
516         /** {@inheritDoc} */
517         @Override
518         public TimeFunction<double[]> getSolidPoleTide(final EOPHistory eopHistory) {
519 
520             // constants from IERS 1996 page 47
521             final double globalFactor = -1.348e-9 / Constants.ARC_SECONDS_TO_RADIANS;
522             final double coupling     =  0.00112;
523 
524             return new TimeFunction<double[]>() {
525                 /** {@inheritDoc} */
526                 @Override
527                 public double[] value(final AbsoluteDate date) {
528                     final PoleCorrection pole = eopHistory.getPoleCorrection(date);
529                     return new double[] {
530                         globalFactor * (pole.getXp() + coupling * pole.getYp()),
531                         globalFactor * (coupling * pole.getXp() - pole.getYp()),
532                     };
533                 }
534             };
535         }
536 
537         /** {@inheritDoc} */
538         @Override
539         public TimeFunction<double[]> getOceanPoleTide(final EOPHistory eopHistory)
540             throws OrekitException {
541 
542             return new TimeFunction<double[]>() {
543                 /** {@inheritDoc} */
544                 @Override
545                 public double[] value(final AbsoluteDate date) {
546                     // there are no model for ocean pole tide prior to conventions 2010
547                     return new double[] {
548                         0.0, 0.0
549                     };
550                 }
551             };
552         }
553 
554     },
555 
556     /** Constant for IERS 2003 conventions. */
557     IERS_2003 {
558 
559         /** Nutation arguments resources. */
560         private static final String NUTATION_ARGUMENTS = IERS_BASE + "2003/nutation-arguments.txt";
561 
562         /** X series resources. */
563         private static final String X_SERIES           = IERS_BASE + "2003/tab5.2a.txt";
564 
565         /** Y series resources. */
566         private static final String Y_SERIES           = IERS_BASE + "2003/tab5.2b.txt";
567 
568         /** S series resources. */
569         private static final String S_SERIES           = IERS_BASE + "2003/tab5.2c.txt";
570 
571         /** Luni-solar series resources. */
572         private static final String LUNI_SOLAR_SERIES  = IERS_BASE + "2003/tab5.3a-first-table.txt";
573 
574         /** Planetary series resources. */
575         private static final String PLANETARY_SERIES   = IERS_BASE + "2003/tab5.3b.txt";
576 
577         /** Greenwhich sidereal time series resources. */
578         private static final String GST_SERIES         = IERS_BASE + "2003/tab5.4.txt";
579 
580         /** Tidal correction for xp, yp series resources. */
581         private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "2003/tab8.2ab.txt";
582 
583         /** Tidal correction for UT1 resources. */
584         private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "2003/tab8.3ab.txt";
585 
586         /** Love numbers resources. */
587         private static final String LOVE_NUMBERS = IERS_BASE + "2003/tab6.1.txt";
588 
589         /** Frequency dependence model for k₂₀. */
590         private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3b.txt";
591 
592         /** Frequency dependence model for k₂₁. */
593         private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3a.txt";
594 
595         /** Frequency dependence model for k₂₂. */
596         private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3c.txt";
597 
598         /** Annual pole table. */
599         private static final String ANNUAL_POLE = IERS_BASE + "2003/annual.pole";
600 
601         /** {@inheritDoc} */
602         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
603             throws OrekitException {
604             return new FundamentalNutationArguments(this, timeScale,
605                                                     getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
606         }
607 
608         /** {@inheritDoc} */
609         @Override
610         public TimeFunction<Double> getMeanObliquityFunction() throws OrekitException {
611 
612             // epsilon 0 value from chapter 5, page 41, other terms from equation 32 page 45
613             final PolynomialNutation<DerivativeStructure> epsilonA =
614                     new PolynomialNutation<DerivativeStructure>(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
615                                                                   -46.84024  * Constants.ARC_SECONDS_TO_RADIANS,
616                                                                    -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
617                                                                     0.001813 * Constants.ARC_SECONDS_TO_RADIANS);
618 
619             return new TimeFunction<Double>() {
620 
621                 /** {@inheritDoc} */
622                 @Override
623                 public Double value(final AbsoluteDate date) {
624                     return epsilonA.value(evaluateTC(date));
625                 }
626 
627             };
628 
629         }
630 
631         /** {@inheritDoc} */
632         @Override
633         public TimeFunction<double[]> getXYSpXY2Function()
634             throws OrekitException {
635 
636             // set up nutation arguments
637             final FundamentalNutationArguments arguments = getNutationArguments(null);
638 
639             // set up Poisson series
640             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
641             final PoissonSeriesParser<DerivativeStructure> parser =
642                     new PoissonSeriesParser<DerivativeStructure>(17).
643                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
644                         withFirstDelaunay(4).
645                         withFirstPlanetary(9).
646                         withSinCos(0, 2, microAS, 3, microAS);
647 
648             final PoissonSeries<DerivativeStructure> xSeries = parser.parse(getStream(X_SERIES), X_SERIES);
649             final PoissonSeries<DerivativeStructure> ySeries = parser.parse(getStream(Y_SERIES), Y_SERIES);
650             final PoissonSeries<DerivativeStructure> sSeries = parser.parse(getStream(S_SERIES), S_SERIES);
651             @SuppressWarnings("unchecked")
652             final PoissonSeries.CompiledSeries<DerivativeStructure> xys = PoissonSeries.compile(xSeries, ySeries, sSeries);
653 
654             // create a function evaluating the series
655             return new TimeFunction<double[]>() {
656 
657                 /** {@inheritDoc} */
658                 @Override
659                 public double[] value(final AbsoluteDate date) {
660                     return xys.value(arguments.evaluateAll(date));
661                 }
662 
663             };
664 
665         }
666 
667 
668         /** {@inheritDoc} */
669         @Override
670         public TimeFunction<double[]> getPrecessionFunction() throws OrekitException {
671 
672             // set up the conventional polynomials
673             // the following values are from equation 32 in IERS 2003 conventions
674             final PolynomialNutation<DerivativeStructure> psiA =
675                     new PolynomialNutation<DerivativeStructure>(    0.0,
676                                                                  5038.47875   * Constants.ARC_SECONDS_TO_RADIANS,
677                                                                    -1.07259   * Constants.ARC_SECONDS_TO_RADIANS,
678                                                                    -0.001147  * Constants.ARC_SECONDS_TO_RADIANS);
679             final PolynomialNutation<DerivativeStructure> omegaA =
680                     new PolynomialNutation<DerivativeStructure>(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
681                                                                 -0.02524   * Constants.ARC_SECONDS_TO_RADIANS,
682                                                                  0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
683                                                                 -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
684             final PolynomialNutation<DerivativeStructure> chiA =
685                     new PolynomialNutation<DerivativeStructure>( 0.0,
686                                                                 10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
687                                                                 -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
688                                                                 -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);
689 
690             return new TimeFunction<double[]>() {
691                 /** {@inheritDoc} */
692                 @Override
693                 public double[] value(final AbsoluteDate date) {
694                     final double tc = evaluateTC(date);
695                     return new double[] {
696                         psiA.value(tc), omegaA.value(tc), chiA.value(tc)
697                     };
698                 }
699             };
700 
701         }
702 
703         /** {@inheritDoc} */
704         @Override
705         public TimeFunction<double[]> getNutationFunction()
706             throws OrekitException {
707 
708             // set up nutation arguments
709             final FundamentalNutationArguments arguments = getNutationArguments(null);
710 
711             // set up Poisson series
712             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
713             final PoissonSeriesParser<DerivativeStructure> luniSolarParser =
714                     new PoissonSeriesParser<DerivativeStructure>(14).withFirstDelaunay(1);
715             final PoissonSeriesParser<DerivativeStructure> luniSolarPsiParser =
716                     luniSolarParser.
717                     withSinCos(0, 7, milliAS, 11, milliAS).
718                     withSinCos(1, 8, milliAS, 12, milliAS);
719             final PoissonSeries<DerivativeStructure> psiLuniSolarSeries =
720                     luniSolarPsiParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);
721             final PoissonSeriesParser<DerivativeStructure> luniSolarEpsilonParser =
722                     luniSolarParser.
723                     withSinCos(0, 13, milliAS, 9, milliAS).
724                     withSinCos(1, 14, milliAS, 10, milliAS);
725             final PoissonSeries<DerivativeStructure> epsilonLuniSolarSeries =
726                     luniSolarEpsilonParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);
727 
728             final PoissonSeriesParser<DerivativeStructure> planetaryParser =
729                     new PoissonSeriesParser<DerivativeStructure>(21).
730                         withFirstDelaunay(2).
731                         withFirstPlanetary(7);
732             final PoissonSeriesParser<DerivativeStructure> planetaryPsiParser =
733                     planetaryParser.withSinCos(0, 17, milliAS, 18, milliAS);
734             final PoissonSeries<DerivativeStructure> psiPlanetarySeries =
735                     planetaryPsiParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);
736             final PoissonSeriesParser<DerivativeStructure> planetaryEpsilonParser =
737                     planetaryParser.withSinCos(0, 19, milliAS, 20, milliAS);
738             final PoissonSeries<DerivativeStructure> epsilonPlanetarySeries =
739                     planetaryEpsilonParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);
740 
741             @SuppressWarnings("unchecked")
742             final PoissonSeries.CompiledSeries<DerivativeStructure> luniSolarSeries =
743                     PoissonSeries.compile(psiLuniSolarSeries, epsilonLuniSolarSeries);
744             @SuppressWarnings("unchecked")
745             final PoissonSeries.CompiledSeries<DerivativeStructure> planetarySeries =
746                     PoissonSeries.compile(psiPlanetarySeries, epsilonPlanetarySeries);
747 
748             return new TimeFunction<double[]>() {
749                 /** {@inheritDoc} */
750                 @Override
751                 public double[] value(final AbsoluteDate date) {
752                     final BodiesElements elements = arguments.evaluateAll(date);
753                     final double[] luniSolar = luniSolarSeries.value(elements);
754                     final double[] planetary = planetarySeries.value(elements);
755                     return new double[] {
756                         luniSolar[0] + planetary[0], luniSolar[1] + planetary[1],
757                         IAU1994ResolutionC7.value(elements)
758                     };
759                 }
760             };
761 
762         }
763 
764         /** {@inheritDoc} */
765         @Override
766         public TimeFunction<DerivativeStructure> getGMSTFunction(final TimeScale ut1)
767             throws OrekitException {
768 
769             // Earth Rotation Angle
770             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);
771 
772             // Polynomial part of the apparent sidereal time series
773             // which is the opposite of Equation of Origins (EO)
774             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
775             final PoissonSeriesParser<DerivativeStructure> parser =
776                     new PoissonSeriesParser<DerivativeStructure>(17).
777                         withFirstDelaunay(4).
778                         withFirstPlanetary(9).
779                         withSinCos(0, 2, microAS, 3, microAS).
780                         withPolynomialPart('t', Unit.ARC_SECONDS);
781             final PolynomialNutation<DerivativeStructure> minusEO =
782                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();
783 
784             // create a function evaluating the series
785             return new TimeFunction<DerivativeStructure>() {
786 
787                 /** {@inheritDoc} */
788                 @Override
789                 public DerivativeStructure value(final AbsoluteDate date) {
790                     return era.value(date).add(minusEO.value(dsEvaluateTC(date)));
791                 }
792 
793             };
794 
795         }
796 
797         /** {@inheritDoc} */
798         @Override
799         public TimeFunction<DerivativeStructure> getGASTFunction(final TimeScale ut1,
800                                                                  final EOPHistory eopHistory)
801             throws OrekitException {
802 
803             // set up nutation arguments
804             final FundamentalNutationArguments arguments = getNutationArguments(null);
805 
806             // mean obliquity function
807             final TimeFunction<Double> epsilon = getMeanObliquityFunction();
808 
809             // set up Poisson series
810             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
811             final PoissonSeriesParser<DerivativeStructure> luniSolarPsiParser =
812                     new PoissonSeriesParser<DerivativeStructure>(14).
813                         withFirstDelaunay(1).
814                         withSinCos(0, 7, milliAS, 11, milliAS).
815                         withSinCos(1, 8, milliAS, 12, milliAS);
816             final PoissonSeries<DerivativeStructure> psiLuniSolarSeries =
817                     luniSolarPsiParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);
818 
819             final PoissonSeriesParser<DerivativeStructure> planetaryPsiParser =
820                     new PoissonSeriesParser<DerivativeStructure>(21).
821                         withFirstDelaunay(2).
822                         withFirstPlanetary(7).
823                         withSinCos(0, 17, milliAS, 18, milliAS);
824             final PoissonSeries<DerivativeStructure> psiPlanetarySeries =
825                     planetaryPsiParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);
826 
827             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
828             final PoissonSeriesParser<DerivativeStructure> gstParser =
829                     new PoissonSeriesParser<DerivativeStructure>(17).
830                         withFirstDelaunay(4).
831                         withFirstPlanetary(9).
832                         withSinCos(0, 2, microAS, 3, microAS).
833                         withPolynomialPart('t', Unit.ARC_SECONDS);
834             final PoissonSeries<DerivativeStructure> gstSeries = gstParser.parse(getStream(GST_SERIES), GST_SERIES);
835             @SuppressWarnings("unchecked")
836             final PoissonSeries.CompiledSeries<DerivativeStructure> psiGstSeries =
837                     PoissonSeries.compile(psiLuniSolarSeries, psiPlanetarySeries, gstSeries);
838 
839             // ERA function
840             final TimeFunction<DerivativeStructure> era = getEarthOrientationAngleFunction(ut1);
841 
842             return new TimeFunction<DerivativeStructure>() {
843 
844                 /** {@inheritDoc} */
845                 @Override
846                 public DerivativeStructure value(final AbsoluteDate date) {
847 
848                     // evaluate equation of origins
849                     final BodiesElements elements = arguments.evaluateAll(date);
850                     final double[] angles = psiGstSeries.value(elements);
851                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
852                     final double deltaPsi = angles[0] + angles[1] + ddPsi;
853                     final double epsilonA = epsilon.value(date);
854 
855                     // subtract equation of origin from EA
856                     // (hence add the series above which have the sign included)
857                     return era.value(date).add(deltaPsi * FastMath.cos(epsilonA) + angles[2]);
858 
859                 }
860 
861             };
862 
863         }
864 
865         /** {@inheritDoc} */
866         @Override
867         public TimeFunction<double[]> getEOPTidalCorrection()
868             throws OrekitException {
869 
870             // set up nutation arguments
871             // BEWARE! Using TT as the time scale here and not UT1 is intentional!
872             // as this correction is used to compute UT1 itself, it is not surprising we cannot use UT1 yet,
873             // however, using the close UTC as would seem logical make the comparison with interp.f from IERS fail
874             // looking in the interp.f code, the same TT scale is used for both Delaunay and gamma argument
875             final FundamentalNutationArguments arguments = getNutationArguments(TimeScalesFactory.getTT());
876 
877             // set up Poisson series
878             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
879             final PoissonSeriesParser<DerivativeStructure> xyParser = new PoissonSeriesParser<DerivativeStructure>(13).
880                     withOptionalColumn(1).
881                     withGamma(2).
882                     withFirstDelaunay(3);
883             final PoissonSeries<DerivativeStructure> xSeries =
884                     xyParser.
885                     withSinCos(0, 10, microAS, 11, microAS).
886                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
887             final PoissonSeries<DerivativeStructure> ySeries =
888                     xyParser.
889                     withSinCos(0, 12, microAS, 13, microAS).
890                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
891 
892             final double microS = 1.0e-6;
893             final PoissonSeriesParser<DerivativeStructure> ut1Parser = new PoissonSeriesParser<DerivativeStructure>(11).
894                     withOptionalColumn(1).
895                     withGamma(2).
896                     withFirstDelaunay(3).
897                     withSinCos(0, 10, microS, 11, microS);
898             final PoissonSeries<DerivativeStructure> ut1Series =
899                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);
900 
901             @SuppressWarnings("unchecked")
902             final PoissonSeries.CompiledSeries<DerivativeStructure> correctionSeries =
903                 PoissonSeries.compile(xSeries, ySeries, ut1Series);
904 
905             return new TimeFunction<double[]>() {
906                 /** {@inheritDoc} */
907                 @Override
908                 public double[] value(final AbsoluteDate date) {
909                     final FieldBodiesElements<DerivativeStructure> elements =
910                             arguments.evaluateDerivative(date);
911                     final DerivativeStructure[] correction = correctionSeries.value(elements);
912                     return new double[] {
913                         correction[0].getValue(),
914                         correction[1].getValue(),
915                         correction[2].getValue(),
916                         -correction[2].getPartialDerivative(1) * Constants.JULIAN_DAY
917                     };
918                 }
919             };
920 
921         }
922 
923         /** {@inheritDoc} */
924         public LoveNumbers getLoveNumbers() throws OrekitException {
925             return loadLoveNumbers(LOVE_NUMBERS);
926         }
927 
928         /** {@inheritDoc} */
929         public TimeFunction<double[]> getTideFrequencyDependenceFunction(final TimeScale ut1)
930             throws OrekitException {
931 
932             // set up nutation arguments
933             final FundamentalNutationArguments arguments = getNutationArguments(ut1);
934 
935             // set up Poisson series
936             final PoissonSeriesParser<DerivativeStructure> k20Parser =
937                     new PoissonSeriesParser<DerivativeStructure>(18).
938                         withOptionalColumn(1).
939                         withDoodson(4, 2).
940                         withFirstDelaunay(10);
941             final PoissonSeriesParser<DerivativeStructure> k21Parser =
942                     new PoissonSeriesParser<DerivativeStructure>(18).
943                         withOptionalColumn(1).
944                         withDoodson(4, 3).
945                         withFirstDelaunay(10);
946             final PoissonSeriesParser<DerivativeStructure> k22Parser =
947                     new PoissonSeriesParser<DerivativeStructure>(16).
948                         withOptionalColumn(1).
949                         withDoodson(4, 2).
950                         withFirstDelaunay(10);
951 
952             final double pico = 1.0e-12;
953             final PoissonSeries<DerivativeStructure> c20Series =
954                     k20Parser.
955                     // TODO: check sign!
956 //                  withSinCos(0, 18, pico, 16, pico).
957                   withSinCos(0, 18, -pico, 16, pico).
958                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
959             final PoissonSeries<DerivativeStructure> c21Series =
960                     k21Parser.
961                     withSinCos(0, 17, pico, 18, pico).
962                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
963             final PoissonSeries<DerivativeStructure> s21Series =
964                     k21Parser.
965                     withSinCos(0, 18, -pico, 17, pico).
966                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
967             final PoissonSeries<DerivativeStructure> c22Series =
968                     k22Parser.
969                     withSinCos(0, -1, pico, 16, pico).
970                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
971             final PoissonSeries<DerivativeStructure> s22Series =
972                     k22Parser.
973                     withSinCos(0, 16, -pico, -1, pico).
974                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
975 
976             @SuppressWarnings("unchecked")
977             final PoissonSeries.CompiledSeries<DerivativeStructure> kSeries =
978                 PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);
979 
980             return new TimeFunction<double[]>() {
981                 /** {@inheritDoc} */
982                 @Override
983                 public double[] value(final AbsoluteDate date) {
984                     return kSeries.value(arguments.evaluateAll(date));
985                 }
986             };
987 
988         }
989 
990         /** {@inheritDoc} */
991         @Override
992         public double getPermanentTide() throws OrekitException {
993             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
994         }
995 
996         /** {@inheritDoc} */
997         @Override
998         public TimeFunction<double[]> getSolidPoleTide(final EOPHistory eopHistory)
999             throws OrekitException {
1000 
1001             // annual pole from ftp://tai.bipm.org/iers/conv2003/chapter7/annual.pole
1002             final TimeScale utc = TimeScalesFactory.getUTC();
1003             final SimpleTimeStampedTableParser.RowConverter<MeanPole> converter =
1004                 new SimpleTimeStampedTableParser.RowConverter<MeanPole>() {
1005                     /** {@inheritDoc} */
1006                     @Override
1007                     public MeanPole convert(final double[] rawFields) throws OrekitException {
1008                         return new MeanPole(new AbsoluteDate((int) rawFields[0], 1, 1, utc),
1009                                             rawFields[1] * Constants.ARC_SECONDS_TO_RADIANS,
1010                                             rawFields[2] * Constants.ARC_SECONDS_TO_RADIANS);
1011                     }
1012                 };
1013             final SimpleTimeStampedTableParser<MeanPole> parser =
1014                     new SimpleTimeStampedTableParser<MeanPole>(3, converter);
1015             final List<MeanPole> annualPoleList = parser.parse(getStream(ANNUAL_POLE), ANNUAL_POLE);
1016             final AbsoluteDate firstAnnualPoleDate = annualPoleList.get(0).getDate();
1017             final AbsoluteDate lastAnnualPoleDate  = annualPoleList.get(annualPoleList.size() - 1).getDate();
1018             final ImmutableTimeStampedCache<MeanPole> annualCache =
1019                     new ImmutableTimeStampedCache<MeanPole>(2, annualPoleList);
1020 
1021             // polynomial extension from IERS 2003, section 7.1.4, equations 23a and 23b
1022             final double xp0    = 0.054   * Constants.ARC_SECONDS_TO_RADIANS;
1023             final double xp0Dot = 0.00083 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
1024             final double yp0    = 0.357   * Constants.ARC_SECONDS_TO_RADIANS;
1025             final double yp0Dot = 0.00395 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
1026 
1027             // constants from IERS 2003, section 6.2
1028             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
1029             final double ratio        =  0.00115;
1030 
1031             return new TimeFunction<double[]>() {
1032                 /** {@inheritDoc} */
1033                 @Override
1034                 public double[] value(final AbsoluteDate date) {
1035 
1036                     // we can't compute anything before the range covered by the annual pole file
1037                     if (date.compareTo(firstAnnualPoleDate) <= 0) {
1038                         return new double[] {
1039                             0.0, 0.0
1040                         };
1041                     }
1042 
1043                     // evaluate mean pole
1044                     double meanPoleX = 0;
1045                     double meanPoleY = 0;
1046                     if (date.compareTo(lastAnnualPoleDate) <= 0) {
1047                         // we are within the range covered by the annual pole file,
1048                         // we interpolate within it
1049                         try {
1050                             final List<MeanPole> neighbors = annualCache.getNeighbors(date);
1051                             final HermiteInterpolator interpolator = new HermiteInterpolator();
1052                             for (final MeanPole neighbor : neighbors) {
1053                                 interpolator.addSamplePoint(neighbor.getDate().durationFrom(date),
1054                                                             new double[] {
1055                                                                 neighbor.getX(), neighbor.getY()
1056                                                             });
1057                             }
1058                             final double[] interpolated = interpolator.value(0);
1059                             meanPoleX = interpolated[0];
1060                             meanPoleY = interpolated[1];
1061                         } catch (TimeStampedCacheException tsce) {
1062                             // this should never happen
1063                             throw OrekitException.createInternalError(tsce);
1064                         }
1065                     } else {
1066 
1067                         // we are after the range covered by the annual pole file,
1068                         // we use the polynomial extension
1069                         final double t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
1070                         meanPoleX = xp0 + t * xp0Dot;
1071                         meanPoleY = yp0 + t * yp0Dot;
1072 
1073                     }
1074 
1075                     // evaluate wobble variables
1076                     final PoleCorrection correction = eopHistory.getPoleCorrection(date);
1077                     final double m1 = correction.getXp() - meanPoleX;
1078                     final double m2 = meanPoleY - correction.getYp();
1079 
1080                     return new double[] {
1081                         // the following correspond to the equations published in IERS 2003 conventions,
1082                         // section 6.2 page 65. In the publication, the equations read:
1083                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
1084                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
1085                         // However, it seems there are sign errors in these equations, which have
1086                         // been fixed in IERS 2010 conventions, section 6.4 page 94. In these newer
1087                         // publication, the equations read:
1088                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
1089                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
1090                         // the newer equations seem more consistent with the premises as the
1091                         // deformation due to the centrifugal potential has the form:
1092                         // −Ω²r²/2 sin 2θ Re [k₂(m₁ − im₂) exp(iλ)] where k₂ is the complex
1093                         // number 0.3077 + 0.0036i, so the real part in the previous equation is:
1094                         // A[Re(k₂) m₁ + Im(k₂) m₂)] cos λ + A[Re(k₂) m₂ - Im(k₂) m₁] sin λ
1095                         // identifying this with ∆C₂₁ cos λ + ∆S₂₁ sin λ we get:
1096                         // ∆C₂₁ = A Re(k₂) [m₁ + Im(k₂)/Re(k₂) m₂)]
1097                         // ∆S₂₁ = A Re(k₂) [m₂ - Im(k₂)/Re(k₂) m₁)]
1098                         // and Im(k₂)/Re(k₂) is very close to +0.00115
1099                         // As the equation as written in the IERS 2003 conventions are used in
1100                         // legacy systems, we have reproduced this alleged error here (and fixed it in
1101                         // the IERS 2010 conventions below) for validation purposes. We don't recommend
1102                         // using the IERS 2003 conventions for solid pole tide computation other than
1103                         // for validation or reproducibility of legacy applications behavior.
1104                         // As solid pole tide is small and as the sign change is on the smallest coefficient,
1105                         // the effect is quite small. A test case on a propagated orbit showed a position change
1106                         // slightly below 0.4m after a 30 days propagation on a Low Earth Orbit
1107                         globalFactor * (m1 - ratio * m2),
1108                         globalFactor * (m2 + ratio * m1),
1109                     };
1110 
1111                 }
1112             };
1113 
1114         }
1115 
1116         /** {@inheritDoc} */
1117         @Override
1118         public TimeFunction<double[]> getOceanPoleTide(final EOPHistory eopHistory)
1119             throws OrekitException {
1120 
1121             return new TimeFunction<double[]>() {
1122                 /** {@inheritDoc} */
1123                 @Override
1124                 public double[] value(final AbsoluteDate date) {
1125                     // there are no model for ocean pole tide prior to conventions 2010
1126                     return new double[] {
1127                         0.0, 0.0
1128                     };
1129                 }
1130             };
1131         }
1132 
1133     },
1134 
1135     /** Constant for IERS 2010 conventions. */
1136     IERS_2010 {
1137 
1138         /** Nutation arguments resources. */
1139         private static final String NUTATION_ARGUMENTS = IERS_BASE + "2010/nutation-arguments.txt";
1140 
1141         /** X series resources. */
1142         private static final String X_SERIES           = IERS_BASE + "2010/tab5.2a.txt";
1143 
1144         /** Y series resources. */
1145         private static final String Y_SERIES           = IERS_BASE + "2010/tab5.2b.txt";
1146 
1147         /** S series resources. */
1148         private static final String S_SERIES           = IERS_BASE + "2010/tab5.2d.txt";
1149 
1150         /** Psi series resources. */
1151         private static final String PSI_SERIES         = IERS_BASE + "2010/tab5.3a.txt";
1152 
1153         /** Epsilon series resources. */
1154         private static final String EPSILON_SERIES     = IERS_BASE + "2010/tab5.3b.txt";
1155 
1156         /** Greenwhich sidereal time series resources. */
1157         private static final String GST_SERIES         = IERS_BASE + "2010/tab5.2e.txt";
1158 
1159         /** Tidal correction for xp, yp series resources. */
1160         private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "2010/tab8.2ab.txt";
1161 
1162         /** Tidal correction for UT1 resources. */
1163         private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "2010/tab8.3ab.txt";
1164 
1165         /** Love numbers resources. */
1166         private static final String LOVE_NUMBERS = IERS_BASE + "2010/tab6.3.txt";
1167 
1168         /** Frequency dependence model for k₂₀. */
1169         private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5b.txt";
1170 
1171         /** Frequency dependence model for k₂₁. */
1172         private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5a.txt";
1173 
1174         /** Frequency dependence model for k₂₂. */
1175         private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5c.txt";
1176 
1177         /** {@inheritDoc} */
1178         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
1179             throws OrekitException {
1180             return new FundamentalNutationArguments(this, timeScale,
1181                                                     getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
1182         }
1183 
1184         /** {@inheritDoc} */
1185         @Override
1186         public TimeFunction<Double> getMeanObliquityFunction() throws OrekitException {
1187 
1188             // epsilon 0 value from chapter 5, page 56, other terms from equation 5.40 page 65
1189             final PolynomialNutation<DerivativeStructure> epsilonA =
1190                     new PolynomialNutation<DerivativeStructure>(84381.406        * Constants.ARC_SECONDS_TO_RADIANS,
1191                                                                   -46.836769     * Constants.ARC_SECONDS_TO_RADIANS,
1192                                                                    -0.0001831    * Constants.ARC_SECONDS_TO_RADIANS,
1193                                                                     0.00200340   * Constants.ARC_SECONDS_TO_RADIANS,
1194                                                                    -0.000000576  * Constants.ARC_SECONDS_TO_RADIANS,
1195                                                                    -0.0000000434 * Constants.ARC_SECONDS_TO_RADIANS);
1196 
1197             return new TimeFunction<Double>() {
1198 
1199                 /** {@inheritDoc} */
1200                 @Override
1201                 public Double value(final AbsoluteDate date) {
1202                     return epsilonA.value(evaluateTC(date));
1203                 }
1204 
1205             };
1206 
1207         }
1208 
1209         /** {@inheritDoc} */
1210         @Override
1211         public TimeFunction<double[]> getXYSpXY2Function() throws OrekitException {
1212 
1213             // set up nutation arguments
1214             final FundamentalNutationArguments arguments = getNutationArguments(null);
1215 
1216             // set up Poisson series
1217             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1218             final PoissonSeriesParser<DerivativeStructure> parser =
1219                     new PoissonSeriesParser<DerivativeStructure>(17).
1220                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
1221                         withFirstDelaunay(4).
1222                         withFirstPlanetary(9).
1223                         withSinCos(0, 2, microAS, 3, microAS);
1224             final PoissonSeries<DerivativeStructure> xSeries = parser.parse(getStream(X_SERIES), X_SERIES);
1225             final PoissonSeries<DerivativeStructure> ySeries = parser.parse(getStream(Y_SERIES), Y_SERIES);
1226             final PoissonSeries<DerivativeStructure> sSeries = parser.parse(getStream(S_SERIES), S_SERIES);
1227             @SuppressWarnings("unchecked")
1228             final PoissonSeries.CompiledSeries<DerivativeStructure> xys = PoissonSeries.compile(xSeries, ySeries, sSeries);
1229 
1230             // create a function evaluating the series
1231             return new TimeFunction<double[]>() {
1232 
1233                 /** {@inheritDoc} */
1234                 @Override
1235                 public double[] value(final AbsoluteDate date) {
1236                     return xys.value(arguments.evaluateAll(date));
1237                 }
1238 
1239             };
1240 
1241         }
1242 
1243         /** {@inheritDoc} */
1244         public LoveNumbers getLoveNumbers() throws OrekitException {
1245             return loadLoveNumbers(LOVE_NUMBERS);
1246         }
1247 
1248         /** {@inheritDoc} */
1249         public TimeFunction<double[]> getTideFrequencyDependenceFunction(final TimeScale ut1)
1250             throws OrekitException {
1251 
1252             // set up nutation arguments
1253             final FundamentalNutationArguments arguments = getNutationArguments(ut1);
1254 
1255             // set up Poisson series
1256             final PoissonSeriesParser<DerivativeStructure> k20Parser =
1257                     new PoissonSeriesParser<DerivativeStructure>(18).
1258                         withOptionalColumn(1).
1259                         withDoodson(4, 2).
1260                         withFirstDelaunay(10);
1261             final PoissonSeriesParser<DerivativeStructure> k21Parser =
1262                     new PoissonSeriesParser<DerivativeStructure>(18).
1263                         withOptionalColumn(1).
1264                         withDoodson(4, 3).
1265                         withFirstDelaunay(10);
1266             final PoissonSeriesParser<DerivativeStructure> k22Parser =
1267                     new PoissonSeriesParser<DerivativeStructure>(16).
1268                         withOptionalColumn(1).
1269                         withDoodson(4, 2).
1270                         withFirstDelaunay(10);
1271 
1272             final double pico = 1.0e-12;
1273             final PoissonSeries<DerivativeStructure> c20Series =
1274                     k20Parser.
1275                     // TODO: check sign!
1276 //                    withSinCos(0, 18, pico, 16, pico).
1277                     withSinCos(0, 18, -pico, 16, pico).
1278                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
1279             final PoissonSeries<DerivativeStructure> c21Series =
1280                     k21Parser.
1281                     withSinCos(0, 17, pico, 18, pico).
1282                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
1283             final PoissonSeries<DerivativeStructure> s21Series =
1284                     k21Parser.
1285                     withSinCos(0, 18, -pico, 17, pico).
1286                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
1287             final PoissonSeries<DerivativeStructure> c22Series =
1288                     k22Parser.
1289                     withSinCos(0, -1, pico, 16, pico).
1290                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
1291             final PoissonSeries<DerivativeStructure> s22Series =
1292                     k22Parser.
1293                     withSinCos(0, 16, -pico, -1, pico).
1294                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
1295 
1296             @SuppressWarnings("unchecked")
1297             final PoissonSeries.CompiledSeries<DerivativeStructure> kSeries =
1298                 PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);
1299 
1300             return new TimeFunction<double[]>() {
1301                 /** {@inheritDoc} */
1302                 @Override
1303                 public double[] value(final AbsoluteDate date) {
1304                     return kSeries.value(arguments.evaluateAll(date));
1305                 }
1306             };
1307 
1308         }
1309 
1310         /** {@inheritDoc} */
1311         @Override
1312         public double getPermanentTide() throws OrekitException {
1313             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
1314         }
1315 
1316         /** Compute pole wobble variables m₁ and m₂.
1317          * @param date current date
1318          * @param eopHistory EOP history
1319          * @return array containing m₁ and m₂
1320          */
1321         private double[] computePoleWobble(final AbsoluteDate date, final EOPHistory eopHistory) {
1322 
1323             // polynomial model from IERS 2010, table 7.7
1324             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
1325             final double f1 = f0 / Constants.JULIAN_YEAR;
1326             final double f2 = f1 / Constants.JULIAN_YEAR;
1327             final double f3 = f2 / Constants.JULIAN_YEAR;
1328             final AbsoluteDate changeDate = new AbsoluteDate(2010, 1, 1, TimeScalesFactory.getTT());
1329 
1330             // evaluate mean pole
1331             final double[] xPolynomial;
1332             final double[] yPolynomial;
1333             if (date.compareTo(changeDate) <= 0) {
1334                 xPolynomial = new double[] {
1335                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
1336                 };
1337                 yPolynomial = new double[] {
1338                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
1339                 };
1340             } else {
1341                 xPolynomial = new double[] {
1342                     23.513 * f0, 7.6141 * f1
1343                 };
1344                 yPolynomial = new double[] {
1345                     358.891 * f0,  -0.6287 * f1
1346                 };
1347             }
1348             double meanPoleX = 0;
1349             double meanPoleY = 0;
1350             final double t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
1351             for (int i = xPolynomial.length - 1; i >= 0; --i) {
1352                 meanPoleX = meanPoleX * t + xPolynomial[i];
1353             }
1354             for (int i = yPolynomial.length - 1; i >= 0; --i) {
1355                 meanPoleY = meanPoleY * t + yPolynomial[i];
1356             }
1357 
1358             // evaluate wobble variables
1359             final PoleCorrection correction = eopHistory.getPoleCorrection(date);
1360             final double m1 = correction.getXp() - meanPoleX;
1361             final double m2 = meanPoleY - correction.getYp();
1362 
1363             return new double[] {
1364                 m1, m2
1365             };
1366 
1367         }
1368 
1369         /** {@inheritDoc} */
1370         @Override
1371         public TimeFunction<double[]> getSolidPoleTide(final EOPHistory eopHistory)
1372             throws OrekitException {
1373 
1374             // constants from IERS 2010, section 6.4
1375             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
1376             final double ratio        =  0.00115;
1377 
1378             return new TimeFunction<double[]>() {
1379                 /** {@inheritDoc} */
1380                 @Override
1381                 public double[] value(final AbsoluteDate date) {
1382 
1383                     // evaluate wobble variables
1384                     final double[] wobbleM = computePoleWobble(date, eopHistory);
1385 
1386                     return new double[] {
1387                         // the following correspond to the equations published in IERS 2010 conventions,
1388                         // section 6.4 page 94. The equations read:
1389                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
1390                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
1391                         // These equations seem to fix what was probably a sign error in IERS 2003
1392                         // conventions section 6.2 page 65. In this older publication, the equations read:
1393                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
1394                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
1395                         globalFactor * (wobbleM[0] + ratio * wobbleM[1]),
1396                         globalFactor * (wobbleM[1] - ratio * wobbleM[0])
1397                     };
1398 
1399                 }
1400             };
1401 
1402         }
1403 
1404         /** {@inheritDoc} */
1405         @Override
1406         public TimeFunction<double[]> getOceanPoleTide(final EOPHistory eopHistory)
1407             throws OrekitException {
1408 
1409             return new TimeFunction<double[]>() {
1410                 /** {@inheritDoc} */
1411                 @Override
1412                 public double[] value(final AbsoluteDate date) {
1413 
1414                     // evaluate wobble variables
1415                     final double[] wobbleM = computePoleWobble(date, eopHistory);
1416 
1417                     return new double[] {
1418                         // the following correspond to the equations published in IERS 2010 conventions,
1419                         // section 6.4 page 94 equation 6.24:
1420                         // ∆C₂₁ = −2.1778 × 10⁻¹⁰ (m₁ − 0.01724m₂)
1421                         // ∆S₂₁ = −1.7232 × 10⁻¹⁰ (m₂ − 0.03365m₁)
1422                         -2.1778e-10 * (wobbleM[0] - 0.01724 * wobbleM[1]) / Constants.ARC_SECONDS_TO_RADIANS,
1423                         -1.7232e-10 * (wobbleM[1] - 0.03365 * wobbleM[0]) / Constants.ARC_SECONDS_TO_RADIANS
1424                     };
1425 
1426                 }
1427             };
1428 
1429         }
1430 
1431         /** {@inheritDoc} */
1432         @Override
1433         public TimeFunction<double[]> getPrecessionFunction() throws OrekitException {
1434 
1435             // set up the conventional polynomials
1436             // the following values are from equation 5.40 in IERS 2010 conventions
1437             final PolynomialNutation<DerivativeStructure> psiA =
1438                     new PolynomialNutation<DerivativeStructure>(   0.0,
1439                                                                 5038.481507     * Constants.ARC_SECONDS_TO_RADIANS,
1440                                                                   -1.0790069    * Constants.ARC_SECONDS_TO_RADIANS,
1441                                                                   -0.00114045   * Constants.ARC_SECONDS_TO_RADIANS,
1442                                                                    0.000132851  * Constants.ARC_SECONDS_TO_RADIANS,
1443                                                                   -0.0000000951 * Constants.ARC_SECONDS_TO_RADIANS);
1444             final PolynomialNutation<DerivativeStructure> omegaA =
1445                     new PolynomialNutation<DerivativeStructure>(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
1446                                                                 -0.025754     * Constants.ARC_SECONDS_TO_RADIANS,
1447                                                                  0.0512623    * Constants.ARC_SECONDS_TO_RADIANS,
1448                                                                 -0.00772503   * Constants.ARC_SECONDS_TO_RADIANS,
1449                                                                 -0.000000467  * Constants.ARC_SECONDS_TO_RADIANS,
1450                                                                  0.0000003337 * Constants.ARC_SECONDS_TO_RADIANS);
1451             final PolynomialNutation<DerivativeStructure> chiA =
1452                     new PolynomialNutation<DerivativeStructure>( 0.0,
1453                                                                 10.556403     * Constants.ARC_SECONDS_TO_RADIANS,
1454                                                                 -2.3814292    * Constants.ARC_SECONDS_TO_RADIANS,
1455                                                                 -0.00121197   * Constants.ARC_SECONDS_TO_RADIANS,
1456                                                                  0.000170663  * Constants.ARC_SECONDS_TO_RADIANS,
1457                                                                 -0.0000000560 * Constants.ARC_SECONDS_TO_RADIANS);
1458 
1459             return new TimeFunction<double[]>() {
1460                 /** {@inheritDoc} */
1461                 @Override
1462                 public double[] value(final AbsoluteDate date) {
1463                     final double tc = evaluateTC(date);
1464                     return new double[] {
1465                         psiA.value(tc), omegaA.value(tc), chiA.value(tc)
1466                     };
1467                 }
1468             };
1469 
1470         }
1471 
1472          /** {@inheritDoc} */
1473         @Override
1474         public TimeFunction<double[]> getNutationFunction()
1475             throws OrekitException {
1476 
1477             // set up nutation arguments
1478             final FundamentalNutationArguments arguments = getNutationArguments(null);
1479 
1480             // set up Poisson series
1481             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1482             final PoissonSeriesParser<DerivativeStructure> parser =
1483                     new PoissonSeriesParser<DerivativeStructure>(17).
1484                         withFirstDelaunay(4).
1485                         withFirstPlanetary(9).
1486                         withSinCos(0, 2, microAS, 3, microAS);
1487             final PoissonSeries<DerivativeStructure> psiSeries     = parser.parse(getStream(PSI_SERIES), PSI_SERIES);
1488             final PoissonSeries<DerivativeStructure> epsilonSeries = parser.parse(getStream(EPSILON_SERIES), EPSILON_SERIES);
1489             @SuppressWarnings("unchecked")
1490             final PoissonSeries.CompiledSeries<DerivativeStructure> psiEpsilonSeries =
1491                     PoissonSeries.compile(psiSeries, epsilonSeries);
1492 
1493             return new TimeFunction<double[]>() {
1494                 /** {@inheritDoc} */
1495                 @Override
1496                 public double[] value(final AbsoluteDate date) {
1497                     final BodiesElements elements = arguments.evaluateAll(date);
1498                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
1499                     return new double[] {
1500                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements)
1501                     };
1502                 }
1503             };
1504 
1505         }
1506 
1507         /** {@inheritDoc} */
1508         @Override
1509         public TimeFunction<DerivativeStructure> getGMSTFunction(final TimeScale ut1) throws OrekitException {
1510 
1511             // Earth Rotation Angle
1512             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);
1513 
1514             // Polynomial part of the apparent sidereal time series
1515             // which is the opposite of Equation of Origins (EO)
1516             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1517             final PoissonSeriesParser<DerivativeStructure> parser =
1518                     new PoissonSeriesParser<DerivativeStructure>(17).
1519                         withFirstDelaunay(4).
1520                         withFirstPlanetary(9).
1521                         withSinCos(0, 2, microAS, 3, microAS).
1522                         withPolynomialPart('t', Unit.ARC_SECONDS);
1523             final PolynomialNutation<DerivativeStructure> minusEO =
1524                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();
1525 
1526             // create a function evaluating the series
1527             return new TimeFunction<DerivativeStructure>() {
1528 
1529                 /** {@inheritDoc} */
1530                 @Override
1531                 public DerivativeStructure value(final AbsoluteDate date) {
1532                     return era.value(date).add(minusEO.value(dsEvaluateTC(date)));
1533                 }
1534 
1535             };
1536 
1537         }
1538 
1539         /** {@inheritDoc} */
1540         @Override
1541         public TimeFunction<DerivativeStructure> getGASTFunction(final TimeScale ut1,
1542                                                                  final EOPHistory eopHistory)
1543             throws OrekitException {
1544 
1545             // set up nutation arguments
1546             final FundamentalNutationArguments arguments = getNutationArguments(null);
1547 
1548             // mean obliquity function
1549             final TimeFunction<Double> epsilon = getMeanObliquityFunction();
1550 
1551             // set up Poisson series
1552             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1553             final PoissonSeriesParser<DerivativeStructure> baseParser =
1554                     new PoissonSeriesParser<DerivativeStructure>(17).
1555                         withFirstDelaunay(4).
1556                         withFirstPlanetary(9).
1557                         withSinCos(0, 2, microAS, 3, microAS);
1558             final PoissonSeriesParser<DerivativeStructure> gstParser  = baseParser.withPolynomialPart('t', Unit.ARC_SECONDS);
1559             final PoissonSeries<DerivativeStructure> psiSeries        = baseParser.parse(getStream(PSI_SERIES), PSI_SERIES);
1560             final PoissonSeries<DerivativeStructure> gstSeries        = gstParser.parse(getStream(GST_SERIES), GST_SERIES);
1561             @SuppressWarnings("unchecked")
1562             final PoissonSeries.CompiledSeries<DerivativeStructure> psiGstSeries =
1563                     PoissonSeries.compile(psiSeries, gstSeries);
1564 
1565             // ERA function
1566             final TimeFunction<DerivativeStructure> era = getEarthOrientationAngleFunction(ut1);
1567 
1568             return new TimeFunction<DerivativeStructure>() {
1569 
1570                 /** {@inheritDoc} */
1571                 @Override
1572                 public DerivativeStructure value(final AbsoluteDate date) {
1573 
1574                     // evaluate equation of origins
1575                     final BodiesElements elements = arguments.evaluateAll(date);
1576                     final double[] angles = psiGstSeries.value(elements);
1577                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
1578                     final double deltaPsi = angles[0] + ddPsi;
1579                     final double epsilonA = epsilon.value(date);
1580 
1581                     // subtract equation of origin from EA
1582                     // (hence add the series above which have the sign included)
1583                     return era.value(date).add(deltaPsi * FastMath.cos(epsilonA) + angles[1]);
1584 
1585                 }
1586 
1587             };
1588 
1589         }
1590 
1591         /** {@inheritDoc} */
1592         @Override
1593         public TimeFunction<double[]> getEOPTidalCorrection()
1594             throws OrekitException {
1595 
1596             // set up nutation arguments
1597             // BEWARE! Using TT as the time scale here and not UT1 is intentional!
1598             // as this correction is used to compute UT1 itself, it is not surprising we cannot use UT1 yet,
1599             // however, using the close UTC as would seem logical make the comparison with interp.f from IERS fail
1600             // looking in the interp.f code, the same TT scale is used for both Delaunay and gamma argument
1601             final FundamentalNutationArguments arguments = getNutationArguments(TimeScalesFactory.getTT());
1602 
1603             // set up Poisson series
1604             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
1605             final PoissonSeriesParser<DerivativeStructure> xyParser = new PoissonSeriesParser<DerivativeStructure>(13).
1606                     withOptionalColumn(1).
1607                     withGamma(2).
1608                     withFirstDelaunay(3);
1609             final PoissonSeries<DerivativeStructure> xSeries =
1610                     xyParser.
1611                     withSinCos(0, 10, microAS, 11, microAS).
1612                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
1613             final PoissonSeries<DerivativeStructure> ySeries =
1614                     xyParser.
1615                     withSinCos(0, 12, microAS, 13, microAS).
1616                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
1617 
1618             final double microS = 1.0e-6;
1619             final PoissonSeriesParser<DerivativeStructure> ut1Parser = new PoissonSeriesParser<DerivativeStructure>(11).
1620                     withOptionalColumn(1).
1621                     withGamma(2).
1622                     withFirstDelaunay(3).
1623                     withSinCos(0, 10, microS, 11, microS);
1624             final PoissonSeries<DerivativeStructure> ut1Series =
1625                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);
1626 
1627             @SuppressWarnings("unchecked")
1628             final PoissonSeries.CompiledSeries<DerivativeStructure> correctionSeries =
1629                     PoissonSeries.compile(xSeries, ySeries, ut1Series);
1630 
1631             return new TimeFunction<double[]>() {
1632                 /** {@inheritDoc} */
1633                 @Override
1634                 public double[] value(final AbsoluteDate date) {
1635                     final FieldBodiesElements<DerivativeStructure> elements =
1636                             arguments.evaluateDerivative(date);
1637                     final DerivativeStructure[] correction = correctionSeries.value(elements);
1638                     return new double[] {
1639                         correction[0].getValue(),
1640                         correction[1].getValue(),
1641                         correction[2].getValue(),
1642                         -correction[2].getPartialDerivative(1) * Constants.JULIAN_DAY
1643                     };
1644                 }
1645             };
1646 
1647         }
1648 
1649     };
1650 
1651     /** IERS conventions resources base directory. */
1652     private static final String IERS_BASE = "/assets/org/orekit/IERS-conventions/";
1653 
1654     /** Get the reference epoch for fundamental nutation arguments.
1655      * @return reference epoch for fundamental nutation arguments
1656      * @since 6.1
1657      */
1658     public AbsoluteDate getNutationReferenceEpoch() {
1659         // IERS 1996, IERS 2003 and IERS 2010 use the same J2000.0 reference date
1660         return AbsoluteDate.J2000_EPOCH;
1661     }
1662 
1663     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
1664      * @param date current date
1665      * @return date offset in Julian centuries
1666      * @since 6.1
1667      */
1668     public double evaluateTC(final AbsoluteDate date) {
1669         return date.durationFrom(getNutationReferenceEpoch()) / Constants.JULIAN_CENTURY;
1670     }
1671 
1672     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
1673      * @param date current date
1674      * @return date offset in Julian centuries
1675      * @since 6.1
1676      */
1677     public DerivativeStructure dsEvaluateTC(final AbsoluteDate date) {
1678         return new DerivativeStructure(1, 1, evaluateTC(date), 1.0 / Constants.JULIAN_CENTURY);
1679     }
1680 
1681     /** Get the fundamental nutation arguments.
1682      * @param timeScale time scale for computing Greenwich Mean Sidereal Time
1683      * (typically {@link TimeScalesFactory#getUT1(IERSConventions, boolean) UT1})
1684      * @return fundamental nutation arguments
1685      * @exception OrekitException if fundamental nutation arguments cannot be loaded
1686      * @since 6.1
1687      */
1688     public abstract FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
1689         throws OrekitException;
1690 
1691     /** Get the function computing mean obliquity of the ecliptic.
1692      * @return function computing mean obliquity of the ecliptic
1693      * @exception OrekitException if table cannot be loaded
1694      * @since 6.1
1695      */
1696     public abstract TimeFunction<Double> getMeanObliquityFunction() throws OrekitException;
1697 
1698     /** Get the function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components.
1699      * <p>
1700      * The returned function computes the two X, Y components of CIP and the S+XY/2 component of the non-rotating CIO.
1701      * </p>
1702      * @return function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components
1703      * @exception OrekitException if table cannot be loaded
1704      * @since 6.1
1705      */
1706     public abstract TimeFunction<double[]> getXYSpXY2Function()
1707         throws OrekitException;
1708 
1709     /** Get the function computing the raw Earth Orientation Angle.
1710      * <p>
1711      * The raw angle does not contain any correction. If for example dTU1 correction
1712      * due to tidal effect is desired, it must be added afterward by the caller.
1713      * The returned value contain the angle as the value and the angular rate as
1714      * the first derivative.
1715      * </p>
1716      * @param ut1 UT1 time scale
1717      * @return function computing the rawEarth Orientation Angle, in the non-rotating origin paradigm,
1718      * the return value containing both the angle and its first time derivative
1719      * @since 6.1
1720      */
1721     public TimeFunction<DerivativeStructure> getEarthOrientationAngleFunction(final TimeScale ut1) {
1722         return new StellarAngleCapitaine(ut1);
1723     }
1724 
1725 
1726     /** Get the function computing the precession angles.
1727      * <p>
1728      * The function returned computes the three precession angles
1729      * ψ<sub>A</sub> (around Z axis), ω<sub>A</sub> (around X axis)
1730      * and χ<sub>A</sub> (around Z axis). The constant angle ε₀
1731      * for the fourth rotation (around X axis) can be retrieved by evaluating the
1732      * function returned by {@link #getMeanObliquityFunction()} at {@link
1733      * #getNutationReferenceEpoch() nutation reference epoch}.
1734      * </p>
1735      * @return function computing the precession angle
1736      * @exception OrekitException if table cannot be loaded
1737      * @since 6.1
1738      */
1739     public abstract TimeFunction<double[]> getPrecessionFunction() throws OrekitException;
1740 
1741     /** Get the function computing the nutation angles.
1742      * <p>
1743      * The function returned computes the two classical angles ΔΨ and Δε,
1744      * and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
1745      * resolution C7 (the correction is forced to 0 before this date)
1746      * </p>
1747      * @return function computing the nutation in longitude ΔΨ and Δε
1748      * and the correction of equation of equinoxes
1749      * @exception OrekitException if table cannot be loaded
1750      * @since 6.1
1751      */
1752     public abstract TimeFunction<double[]> getNutationFunction()
1753         throws OrekitException;
1754 
1755     /** Get the function computing Greenwich mean sidereal time, in radians.
1756      * @param ut1 UT1 time scale
1757      * @return function computing Greenwich mean sidereal time,
1758      * the return value containing both the angle and its first time derivative
1759      * @exception OrekitException if table cannot be loaded
1760      * @since 6.1
1761      */
1762     public abstract TimeFunction<DerivativeStructure> getGMSTFunction(TimeScale ut1)
1763         throws OrekitException;
1764 
1765     /** Get the function computing Greenwich apparent sidereal time, in radians.
1766      * @param ut1 UT1 time scale
1767      * @param eopHistory EOP history
1768      * @return function computing Greenwich apparent sidereal time,
1769      * the return value containing both the angle and its first time derivative
1770      * @exception OrekitException if table cannot be loaded
1771      * @since 6.1
1772      */
1773     public abstract TimeFunction<DerivativeStructure> getGASTFunction(TimeScale ut1,
1774                                                                       EOPHistory eopHistory)
1775         throws OrekitException;
1776 
1777     /** Get the function computing tidal corrections for Earth Orientation Parameters.
1778      * @return function computing tidal corrections for Earth Orientation Parameters,
1779      * for xp, yp, ut1 and lod respectively
1780      * @exception OrekitException if table cannot be loaded
1781      * @since 6.1
1782      */
1783     public abstract TimeFunction<double[]> getEOPTidalCorrection()
1784         throws OrekitException;
1785 
1786     /** Get the Love numbers.
1787      * @return Love numbers
1788      * @exception OrekitException if table cannot be loaded
1789      * @since 6.1
1790      */
1791     public abstract LoveNumbers getLoveNumbers()
1792         throws OrekitException;
1793 
1794     /** Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
1795      * @param ut1 UT1 time scale
1796      * @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
1797      * @exception OrekitException if table cannot be loaded
1798      * @since 6.1
1799      */
1800     public abstract TimeFunction<double[]> getTideFrequencyDependenceFunction(TimeScale ut1)
1801         throws OrekitException;
1802 
1803     /** Get the permanent tide to be <em>removed</em> from ΔC₂₀ when zero-tide potentials are used.
1804      * @return permanent tide to remove
1805      * @exception OrekitException if table cannot be loaded
1806      */
1807     public abstract double getPermanentTide() throws OrekitException;
1808 
1809     /** Get the function computing solid pole tide (ΔC₂₁, ΔS₂₁).
1810      * @param eopHistory EOP history
1811      * @return model for solid pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
1812      * @exception OrekitException if table cannot be loaded
1813      * @since 6.1
1814      */
1815     public abstract TimeFunction<double[]> getSolidPoleTide(EOPHistory eopHistory)
1816         throws OrekitException;
1817 
1818     /** Get the function computing ocean pole tide (ΔC₂₁, ΔS₂₁).
1819      * @param eopHistory EOP history
1820      * @return model for ocean pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
1821      * @exception OrekitException if table cannot be loaded
1822      * @since 6.1
1823      */
1824     public abstract TimeFunction<double[]> getOceanPoleTide(EOPHistory eopHistory)
1825         throws OrekitException;
1826 
1827     /** Interface for functions converting nutation corrections between
1828      * δΔψ/δΔε to δX/δY.
1829      * <ul>
1830      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
1831      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
1832      * </ul>
1833      * @since 6.1
1834      */
1835     public interface NutationCorrectionConverter {
1836 
1837         /** Convert nutation corrections.
1838          * @param date current date
1839          * @param ddPsi δΔψ part of the nutation correction
1840          * @param ddEpsilon δΔε part of the nutation correction
1841          * @return array containing δX and δY
1842          * @exception OrekitException if correction cannot be converted
1843          */
1844         double[] toNonRotating(AbsoluteDate date, double ddPsi, double ddEpsilon)
1845             throws OrekitException;
1846 
1847         /** Convert nutation corrections.
1848          * @param date current date
1849          * @param dX δX part of the nutation correction
1850          * @param dY δY part of the nutation correction
1851          * @return array containing δΔψ and δΔε
1852          * @exception OrekitException if correction cannot be converted
1853          */
1854         double[] toEquinox(AbsoluteDate date, double dX, double dY)
1855             throws OrekitException;
1856 
1857     }
1858 
1859     /** Create a function converting nutation corrections between
1860      * δX/δY and δΔψ/δΔε.
1861      * <ul>
1862      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
1863      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
1864      * </ul>
1865      * @return a new converter
1866      * @exception OrekitException if some convention table cannot be loaded
1867      * @since 6.1
1868      */
1869     public NutationCorrectionConverter getNutationCorrectionConverter()
1870         throws OrekitException {
1871 
1872         // get models parameters
1873         final TimeFunction<double[]> precessionFunction = getPrecessionFunction();
1874         final TimeFunction<Double> epsilonAFunction = getMeanObliquityFunction();
1875         final AbsoluteDate date0 = getNutationReferenceEpoch();
1876         final double cosE0 = FastMath.cos(epsilonAFunction.value(date0));
1877 
1878         return new NutationCorrectionConverter() {
1879 
1880             /** {@inheritDoc} */
1881             @Override
1882             public double[] toNonRotating(final AbsoluteDate date,
1883                                           final double ddPsi, final double ddEpsilon)
1884                 throws OrekitException {
1885                 // compute precession angles psiA, omegaA and chiA
1886                 final double[] angles = precessionFunction.value(date);
1887 
1888                 // conversion coefficients
1889                 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
1890                 final double c     = angles[0] * cosE0 - angles[2];
1891 
1892                 // convert nutation corrections (equation 23/IERS-2003 or 5.25/IERS-2010)
1893                 return new double[] {
1894                     sinEA * ddPsi + c * ddEpsilon,
1895                     ddEpsilon - c * sinEA * ddPsi
1896                 };
1897 
1898             }
1899 
1900             /** {@inheritDoc} */
1901             @Override
1902             public double[] toEquinox(final AbsoluteDate date,
1903                                       final double dX, final double dY)
1904                 throws OrekitException {
1905                 // compute precession angles psiA, omegaA and chiA
1906                 final double[] angles   = precessionFunction.value(date);
1907 
1908                 // conversion coefficients
1909                 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
1910                 final double c     = angles[0] * cosE0 - angles[2];
1911                 final double opC2  = 1 + c * c;
1912 
1913                 // convert nutation corrections (inverse of equation 23/IERS-2003 or 5.25/IERS-2010)
1914                 return new double[] {
1915                     (dX - c * dY) / (sinEA * opC2),
1916                     (dY + c * dX) / opC2
1917                 };
1918             }
1919 
1920         };
1921 
1922     }
1923 
1924     /** Load the Love numbers.
1925      * @param nameLove name of the Love number resource
1926      * @return Love numbers
1927      * @exception OrekitException if the Love numbers embedded in the
1928      * library cannot be read
1929      */
1930     protected LoveNumbers loadLoveNumbers(final String nameLove) throws OrekitException {
1931         InputStream stream = null;
1932         BufferedReader reader = null;
1933         try {
1934 
1935             // allocate the three triangular arrays for real, imaginary and time-dependent numbers
1936             final double[][] real      = new double[4][];
1937             final double[][] imaginary = new double[4][];
1938             final double[][] plus      = new double[4][];
1939             for (int i = 0; i < real.length; ++i) {
1940                 real[i]      = new double[i + 1];
1941                 imaginary[i] = new double[i + 1];
1942                 plus[i]      = new double[i + 1];
1943             }
1944 
1945             stream = IERSConventions.class.getResourceAsStream(nameLove);
1946             if (stream == null) {
1947                 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, nameLove);
1948             }
1949 
1950             // setup the reader
1951             reader = new BufferedReader(new InputStreamReader(stream, "UTF-8"));
1952 
1953             String line = reader.readLine();
1954             int lineNumber = 1;
1955 
1956             // look for the Love numbers
1957             while (line != null) {
1958 
1959                 line = line.trim();
1960                 if (!(line.isEmpty() || line.startsWith("#"))) {
1961                     final String[] fields = line.split("\\p{Space}+");
1962                     if (fields.length != 5) {
1963                         // this should never happen with files embedded within Orekit
1964                         throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
1965                                                   lineNumber, nameLove, line);
1966                     }
1967                     final int n = Integer.parseInt(fields[0]);
1968                     final int m = Integer.parseInt(fields[1]);
1969                     if ((n < 2) || (n > 3) || (m < 0) || (m > n)) {
1970                         // this should never happen with files embedded within Orekit
1971                         throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
1972                                                   lineNumber, nameLove, line);
1973 
1974                     }
1975                     real[n][m]      = Double.parseDouble(fields[2]);
1976                     imaginary[n][m] = Double.parseDouble(fields[3]);
1977                     plus[n][m]      = Double.parseDouble(fields[4]);
1978                 }
1979 
1980                 // next line
1981                 lineNumber++;
1982                 line = reader.readLine();
1983 
1984             }
1985 
1986             return new LoveNumbers(real, imaginary, plus);
1987 
1988         } catch (IOException ioe) {
1989             throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, nameLove);
1990         } finally {
1991             try {
1992                 if (stream != null) {
1993                     stream.close();
1994                 }
1995                 if (reader != null) {
1996                     reader.close();
1997                 }
1998             } catch (IOException ioe) {
1999                 // ignored here
2000             }
2001         }
2002     }
2003 
2004     /** Get a data stream.
2005      * @param name file name of the resource stream
2006      * @return stream
2007      */
2008     private static InputStream getStream(final String name) {
2009         return IERSConventions.class.getResourceAsStream(name);
2010     }
2011 
2012     /** Correction to equation of equinoxes.
2013      * <p>IAU 1994 resolution C7 added two terms to the equation of equinoxes
2014      * taking effect since 1997-02-27 for continuity
2015      * </p>
2016      */
2017     private static class IAU1994ResolutionC7 {
2018 
2019         /** First Moon correction term for the Equation of the Equinoxes. */
2020         private static final double EQE1 =     0.00264  * Constants.ARC_SECONDS_TO_RADIANS;
2021 
2022         /** Second Moon correction term for the Equation of the Equinoxes. */
2023         private static final double EQE2 =     0.000063 * Constants.ARC_SECONDS_TO_RADIANS;
2024 
2025         /** Start date for applying Moon corrections to the equation of the equinoxes.
2026          * This date corresponds to 1997-02-27T00:00:00 UTC, hence the 30s offset from TAI.
2027          */
2028         private static final AbsoluteDate MODEL_START =
2029             new AbsoluteDate(1997, 2, 27, 0, 0, 30, TimeScalesFactory.getTAI());
2030 
2031         /** Evaluate the correction.
2032          * @param arguments Delaunay for nutation
2033          * @return correction value (0 before 1997-02-27)
2034          */
2035         public static double value(final DelaunayArguments arguments) {
2036             if (arguments.getDate().compareTo(MODEL_START) >= 0) {
2037 
2038                 // IAU 1994 resolution C7 added two terms to the equation of equinoxes
2039                 // taking effect since 1997-02-27 for continuity
2040 
2041                 // Mean longitude of the ascending node of the Moon
2042                 final double om = arguments.getOmega();
2043 
2044                 // add the two correction terms
2045                 return EQE1 * FastMath.sin(om) + EQE2 * FastMath.sin(om + om);
2046 
2047             } else {
2048                 return 0.0;
2049             }
2050         }
2051 
2052     };
2053 
2054     /** Stellar angle model.
2055      * <p>
2056      * The stellar angle computed here has been defined in the paper "A non-rotating origin on the
2057      * instantaneous equator: Definition, properties and use", N. Capitaine, Guinot B., and Souchay J.,
2058      * Celestial Mechanics, Volume 39, Issue 3, pp 283-307. It has been proposed as a conventional
2059      * conventional relationship between UT1 and Earth rotation in the paper "Definition of the Celestial
2060      * Ephemeris origin and of UT1 in the International Celestial Reference Frame", Capitaine, N.,
2061      * Guinot, B., and McCarthy, D. D., 2000, “,” Astronomy and Astrophysics, 355(1), pp. 398–405.
2062      * </p>
2063      * <p>
2064      * It is presented simply as stellar angle in IERS conventions 1996 but as since been adopted as
2065      * the conventional relationship defining UT1 from ICRF and is called Earth Rotation Angle in
2066      * IERS conventions 2003 and 2010.
2067      * </p>
2068      */
2069     private static class StellarAngleCapitaine implements TimeFunction<DerivativeStructure> {
2070 
2071         /** Reference date of Capitaine's Earth Rotation Angle model. */
2072         private static final AbsoluteDate REFERENCE_DATE = new AbsoluteDate(DateComponents.J2000_EPOCH,
2073                                                                             TimeComponents.H12,
2074                                                                             TimeScalesFactory.getTAI());
2075 
2076         /** Constant term of Capitaine's Earth Rotation Angle model. */
2077         private static final double ERA_0   = MathUtils.TWO_PI * 0.7790572732640;
2078 
2079         /** Rate term of Capitaine's Earth Rotation Angle model.
2080          * (radians per day, main part) */
2081         private static final double ERA_1A  = MathUtils.TWO_PI / Constants.JULIAN_DAY;
2082 
2083         /** Rate term of Capitaine's Earth Rotation Angle model.
2084          * (radians per day, fractional part) */
2085         private static final double ERA_1B  = ERA_1A * 0.00273781191135448;
2086 
2087         /** Total rate term of Capitaine's Earth Rotation Angle model. */
2088         private static final double ERA_1AB = ERA_1A + ERA_1B;
2089 
2090         /** UT1 time scale. */
2091         private final TimeScale ut1;
2092 
2093         /** Simple constructor.
2094          * @param ut1 UT1 time scale
2095          */
2096         public StellarAngleCapitaine(final TimeScale ut1) {
2097             this.ut1 = ut1;
2098         }
2099 
2100         /** {@inheritDoc} */
2101         @Override
2102         public DerivativeStructure value(final AbsoluteDate date) {
2103 
2104             // split the date offset as a full number of days plus a smaller part
2105             final int secondsInDay = 86400;
2106             final double dt  = date.durationFrom(REFERENCE_DATE);
2107             final long days  = ((long) dt) / secondsInDay;
2108             final double dtA = secondsInDay * days;
2109             final double dtB = (dt - dtA) + ut1.offsetFromTAI(date);
2110 
2111             return new DerivativeStructure(1, 1,
2112                                            ERA_0 + ERA_1A * dtB + ERA_1B * (dtA + dtB),
2113                                            ERA_1AB);
2114 
2115         }
2116 
2117     }
2118 
2119     /** Mean pole. */
2120     private static class MeanPole implements TimeStamped, Serializable {
2121 
2122         /** Serializable UID. */
2123         private static final long serialVersionUID = 20131028l;
2124 
2125         /** Date. */
2126         private final AbsoluteDate date;
2127 
2128         /** X coordinate. */
2129         private double x;
2130 
2131         /** Y coordinate. */
2132         private double y;
2133 
2134         /** Simple constructor.
2135          * @param date date
2136          * @param x x coordinate
2137          * @param y y coordinate
2138          */
2139         public MeanPole(final AbsoluteDate date, final double x, final double y) {
2140             this.date = date;
2141             this.x    = x;
2142             this.y    = y;
2143         }
2144 
2145         /** {@inheritDoc} */
2146         @Override
2147         public AbsoluteDate getDate() {
2148             return date;
2149         }
2150 
2151         /** Get x coordinate.
2152          * @return x coordinate
2153          */
2154         public double getX() {
2155             return x;
2156         }
2157 
2158         /** Get y coordinate.
2159          * @return y coordinate
2160          */
2161         public double getY() {
2162             return y;
2163         }
2164 
2165     }
2166 }