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