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