IERSConventions.java

  1. /* Copyright 2002-2020 CS Group
  2.  * Licensed to CS Group (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.utils;

  18. import java.io.BufferedReader;
  19. import java.io.IOException;
  20. import java.io.InputStream;
  21. import java.io.InputStreamReader;
  22. import java.io.Serializable;
  23. import java.nio.charset.StandardCharsets;
  24. import java.util.List;

  25. import org.hipparchus.RealFieldElement;
  26. import org.hipparchus.analysis.interpolation.FieldHermiteInterpolator;
  27. import org.hipparchus.analysis.interpolation.HermiteInterpolator;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.MathArrays;
  30. import org.hipparchus.util.MathUtils;
  31. import org.orekit.annotation.DefaultDataContext;
  32. import org.orekit.data.BodiesElements;
  33. import org.orekit.data.DataContext;
  34. import org.orekit.data.DelaunayArguments;
  35. import org.orekit.data.FieldBodiesElements;
  36. import org.orekit.data.FieldDelaunayArguments;
  37. import org.orekit.data.FundamentalNutationArguments;
  38. import org.orekit.data.PoissonSeries;
  39. import org.orekit.data.PoissonSeries.CompiledSeries;
  40. import org.orekit.data.PoissonSeriesParser;
  41. import org.orekit.data.PolynomialNutation;
  42. import org.orekit.data.PolynomialParser;
  43. import org.orekit.data.PolynomialParser.Unit;
  44. import org.orekit.data.SimpleTimeStampedTableParser;
  45. import org.orekit.errors.OrekitException;
  46. import org.orekit.errors.OrekitInternalError;
  47. import org.orekit.errors.OrekitMessages;
  48. import org.orekit.errors.TimeStampedCacheException;
  49. import org.orekit.frames.EOPHistory;
  50. import org.orekit.frames.FieldPoleCorrection;
  51. import org.orekit.frames.PoleCorrection;
  52. import org.orekit.time.AbsoluteDate;
  53. import org.orekit.time.DateComponents;
  54. import org.orekit.time.FieldAbsoluteDate;
  55. import org.orekit.time.TimeComponents;
  56. import org.orekit.time.TimeScalarFunction;
  57. import org.orekit.time.TimeScale;
  58. import org.orekit.time.TimeScales;
  59. import org.orekit.time.TimeStamped;
  60. import org.orekit.time.TimeVectorFunction;


  61. /** Supported IERS conventions.
  62.  * @since 6.0
  63.  * @author Luc Maisonobe
  64.  */
  65. public enum IERSConventions {

  66.     /** Constant for IERS 1996 conventions. */
  67.     IERS_1996 {

  68.         /** Nutation arguments resources. */
  69.         private static final String NUTATION_ARGUMENTS = IERS_BASE + "1996/nutation-arguments.txt";

  70.         /** X series resources. */
  71.         private static final String X_Y_SERIES         = IERS_BASE + "1996/tab5.4.txt";

  72.         /** Psi series resources. */
  73.         private static final String PSI_EPSILON_SERIES = IERS_BASE + "1996/tab5.1.txt";

  74.         /** Tidal correction for xp, yp series resources. */
  75.         private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "1996/tab8.4.txt";

  76.         /** Tidal correction for UT1 resources. */
  77.         private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "1996/tab8.3.txt";

  78.         /** Love numbers resources. */
  79.         private static final String LOVE_NUMBERS = IERS_BASE + "1996/tab6.1.txt";

  80.         /** Frequency dependence model for k₂₀. */
  81.         private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2b.txt";

  82.         /** Frequency dependence model for k₂₁. */
  83.         private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2a.txt";

  84.         /** Frequency dependence model for k₂₂. */
  85.         private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2c.txt";

  86.         /** Tidal displacement frequency correction for diurnal tides. */
  87.         private static final String TIDAL_DISPLACEMENT_CORRECTION_DIURNAL = IERS_BASE + "1996/tab7.3a.txt";

  88.         /** Tidal displacement frequency correction for zonal tides. */
  89.         private static final String TIDAL_DISPLACEMENT_CORRECTION_ZONAL = IERS_BASE + "1996/tab7.3b.txt";

  90.         /** {@inheritDoc} */
  91.         @Override
  92.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
  93.                                                                  final TimeScales timeScales) {
  94.             try (InputStream in = getStream(NUTATION_ARGUMENTS)) {
  95.                 return new FundamentalNutationArguments(this, timeScale,
  96.                         in, NUTATION_ARGUMENTS, timeScales);
  97.             } catch (IOException e) {
  98.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  99.             }
  100.         }

  101.         /** {@inheritDoc} */
  102.         @Override
  103.         public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {

  104.             // value from chapter 5, page 22
  105.             final PolynomialNutation epsilonA =
  106.                     new PolynomialNutation(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  107.                                              -46.8150   * Constants.ARC_SECONDS_TO_RADIANS,
  108.                                               -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  109.                                                0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  110.             return new TimeScalarFunction() {

  111.                 /** {@inheritDoc} */
  112.                 @Override
  113.                 public double value(final AbsoluteDate date) {
  114.                     return epsilonA.value(evaluateTC(date, timeScales));
  115.                 }

  116.                 /** {@inheritDoc} */
  117.                 @Override
  118.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  119.                     return epsilonA.value(evaluateTC(date, timeScales));
  120.                 }

  121.             };

  122.         }

  123.         /** {@inheritDoc} */
  124.         @Override
  125.         public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {

  126.             // set up nutation arguments
  127.             final FundamentalNutationArguments arguments =
  128.                     getNutationArguments(timeScales);

  129.             // X = 2004.3109″t - 0.42665″t² - 0.198656″t³ + 0.0000140″t⁴
  130.             //     + 0.00006″t² cos Ω + sin ε0 { Σ [(Ai + Ai' t) sin(ARGUMENT) + Ai'' t cos(ARGUMENT)]}
  131.             //     + 0.00204″t² sin Ω + 0.00016″t² sin 2(F - D + Ω),
  132.             final PolynomialNutation xPolynomial =
  133.                     new PolynomialNutation(0,
  134.                                            2004.3109 * Constants.ARC_SECONDS_TO_RADIANS,
  135.                                            -0.42665  * Constants.ARC_SECONDS_TO_RADIANS,
  136.                                            -0.198656 * Constants.ARC_SECONDS_TO_RADIANS,
  137.                                            0.0000140 * Constants.ARC_SECONDS_TO_RADIANS);

  138.             final double fXCosOm    = 0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
  139.             final double fXSinOm    = 0.00204 * Constants.ARC_SECONDS_TO_RADIANS;
  140.             final double fXSin2FDOm = 0.00016 * Constants.ARC_SECONDS_TO_RADIANS;
  141.             final double sinEps0   = FastMath.sin(getMeanObliquityFunction(timeScales)
  142.                     .value(getNutationReferenceEpoch(timeScales)));

  143.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  144.             final PoissonSeriesParser baseParser =
  145.                     new PoissonSeriesParser(12).withFirstDelaunay(1);

  146.             final PoissonSeriesParser xParser =
  147.                     baseParser.
  148.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  149.                     withSinCos(1, 8, deciMilliAS,  9, deciMilliAS);
  150.             final PoissonSeries xSum;
  151.             try (InputStream in = getStream(X_Y_SERIES)) {
  152.                 xSum = xParser.parse(in, X_Y_SERIES);
  153.             } catch (IOException e) {
  154.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  155.             }

  156.             // Y = -0.00013″ - 22.40992″t² + 0.001836″t³ + 0.0011130″t⁴
  157.             //     + Σ [(Bi + Bi' t) cos(ARGUMENT) + Bi'' t sin(ARGUMENT)]
  158.             //    - 0.00231″t² cos Ω − 0.00014″t² cos 2(F - D + Ω)
  159.             final PolynomialNutation yPolynomial =
  160.                     new PolynomialNutation(-0.00013  * Constants.ARC_SECONDS_TO_RADIANS,
  161.                                            0.0,
  162.                                            -22.40992 * Constants.ARC_SECONDS_TO_RADIANS,
  163.                                            0.001836  * Constants.ARC_SECONDS_TO_RADIANS,
  164.                                            0.0011130 * Constants.ARC_SECONDS_TO_RADIANS);

  165.             final double fYCosOm    = -0.00231 * Constants.ARC_SECONDS_TO_RADIANS;
  166.             final double fYCos2FDOm = -0.00014 * Constants.ARC_SECONDS_TO_RADIANS;

  167.             final PoissonSeriesParser yParser =
  168.                     baseParser.
  169.                     withSinCos(0, -1, deciMilliAS, 10, deciMilliAS).
  170.                     withSinCos(1, 12, deciMilliAS, 11, deciMilliAS);
  171.             final PoissonSeries ySum;
  172.             try (InputStream in = getStream(X_Y_SERIES)) {
  173.                 ySum = yParser.parse(in, X_Y_SERIES);
  174.             } catch (IOException e) {
  175.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  176.             }

  177.             final PoissonSeries.CompiledSeries xySum =
  178.                     PoissonSeries.compile(xSum, ySum);

  179.             // s = -XY/2 + 0.00385″t - 0.07259″t³ - 0.00264″ sin Ω - 0.00006″ sin 2Ω
  180.             //     + 0.00074″t² sin Ω + 0.00006″t² sin 2(F - D + Ω)
  181.             final double fST          =  0.00385 * Constants.ARC_SECONDS_TO_RADIANS;
  182.             final double fST3         = -0.07259 * Constants.ARC_SECONDS_TO_RADIANS;
  183.             final double fSSinOm      = -0.00264 * Constants.ARC_SECONDS_TO_RADIANS;
  184.             final double fSSin2Om     = -0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
  185.             final double fST2SinOm    =  0.00074 * Constants.ARC_SECONDS_TO_RADIANS;
  186.             final double fST2Sin2FDOm =  0.00006 * Constants.ARC_SECONDS_TO_RADIANS;

  187.             return new TimeVectorFunction() {

  188.                 /** {@inheritDoc} */
  189.                 @Override
  190.                 public double[] value(final AbsoluteDate date) {

  191.                     final BodiesElements elements = arguments.evaluateAll(date);
  192.                     final double[] xy             = xySum.value(elements);

  193.                     final double omega     = elements.getOmega();
  194.                     final double f         = elements.getF();
  195.                     final double d         = elements.getD();
  196.                     final double t         = elements.getTC();

  197.                     final double cosOmega  = FastMath.cos(omega);
  198.                     final double sinOmega  = FastMath.sin(omega);
  199.                     final double sin2Omega = FastMath.sin(2 * omega);
  200.                     final double cos2FDOm  = FastMath.cos(2 * (f - d + omega));
  201.                     final double sin2FDOm  = FastMath.sin(2 * (f - d + omega));

  202.                     final double x = xPolynomial.value(t) + sinEps0 * xy[0] +
  203.                             t * t * (fXCosOm * cosOmega + fXSinOm * sinOmega + fXSin2FDOm * cos2FDOm);
  204.                     final double y = yPolynomial.value(t) + xy[1] +
  205.                             t * t * (fYCosOm * cosOmega + fYCos2FDOm * cos2FDOm);
  206.                     final double sPxy2 = fSSinOm * sinOmega + fSSin2Om * sin2Omega +
  207.                             t * (fST + t * (fST2SinOm * sinOmega + fST2Sin2FDOm * sin2FDOm + t * fST3));

  208.                     return new double[] {
  209.                         x, y, sPxy2
  210.                     };

  211.                 }

  212.                 /** {@inheritDoc} */
  213.                 @Override
  214.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

  215.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  216.                     final T[] xy             = xySum.value(elements);

  217.                     final T omega     = elements.getOmega();
  218.                     final T f         = elements.getF();
  219.                     final T d         = elements.getD();
  220.                     final T t         = elements.getTC();
  221.                     final T t2        = t.multiply(t);

  222.                     final T cosOmega  = omega.cos();
  223.                     final T sinOmega  = omega.sin();
  224.                     final T sin2Omega = omega.multiply(2).sin();
  225.                     final T fMDPO2 = f.subtract(d).add(omega).multiply(2);
  226.                     final T cos2FDOm  = fMDPO2.cos();
  227.                     final T sin2FDOm  = fMDPO2.sin();

  228.                     final T x = xPolynomial.value(t).
  229.                                 add(xy[0].multiply(sinEps0)).
  230.                                 add(t2.multiply(cosOmega.multiply(fXCosOm).add(sinOmega.multiply(fXSinOm)).add(cos2FDOm.multiply(fXSin2FDOm))));
  231.                     final T y = yPolynomial.value(t).
  232.                                 add(xy[1]).
  233.                                 add(t2.multiply(cosOmega.multiply(fYCosOm).add(cos2FDOm.multiply(fYCos2FDOm))));
  234.                     final T sPxy2 = sinOmega.multiply(fSSinOm).
  235.                                     add(sin2Omega.multiply(fSSin2Om)).
  236.                                     add(t.multiply(fST3).add(sinOmega.multiply(fST2SinOm)).add(sin2FDOm.multiply(fST2Sin2FDOm)).multiply(t).add(fST).multiply(t));

  237.                     final T[] a = MathArrays.buildArray(date.getField(), 3);
  238.                     a[0] = x;
  239.                     a[1] = y;
  240.                     a[2] = sPxy2;
  241.                     return a;

  242.                 }

  243.             };

  244.         }

  245.         /** {@inheritDoc} */
  246.         @Override
  247.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

  248.             // set up the conventional polynomials
  249.             // the following values are from Lieske et al. paper:
  250.             // Expressions for the precession quantities based upon the IAU(1976) system of astronomical constants
  251.             // http://articles.adsabs.harvard.edu/full/1977A%26A....58....1L
  252.             // also available as equation 30 in IERS 2003 conventions
  253.             final PolynomialNutation psiA =
  254.                     new PolynomialNutation(   0.0,
  255.                                            5038.7784   * Constants.ARC_SECONDS_TO_RADIANS,
  256.                                              -1.07259  * Constants.ARC_SECONDS_TO_RADIANS,
  257.                                              -0.001147 * Constants.ARC_SECONDS_TO_RADIANS);
  258.             final PolynomialNutation omegaA =
  259.                     new PolynomialNutation(getMeanObliquityFunction(timeScales)
  260.                             .value(getNutationReferenceEpoch(timeScales)),
  261.                                             0.0,
  262.                                             0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
  263.                                            -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
  264.             final PolynomialNutation chiA =
  265.                     new PolynomialNutation( 0.0,
  266.                                            10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
  267.                                            -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
  268.                                            -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);

  269.             return new PrecessionFunction(psiA, omegaA, chiA, timeScales);

  270.         }

  271.         /** {@inheritDoc} */
  272.         @Override
  273.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  274.             // set up nutation arguments
  275.             final FundamentalNutationArguments arguments =
  276.                     getNutationArguments(timeScales);

  277.             // set up Poisson series
  278.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  279.             final PoissonSeriesParser baseParser =
  280.                     new PoissonSeriesParser(10).withFirstDelaunay(1);

  281.             final PoissonSeriesParser psiParser =
  282.                     baseParser.
  283.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  284.                     withSinCos(1, 8, deciMilliAS, -1, deciMilliAS);
  285.             final PoissonSeries psiSeries;
  286.             try (InputStream in = getStream(PSI_EPSILON_SERIES)) {
  287.                 psiSeries = psiParser.parse(in, PSI_EPSILON_SERIES);
  288.             } catch (IOException e) {
  289.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  290.             }

  291.             final PoissonSeriesParser epsilonParser =
  292.                     baseParser.
  293.                     withSinCos(0, -1, deciMilliAS, 9, deciMilliAS).
  294.                     withSinCos(1, -1, deciMilliAS, 10, deciMilliAS);
  295.             final PoissonSeries epsilonSeries;
  296.             try (InputStream in = getStream(PSI_EPSILON_SERIES)) {
  297.                 epsilonSeries = epsilonParser.parse(in, PSI_EPSILON_SERIES);
  298.             } catch (IOException e) {
  299.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  300.             }

  301.             final PoissonSeries.CompiledSeries psiEpsilonSeries =
  302.                     PoissonSeries.compile(psiSeries, epsilonSeries);

  303.             return new TimeVectorFunction() {

  304.                 /** {@inheritDoc} */
  305.                 @Override
  306.                 public double[] value(final AbsoluteDate date) {
  307.                     final BodiesElements elements = arguments.evaluateAll(date);
  308.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  309.                     return new double[] {
  310.                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  311.                     };
  312.                 }

  313.                 /** {@inheritDoc} */
  314.                 @Override
  315.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  316.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  317.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  318.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  319.                     result[0] = psiEpsilon[0];
  320.                     result[1] = psiEpsilon[1];
  321.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  322.                     return result;
  323.                 }

  324.             };

  325.         }

  326.         /** {@inheritDoc} */
  327.         @Override
  328.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  329.                                                   final TimeScales timeScales) {

  330.             // Radians per second of time
  331.             final double radiansPerSecond = MathUtils.TWO_PI / Constants.JULIAN_DAY;

  332.             // constants from IERS 1996 page 21
  333.             // the underlying model is IAU 1982 GMST-UT1
  334.             final AbsoluteDate gmstReference = new AbsoluteDate(
  335.                     DateComponents.J2000_EPOCH, TimeComponents.H12, timeScales.getTAI());
  336.             final double gmst0 = 24110.54841;
  337.             final double gmst1 = 8640184.812866;
  338.             final double gmst2 = 0.093104;
  339.             final double gmst3 = -6.2e-6;

  340.             return new TimeScalarFunction() {

  341.                 /** {@inheritDoc} */
  342.                 @Override
  343.                 public double value(final AbsoluteDate date) {

  344.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  345.                     final double dtai = date.durationFrom(gmstReference);
  346.                     final double tut1 = dtai + ut1.offsetFromTAI(date);
  347.                     final double tt   = tut1 / Constants.JULIAN_CENTURY;

  348.                     // Seconds in the day, adjusted by 12 hours because the
  349.                     // UT1 is supplied as a Julian date beginning at noon.
  350.                     final double sd = FastMath.IEEEremainder(tut1 + Constants.JULIAN_DAY / 2, Constants.JULIAN_DAY);

  351.                     // compute Greenwich mean sidereal time, in radians
  352.                     return ((((((tt * gmst3 + gmst2) * tt) + gmst1) * tt) + gmst0) + sd) * radiansPerSecond;

  353.                 }

  354.                 /** {@inheritDoc} */
  355.                 @Override
  356.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  357.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  358.                     final T dtai = date.durationFrom(gmstReference);
  359.                     final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
  360.                     final T tt   = tut1.divide(Constants.JULIAN_CENTURY);

  361.                     // Seconds in the day, adjusted by 12 hours because the
  362.                     // UT1 is supplied as a Julian date beginning at noon.
  363.                     final T sd = tut1.add(Constants.JULIAN_DAY / 2).remainder(Constants.JULIAN_DAY);

  364.                     // compute Greenwich mean sidereal time, in radians
  365.                     return tt.multiply(gmst3).add(gmst2).multiply(tt).add(gmst1).multiply(tt).add(gmst0).add(sd).multiply(radiansPerSecond);

  366.                 }

  367.             };

  368.         }

  369.         /** {@inheritDoc} */
  370.         @Override
  371.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  372.                                                       final TimeScales timeScales) {

  373.             // Radians per second of time
  374.             final double radiansPerSecond = MathUtils.TWO_PI / Constants.JULIAN_DAY;

  375.             // constants from IERS 1996 page 21
  376.             // the underlying model is IAU 1982 GMST-UT1
  377.             final AbsoluteDate gmstReference = new AbsoluteDate(
  378.                     DateComponents.J2000_EPOCH, TimeComponents.H12, timeScales.getTAI());
  379.             final double gmst1 = 8640184.812866;
  380.             final double gmst2 = 0.093104;
  381.             final double gmst3 = -6.2e-6;

  382.             return new TimeScalarFunction() {

  383.                 /** {@inheritDoc} */
  384.                 @Override
  385.                 public double value(final AbsoluteDate date) {

  386.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  387.                     final double dtai = date.durationFrom(gmstReference);
  388.                     final double tut1 = dtai + ut1.offsetFromTAI(date);
  389.                     final double tt   = tut1 / Constants.JULIAN_CENTURY;

  390.                     // compute Greenwich mean sidereal time rate, in radians per second
  391.                     return ((((tt * 3 * gmst3 + 2 * gmst2) * tt) + gmst1) / Constants.JULIAN_CENTURY + 1) * radiansPerSecond;

  392.                 }

  393.                 /** {@inheritDoc} */
  394.                 @Override
  395.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  396.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  397.                     final T dtai = date.durationFrom(gmstReference);
  398.                     final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
  399.                     final T tt   = tut1.divide(Constants.JULIAN_CENTURY);

  400.                     // compute Greenwich mean sidereal time, in radians
  401.                     return tt.multiply(3 * gmst3).add(2 * gmst2).multiply(tt).add(gmst1).divide(Constants.JULIAN_CENTURY).add(1).multiply(radiansPerSecond);

  402.                 }

  403.             };

  404.         }

  405.         /** {@inheritDoc} */
  406.         @Override
  407.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  408.                                                   final EOPHistory eopHistory,
  409.                                                   final TimeScales timeScales) {

  410.             // obliquity
  411.             final TimeScalarFunction epsilonA = getMeanObliquityFunction(timeScales);

  412.             // GMST function
  413.             final TimeScalarFunction gmst = getGMSTFunction(ut1, timeScales);

  414.             // nutation function
  415.             final TimeVectorFunction nutation = getNutationFunction(timeScales);

  416.             return new TimeScalarFunction() {

  417.                 /** {@inheritDoc} */
  418.                 @Override
  419.                 public double value(final AbsoluteDate date) {

  420.                     // compute equation of equinoxes
  421.                     final double[] angles = nutation.value(date);
  422.                     double deltaPsi = angles[0];
  423.                     if (eopHistory != null) {
  424.                         deltaPsi += eopHistory.getEquinoxNutationCorrection(date)[0];
  425.                     }
  426.                     final double eqe = deltaPsi  * FastMath.cos(epsilonA.value(date)) + angles[2];

  427.                     // add mean sidereal time and equation of equinoxes
  428.                     return gmst.value(date) + eqe;

  429.                 }

  430.                 /** {@inheritDoc} */
  431.                 @Override
  432.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  433.                     // compute equation of equinoxes
  434.                     final T[] angles = nutation.value(date);
  435.                     T deltaPsi = angles[0];
  436.                     if (eopHistory != null) {
  437.                         deltaPsi = deltaPsi.add(eopHistory.getEquinoxNutationCorrection(date)[0]);
  438.                     }
  439.                     final T eqe = deltaPsi.multiply(epsilonA.value(date).cos()).add(angles[2]);

  440.                     // add mean sidereal time and equation of equinoxes
  441.                     return gmst.value(date).add(eqe);

  442.                 }

  443.             };

  444.         }

  445.         /** {@inheritDoc} */
  446.         @Override
  447.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

  448.             // set up nutation arguments
  449.             // BEWARE! Using TT as the time scale here and not UT1 is intentional!
  450.             // as this correction is used to compute UT1 itself, it is not surprising we cannot use UT1 yet,
  451.             // however, using the close UTC as would seem logical make the comparison with interp.f from IERS fail
  452.             // looking in the interp.f code, the same TT scale is used for both Delaunay and gamma argument
  453.             final FundamentalNutationArguments arguments =
  454.                     getNutationArguments(timeScales.getTT(), timeScales);

  455.             // set up Poisson series
  456.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  457.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(17).
  458.                     withOptionalColumn(1).
  459.                     withGamma(7).
  460.                     withFirstDelaunay(2);
  461.             final PoissonSeries xSeries;
  462.             try (InputStream in = getStream(TIDAL_CORRECTION_XP_YP_SERIES)) {
  463.                 xSeries = xyParser.
  464.                         withSinCos(0, 14, milliAS, 15, milliAS).
  465.                         parse(in, TIDAL_CORRECTION_XP_YP_SERIES);
  466.             } catch (IOException e) {
  467.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  468.             }
  469.             final PoissonSeries ySeries;
  470.             try (InputStream in = getStream(TIDAL_CORRECTION_XP_YP_SERIES)) {
  471.                 ySeries = xyParser.
  472.                         withSinCos(0, 16, milliAS, 17, milliAS).
  473.                         parse(in, TIDAL_CORRECTION_XP_YP_SERIES);
  474.             } catch (IOException e) {
  475.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  476.             }

  477.             final double deciMilliS = 1.0e-4;
  478.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(17).
  479.                     withOptionalColumn(1).
  480.                     withGamma(7).
  481.                     withFirstDelaunay(2).
  482.                     withSinCos(0, 16, deciMilliS, 17, deciMilliS);
  483.             final PoissonSeries ut1Series;
  484.             try (InputStream in = getStream(TIDAL_CORRECTION_UT1_SERIES)) {
  485.                 ut1Series = ut1Parser.parse(in, TIDAL_CORRECTION_UT1_SERIES);
  486.             } catch (IOException e) {
  487.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  488.             }

  489.             return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);

  490.         }

  491.         /** {@inheritDoc} */
  492.         public LoveNumbers getLoveNumbers() {
  493.             return loadLoveNumbers(LOVE_NUMBERS);
  494.         }

  495.         /** {@inheritDoc} */
  496.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  497.                                                                      final TimeScales timeScales) {

  498.             // set up nutation arguments
  499.             final FundamentalNutationArguments arguments = getNutationArguments(ut1, timeScales);

  500.             // set up Poisson series
  501.             final PoissonSeriesParser k20Parser =
  502.                     new PoissonSeriesParser(18).
  503.                         withOptionalColumn(1).
  504.                         withDoodson(4, 2).
  505.                         withFirstDelaunay(10);
  506.             final PoissonSeriesParser k21Parser =
  507.                     new PoissonSeriesParser(18).
  508.                         withOptionalColumn(1).
  509.                         withDoodson(4, 3).
  510.                         withFirstDelaunay(10);
  511.             final PoissonSeriesParser k22Parser =
  512.                     new PoissonSeriesParser(16).
  513.                         withOptionalColumn(1).
  514.                         withDoodson(4, 2).
  515.                         withFirstDelaunay(10);

  516.             final double pico = 1.0e-12;
  517.             try (InputStream k20 = getStream(K20_FREQUENCY_DEPENDENCE);
  518.                  InputStream k21In1 = getStream(K21_FREQUENCY_DEPENDENCE);
  519.                  InputStream k21In2 = getStream(K21_FREQUENCY_DEPENDENCE);
  520.                  InputStream k22In1 = getStream(K22_FREQUENCY_DEPENDENCE);
  521.                  InputStream k22In2 = getStream(K22_FREQUENCY_DEPENDENCE)) {
  522.                 final PoissonSeries c20Series =
  523.                         k20Parser.
  524.                         withSinCos(0, 18, -pico, 16, pico).
  525.                         parse(k20, K20_FREQUENCY_DEPENDENCE);
  526.                 final PoissonSeries c21Series =
  527.                         k21Parser.
  528.                         withSinCos(0, 17, pico, 18, pico).
  529.                         parse(k21In1, K21_FREQUENCY_DEPENDENCE);
  530.                 final PoissonSeries s21Series =
  531.                         k21Parser.
  532.                         withSinCos(0, 18, -pico, 17, pico).
  533.                         parse(k21In2, K21_FREQUENCY_DEPENDENCE);
  534.                 final PoissonSeries c22Series =
  535.                         k22Parser.
  536.                         withSinCos(0, -1, pico, 16, pico).
  537.                         parse(k22In1, K22_FREQUENCY_DEPENDENCE);
  538.                 final PoissonSeries s22Series =
  539.                         k22Parser.
  540.                         withSinCos(0, 16, -pico, -1, pico).
  541.                         parse(k22In2, K22_FREQUENCY_DEPENDENCE);

  542.                 return new TideFrequencyDependenceFunction(arguments,
  543.                                                            c20Series,
  544.                                                            c21Series, s21Series,
  545.                                                            c22Series, s22Series);
  546.             } catch (IOException e) {
  547.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  548.             }

  549.         }

  550.         /** {@inheritDoc} */
  551.         @Override
  552.         public double getPermanentTide() {
  553.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  554.         }

  555.         /** {@inheritDoc} */
  556.         @Override
  557.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  558.             // constants from IERS 1996 page 47
  559.             final double globalFactor = -1.348e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  560.             final double coupling     =  0.00112;

  561.             return new TimeVectorFunction() {

  562.                 /** {@inheritDoc} */
  563.                 @Override
  564.                 public double[] value(final AbsoluteDate date) {
  565.                     final PoleCorrection pole = eopHistory.getPoleCorrection(date);
  566.                     return new double[] {
  567.                         globalFactor * (pole.getXp() + coupling * pole.getYp()),
  568.                         globalFactor * (coupling * pole.getXp() - pole.getYp()),
  569.                     };
  570.                 }

  571.                 /** {@inheritDoc} */
  572.                 @Override
  573.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  574.                     final FieldPoleCorrection<T> pole = eopHistory.getPoleCorrection(date);
  575.                     final T[] a = MathArrays.buildArray(date.getField(), 2);
  576.                     a[0] = pole.getXp().add(pole.getYp().multiply(coupling)).multiply(globalFactor);
  577.                     a[1] = pole.getXp().multiply(coupling).subtract(pole.getYp()).multiply(globalFactor);
  578.                     return a;
  579.                 }

  580.             };
  581.         }

  582.         /** {@inheritDoc} */
  583.         @Override
  584.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {

  585.             return new TimeVectorFunction() {

  586.                 /** {@inheritDoc} */
  587.                 @Override
  588.                 public double[] value(final AbsoluteDate date) {
  589.                     // there are no model for ocean pole tide prior to conventions 2010
  590.                     return new double[] {
  591.                         0.0, 0.0
  592.                     };
  593.                 }

  594.                 /** {@inheritDoc} */
  595.                 @Override
  596.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  597.                     // there are no model for ocean pole tide prior to conventions 2010
  598.                     return MathArrays.buildArray(date.getField(), 2);
  599.                 }

  600.             };
  601.         }

  602.         /** {@inheritDoc} */
  603.         @Override
  604.         public double[] getNominalTidalDisplacement() {

  605.             //  // elastic Earth values
  606.             //  return new double[] {
  607.             //      // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  608.             //      0.6026,                                  -0.0006, 0.292, -0.0025, -0.0022,
  609.             //      // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  610.             //      0.0831,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  611.             //      // H₀
  612.             //      -0.31460
  613.             //  };

  614.             // anelastic Earth values
  615.             return new double[] {
  616.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  617.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  618.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  619.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  620.                 // H₀
  621.                 -0.31460
  622.             };

  623.         }

  624.         /** {@inheritDoc} */
  625.         @Override
  626.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
  627.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  628.                                                                   18, 17, -1, 18, -1);
  629.         }

  630.         /** {@inheritDoc} */
  631.         @Override
  632.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  633.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  634.                                                                 20, 17, 19, 18, 20);
  635.         }

  636.     },

  637.     /** Constant for IERS 2003 conventions. */
  638.     IERS_2003 {

  639.         /** Nutation arguments resources. */
  640.         private static final String NUTATION_ARGUMENTS = IERS_BASE + "2003/nutation-arguments.txt";

  641.         /** X series resources. */
  642.         private static final String X_SERIES           = IERS_BASE + "2003/tab5.2a.txt";

  643.         /** Y series resources. */
  644.         private static final String Y_SERIES           = IERS_BASE + "2003/tab5.2b.txt";

  645.         /** S series resources. */
  646.         private static final String S_SERIES           = IERS_BASE + "2003/tab5.2c.txt";

  647.         /** Luni-solar series resources. */
  648.         private static final String LUNI_SOLAR_SERIES  = IERS_BASE + "2003/tab5.3a-first-table.txt";

  649.         /** Planetary series resources. */
  650.         private static final String PLANETARY_SERIES   = IERS_BASE + "2003/tab5.3b.txt";

  651.         /** Greenwhich sidereal time series resources. */
  652.         private static final String GST_SERIES         = IERS_BASE + "2003/tab5.4.txt";

  653.         /** Tidal correction for xp, yp series resources. */
  654.         private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "2003/tab8.2ab.txt";

  655.         /** Tidal correction for UT1 resources. */
  656.         private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "2003/tab8.3ab.txt";

  657.         /** Love numbers resources. */
  658.         private static final String LOVE_NUMBERS = IERS_BASE + "2003/tab6.1.txt";

  659.         /** Frequency dependence model for k₂₀. */
  660.         private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3b.txt";

  661.         /** Frequency dependence model for k₂₁. */
  662.         private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3a.txt";

  663.         /** Frequency dependence model for k₂₂. */
  664.         private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3c.txt";

  665.         /** Annual pole table. */
  666.         private static final String ANNUAL_POLE = IERS_BASE + "2003/annual.pole";

  667.         /** Tidal displacement frequency correction for diurnal tides. */
  668.         private static final String TIDAL_DISPLACEMENT_CORRECTION_DIURNAL = IERS_BASE + "2003/tab7.5a.txt";

  669.         /** Tidal displacement frequency correction for zonal tides. */
  670.         private static final String TIDAL_DISPLACEMENT_CORRECTION_ZONAL = IERS_BASE + "2003/tab7.5b.txt";

  671.         /** {@inheritDoc} */
  672.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
  673.                                                                  final TimeScales timeScales) {
  674.             try (InputStream in = getStream(NUTATION_ARGUMENTS)) {
  675.                 return new FundamentalNutationArguments(this, timeScale,
  676.                         in, NUTATION_ARGUMENTS, timeScales);
  677.             } catch (IOException e) {
  678.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  679.             }
  680.         }

  681.         /** {@inheritDoc} */
  682.         @Override
  683.         public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {

  684.             // epsilon 0 value from chapter 5, page 41, other terms from equation 32 page 45
  685.             final PolynomialNutation epsilonA =
  686.                     new PolynomialNutation(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  687.                                              -46.84024  * Constants.ARC_SECONDS_TO_RADIANS,
  688.                                               -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  689.                                                0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  690.             return new TimeScalarFunction() {

  691.                 /** {@inheritDoc} */
  692.                 @Override
  693.                 public double value(final AbsoluteDate date) {
  694.                     return epsilonA.value(evaluateTC(date, timeScales));
  695.                 }

  696.                 /** {@inheritDoc} */
  697.                 @Override
  698.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  699.                     return epsilonA.value(evaluateTC(date, timeScales));
  700.                 }

  701.             };

  702.         }

  703.         /** {@inheritDoc} */
  704.         @Override
  705.         public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {

  706.             // set up nutation arguments
  707.             final FundamentalNutationArguments arguments =
  708.                     getNutationArguments(timeScales);

  709.             // set up Poisson series
  710.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  711.             final PoissonSeriesParser parser =
  712.                     new PoissonSeriesParser(17).
  713.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  714.                         withFirstDelaunay(4).
  715.                         withFirstPlanetary(9).
  716.                         withSinCos(0, 2, microAS, 3, microAS);

  717.             final PoissonSeries.CompiledSeries xys;
  718.             try (InputStream xIn = getStream(X_SERIES);
  719.                  InputStream yIn = getStream(Y_SERIES);
  720.                  InputStream sIn = getStream(S_SERIES)) {
  721.                 final PoissonSeries xSeries = parser.parse(xIn, X_SERIES);
  722.                 final PoissonSeries ySeries = parser.parse(yIn, Y_SERIES);
  723.                 final PoissonSeries sSeries = parser.parse(sIn, S_SERIES);
  724.                 xys = PoissonSeries.compile(xSeries, ySeries, sSeries);
  725.             } catch (IOException e) {
  726.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  727.             }

  728.             // create a function evaluating the series
  729.             return new TimeVectorFunction() {

  730.                 /** {@inheritDoc} */
  731.                 @Override
  732.                 public double[] value(final AbsoluteDate date) {
  733.                     return xys.value(arguments.evaluateAll(date));
  734.                 }

  735.                 /** {@inheritDoc} */
  736.                 @Override
  737.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  738.                     return xys.value(arguments.evaluateAll(date));
  739.                 }

  740.             };

  741.         }


  742.         /** {@inheritDoc} */
  743.         @Override
  744.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

  745.             // set up the conventional polynomials
  746.             // the following values are from equation 32 in IERS 2003 conventions
  747.             final PolynomialNutation psiA =
  748.                     new PolynomialNutation(    0.0,
  749.                                             5038.47875   * Constants.ARC_SECONDS_TO_RADIANS,
  750.                                               -1.07259   * Constants.ARC_SECONDS_TO_RADIANS,
  751.                                               -0.001147  * Constants.ARC_SECONDS_TO_RADIANS);
  752.             final PolynomialNutation omegaA =
  753.                     new PolynomialNutation(getMeanObliquityFunction(timeScales)
  754.                             .value(getNutationReferenceEpoch(timeScales)),
  755.                                            -0.02524   * Constants.ARC_SECONDS_TO_RADIANS,
  756.                                             0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
  757.                                            -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
  758.             final PolynomialNutation chiA =
  759.                     new PolynomialNutation( 0.0,
  760.                                            10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
  761.                                            -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
  762.                                            -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);

  763.             return new PrecessionFunction(psiA, omegaA, chiA, timeScales);

  764.         }

  765.         /** {@inheritDoc} */
  766.         @Override
  767.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  768.             // set up nutation arguments
  769.             final FundamentalNutationArguments arguments =
  770.                     getNutationArguments(timeScales);

  771.             // set up Poisson series
  772.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  773.             final PoissonSeriesParser luniSolarParser =
  774.                     new PoissonSeriesParser(14).withFirstDelaunay(1);
  775.             final PoissonSeriesParser luniSolarPsiParser =
  776.                     luniSolarParser.
  777.                     withSinCos(0, 7, milliAS, 11, milliAS).
  778.                     withSinCos(1, 8, milliAS, 12, milliAS);
  779.             final PoissonSeries psiLuniSolarSeries;
  780.             try (InputStream in = getStream(LUNI_SOLAR_SERIES)) {
  781.                 psiLuniSolarSeries = luniSolarPsiParser.parse(in, LUNI_SOLAR_SERIES);
  782.             } catch (IOException e) {
  783.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  784.             }
  785.             final PoissonSeriesParser luniSolarEpsilonParser =
  786.                     luniSolarParser.
  787.                     withSinCos(0, 13, milliAS, 9, milliAS).
  788.                     withSinCos(1, 14, milliAS, 10, milliAS);
  789.             final PoissonSeries epsilonLuniSolarSeries;
  790.             try (InputStream in = getStream(LUNI_SOLAR_SERIES)) {
  791.                 epsilonLuniSolarSeries = luniSolarEpsilonParser.parse(in, LUNI_SOLAR_SERIES);
  792.             } catch (IOException e) {
  793.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  794.             }

  795.             final PoissonSeriesParser planetaryParser =
  796.                     new PoissonSeriesParser(21).
  797.                         withFirstDelaunay(2).
  798.                         withFirstPlanetary(7);
  799.             final PoissonSeriesParser planetaryPsiParser =
  800.                     planetaryParser.withSinCos(0, 17, milliAS, 18, milliAS);
  801.             final PoissonSeries psiPlanetarySeries;
  802.             try (InputStream in = getStream(PLANETARY_SERIES)) {
  803.                 psiPlanetarySeries = planetaryPsiParser.parse(in, PLANETARY_SERIES);
  804.             } catch (IOException e) {
  805.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  806.             }
  807.             final PoissonSeriesParser planetaryEpsilonParser =
  808.                     planetaryParser.withSinCos(0, 19, milliAS, 20, milliAS);
  809.             final PoissonSeries epsilonPlanetarySeries;
  810.             try (InputStream in = getStream(PLANETARY_SERIES)) {
  811.                 epsilonPlanetarySeries = planetaryEpsilonParser.parse(in, PLANETARY_SERIES);
  812.             } catch (IOException e) {
  813.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  814.             }

  815.             final PoissonSeries.CompiledSeries luniSolarSeries =
  816.                     PoissonSeries.compile(psiLuniSolarSeries, epsilonLuniSolarSeries);
  817.             final PoissonSeries.CompiledSeries planetarySeries =
  818.                     PoissonSeries.compile(psiPlanetarySeries, epsilonPlanetarySeries);

  819.             return new TimeVectorFunction() {

  820.                 /** {@inheritDoc} */
  821.                 @Override
  822.                 public double[] value(final AbsoluteDate date) {
  823.                     final BodiesElements elements = arguments.evaluateAll(date);
  824.                     final double[] luniSolar = luniSolarSeries.value(elements);
  825.                     final double[] planetary = planetarySeries.value(elements);
  826.                     return new double[] {
  827.                         luniSolar[0] + planetary[0], luniSolar[1] + planetary[1],
  828.                         IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  829.                     };
  830.                 }

  831.                 /** {@inheritDoc} */
  832.                 @Override
  833.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  834.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  835.                     final T[] luniSolar = luniSolarSeries.value(elements);
  836.                     final T[] planetary = planetarySeries.value(elements);
  837.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  838.                     result[0] = luniSolar[0].add(planetary[0]);
  839.                     result[1] = luniSolar[1].add(planetary[1]);
  840.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  841.                     return result;
  842.                 }

  843.             };

  844.         }

  845.         /** {@inheritDoc} */
  846.         @Override
  847.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  848.                                                   final TimeScales timeScales) {

  849.             // Earth Rotation Angle
  850.             final StellarAngleCapitaine era =
  851.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  852.             // Polynomial part of the apparent sidereal time series
  853.             // which is the opposite of Equation of Origins (EO)
  854.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  855.             final PoissonSeriesParser parser =
  856.                     new PoissonSeriesParser(17).
  857.                         withFirstDelaunay(4).
  858.                         withFirstPlanetary(9).
  859.                         withSinCos(0, 2, microAS, 3, microAS).
  860.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  861.             final PolynomialNutation minusEO;
  862.             try (InputStream in = getStream(GST_SERIES)) {
  863.                 minusEO = parser.parse(in, GST_SERIES).getPolynomial();
  864.             } catch (IOException e) {
  865.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  866.             }

  867.             // create a function evaluating the series
  868.             return new TimeScalarFunction() {

  869.                 /** {@inheritDoc} */
  870.                 @Override
  871.                 public double value(final AbsoluteDate date) {
  872.                     return era.value(date) + minusEO.value(evaluateTC(date, timeScales));
  873.                 }

  874.                 /** {@inheritDoc} */
  875.                 @Override
  876.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  877.                     return era.value(date).add(minusEO.value(evaluateTC(date, timeScales)));
  878.                 }

  879.             };

  880.         }

  881.         /** {@inheritDoc} */
  882.         @Override
  883.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  884.                                                       final TimeScales timeScales) {

  885.             // Earth Rotation Angle
  886.             final StellarAngleCapitaine era =
  887.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  888.             // Polynomial part of the apparent sidereal time series
  889.             // which is the opposite of Equation of Origins (EO)
  890.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  891.             final PoissonSeriesParser parser =
  892.                     new PoissonSeriesParser(17).
  893.                         withFirstDelaunay(4).
  894.                         withFirstPlanetary(9).
  895.                         withSinCos(0, 2, microAS, 3, microAS).
  896.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  897.             final PolynomialNutation minusEO;
  898.             try (InputStream in = getStream(GST_SERIES)) {
  899.                 minusEO = parser.parse(in, GST_SERIES).getPolynomial();
  900.             } catch (IOException e) {
  901.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  902.             }

  903.             // create a function evaluating the series
  904.             return new TimeScalarFunction() {

  905.                 /** {@inheritDoc} */
  906.                 @Override
  907.                 public double value(final AbsoluteDate date) {
  908.                     return era.getRate() + minusEO.derivative(evaluateTC(date, timeScales));
  909.                 }

  910.                 /** {@inheritDoc} */
  911.                 @Override
  912.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  913.                     return minusEO.derivative(evaluateTC(date, timeScales)).add(era.getRate());
  914.                 }

  915.             };

  916.         }

  917.         /** {@inheritDoc} */
  918.         @Override
  919.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  920.                                                   final EOPHistory eopHistory,
  921.                                                   final TimeScales timeScales) {

  922.             // set up nutation arguments
  923.             final FundamentalNutationArguments arguments =
  924.                     getNutationArguments(timeScales);

  925.             // mean obliquity function
  926.             final TimeScalarFunction epsilon = getMeanObliquityFunction(timeScales);

  927.             // set up Poisson series
  928.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  929.             final PoissonSeriesParser luniSolarPsiParser =
  930.                     new PoissonSeriesParser(14).
  931.                         withFirstDelaunay(1).
  932.                         withSinCos(0, 7, milliAS, 11, milliAS).
  933.                         withSinCos(1, 8, milliAS, 12, milliAS);
  934.             final PoissonSeries psiLuniSolarSeries;
  935.             try (InputStream in = getStream(LUNI_SOLAR_SERIES)) {
  936.                 psiLuniSolarSeries = luniSolarPsiParser.parse(in, LUNI_SOLAR_SERIES);
  937.             } catch (IOException e) {
  938.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  939.             }

  940.             final PoissonSeriesParser planetaryPsiParser =
  941.                     new PoissonSeriesParser(21).
  942.                         withFirstDelaunay(2).
  943.                         withFirstPlanetary(7).
  944.                         withSinCos(0, 17, milliAS, 18, milliAS);
  945.             final PoissonSeries psiPlanetarySeries;
  946.             try (InputStream in = getStream(PLANETARY_SERIES)) {
  947.                 psiPlanetarySeries = planetaryPsiParser.parse(in, PLANETARY_SERIES);
  948.             } catch (IOException e) {
  949.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  950.             }


  951.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  952.             final PoissonSeriesParser gstParser =
  953.                     new PoissonSeriesParser(17).
  954.                         withFirstDelaunay(4).
  955.                         withFirstPlanetary(9).
  956.                         withSinCos(0, 2, microAS, 3, microAS).
  957.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  958.             final PoissonSeries gstSeries;
  959.             try (InputStream in = getStream(GST_SERIES)) {
  960.                 gstSeries = gstParser.parse(in, GST_SERIES);
  961.             } catch (IOException e) {
  962.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  963.             }
  964.             final PoissonSeries.CompiledSeries psiGstSeries =
  965.                     PoissonSeries.compile(psiLuniSolarSeries, psiPlanetarySeries, gstSeries);

  966.             // ERA function
  967.             final TimeScalarFunction era =
  968.                     getEarthOrientationAngleFunction(ut1, timeScales.getTAI());

  969.             return new TimeScalarFunction() {

  970.                 /** {@inheritDoc} */
  971.                 @Override
  972.                 public double value(final AbsoluteDate date) {

  973.                     // evaluate equation of origins
  974.                     final BodiesElements elements = arguments.evaluateAll(date);
  975.                     final double[] angles = psiGstSeries.value(elements);
  976.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  977.                     final double deltaPsi = angles[0] + angles[1] + ddPsi;
  978.                     final double epsilonA = epsilon.value(date);

  979.                     // subtract equation of origin from EA
  980.                     // (hence add the series above which have the sign included)
  981.                     return era.value(date) + deltaPsi * FastMath.cos(epsilonA) + angles[2];

  982.                 }

  983.                 /** {@inheritDoc} */
  984.                 @Override
  985.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  986.                     // evaluate equation of origins
  987.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  988.                     final T[] angles = psiGstSeries.value(elements);
  989.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  990.                     final T deltaPsi = angles[0].add(angles[1]).add(ddPsi);
  991.                     final T epsilonA = epsilon.value(date);

  992.                     // subtract equation of origin from EA
  993.                     // (hence add the series above which have the sign included)
  994.                     return era.value(date).add(deltaPsi.multiply(epsilonA.cos())).add(angles[2]);

  995.                 }

  996.             };

  997.         }

  998.         /** {@inheritDoc} */
  999.         @Override
  1000.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

  1001.             // set up nutation arguments
  1002.             // BEWARE! Using TT as the time scale here and not UT1 is intentional!
  1003.             // as this correction is used to compute UT1 itself, it is not surprising we cannot use UT1 yet,
  1004.             // however, using the close UTC as would seem logical make the comparison with interp.f from IERS fail
  1005.             // looking in the interp.f code, the same TT scale is used for both Delaunay and gamma argument
  1006.             final FundamentalNutationArguments arguments =
  1007.                     getNutationArguments(timeScales.getTT(), timeScales);

  1008.             // set up Poisson series
  1009.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1010.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
  1011.                     withOptionalColumn(1).
  1012.                     withGamma(2).
  1013.                     withFirstDelaunay(3);
  1014.             try (InputStream tidalIn1 = getStream(TIDAL_CORRECTION_XP_YP_SERIES);
  1015.                  InputStream tidalIn2 = getStream(TIDAL_CORRECTION_XP_YP_SERIES);
  1016.                  InputStream tidalUt1 = getStream(TIDAL_CORRECTION_UT1_SERIES)) {
  1017.                 final PoissonSeries xSeries =
  1018.                         xyParser.
  1019.                         withSinCos(0, 10, microAS, 11, microAS).
  1020.                         parse(tidalIn1, TIDAL_CORRECTION_XP_YP_SERIES);
  1021.                 final PoissonSeries ySeries =
  1022.                         xyParser.
  1023.                         withSinCos(0, 12, microAS, 13, microAS).
  1024.                         parse(tidalIn2, TIDAL_CORRECTION_XP_YP_SERIES);

  1025.                 final double microS = 1.0e-6;
  1026.                 final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  1027.                         withOptionalColumn(1).
  1028.                         withGamma(2).
  1029.                         withFirstDelaunay(3).
  1030.                         withSinCos(0, 10, microS, 11, microS);
  1031.                 final PoissonSeries ut1Series =
  1032.                         ut1Parser.parse(tidalUt1, TIDAL_CORRECTION_UT1_SERIES);

  1033.                 return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
  1034.             } catch (IOException e) {
  1035.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1036.             }

  1037.         }

  1038.         /** {@inheritDoc} */
  1039.         public LoveNumbers getLoveNumbers() {
  1040.             return loadLoveNumbers(LOVE_NUMBERS);
  1041.         }

  1042.         /** {@inheritDoc} */
  1043.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  1044.                                                                      final TimeScales timeScales) {

  1045.             // set up nutation arguments
  1046.             final FundamentalNutationArguments arguments = getNutationArguments(ut1, timeScales);

  1047.             // set up Poisson series
  1048.             final PoissonSeriesParser k20Parser =
  1049.                     new PoissonSeriesParser(18).
  1050.                         withOptionalColumn(1).
  1051.                         withDoodson(4, 2).
  1052.                         withFirstDelaunay(10);
  1053.             final PoissonSeriesParser k21Parser =
  1054.                     new PoissonSeriesParser(18).
  1055.                         withOptionalColumn(1).
  1056.                         withDoodson(4, 3).
  1057.                         withFirstDelaunay(10);
  1058.             final PoissonSeriesParser k22Parser =
  1059.                     new PoissonSeriesParser(16).
  1060.                         withOptionalColumn(1).
  1061.                         withDoodson(4, 2).
  1062.                         withFirstDelaunay(10);

  1063.             final double pico = 1.0e-12;
  1064.             try (InputStream k20In = getStream(K20_FREQUENCY_DEPENDENCE);
  1065.                  InputStream k21In1 = getStream(K21_FREQUENCY_DEPENDENCE);
  1066.                  InputStream k21In2 = getStream(K21_FREQUENCY_DEPENDENCE);
  1067.                  InputStream k22In1 = getStream(K22_FREQUENCY_DEPENDENCE);
  1068.                  InputStream k22In2 = getStream(K22_FREQUENCY_DEPENDENCE)) {
  1069.                 final PoissonSeries c20Series =
  1070.                         k20Parser.
  1071.                         withSinCos(0, 18, -pico, 16, pico).
  1072.                         parse(k20In, K20_FREQUENCY_DEPENDENCE);
  1073.                 final PoissonSeries c21Series =
  1074.                         k21Parser.
  1075.                         withSinCos(0, 17, pico, 18, pico).
  1076.                         parse(k21In1, K21_FREQUENCY_DEPENDENCE);
  1077.                 final PoissonSeries s21Series =
  1078.                         k21Parser.
  1079.                         withSinCos(0, 18, -pico, 17, pico).
  1080.                         parse(k21In2, K21_FREQUENCY_DEPENDENCE);
  1081.                 final PoissonSeries c22Series =
  1082.                         k22Parser.
  1083.                         withSinCos(0, -1, pico, 16, pico).
  1084.                         parse(k22In1, K22_FREQUENCY_DEPENDENCE);
  1085.                 final PoissonSeries s22Series =
  1086.                         k22Parser.
  1087.                         withSinCos(0, 16, -pico, -1, pico).
  1088.                         parse(k22In2, K22_FREQUENCY_DEPENDENCE);

  1089.                 return new TideFrequencyDependenceFunction(arguments,
  1090.                                                            c20Series,
  1091.                                                            c21Series, s21Series,
  1092.                                                            c22Series, s22Series);
  1093.             } catch (IOException e) {
  1094.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1095.             }

  1096.         }

  1097.         /** {@inheritDoc} */
  1098.         @Override
  1099.         public double getPermanentTide() {
  1100.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1101.         }

  1102.         /** {@inheritDoc} */
  1103.         @Override
  1104.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  1105.             // annual pole from ftp://tai.bipm.org/iers/conv2003/chapter7/annual.pole
  1106.             final TimeScale utc = eopHistory.getTimeScales().getUTC();
  1107.             final SimpleTimeStampedTableParser.RowConverter<MeanPole> converter =
  1108.                 new SimpleTimeStampedTableParser.RowConverter<MeanPole>() {
  1109.                     /** {@inheritDoc} */
  1110.                     @Override
  1111.                     public MeanPole convert(final double[] rawFields) {
  1112.                         return new MeanPole(new AbsoluteDate((int) rawFields[0], 1, 1, utc),
  1113.                                             rawFields[1] * Constants.ARC_SECONDS_TO_RADIANS,
  1114.                                             rawFields[2] * Constants.ARC_SECONDS_TO_RADIANS);
  1115.                     }
  1116.                 };
  1117.             final SimpleTimeStampedTableParser<MeanPole> parser =
  1118.                     new SimpleTimeStampedTableParser<MeanPole>(3, converter);
  1119.             final List<MeanPole> annualPoleList;
  1120.             try (InputStream in = getStream(ANNUAL_POLE)) {
  1121.                 annualPoleList = parser.parse(in, ANNUAL_POLE);
  1122.             } catch (IOException e) {
  1123.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1124.             }
  1125.             final AbsoluteDate firstAnnualPoleDate = annualPoleList.get(0).getDate();
  1126.             final AbsoluteDate lastAnnualPoleDate  = annualPoleList.get(annualPoleList.size() - 1).getDate();
  1127.             final ImmutableTimeStampedCache<MeanPole> annualCache =
  1128.                     new ImmutableTimeStampedCache<MeanPole>(2, annualPoleList);

  1129.             // polynomial extension from IERS 2003, section 7.1.4, equations 23a and 23b
  1130.             final double xp0    = 0.054   * Constants.ARC_SECONDS_TO_RADIANS;
  1131.             final double xp0Dot = 0.00083 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
  1132.             final double yp0    = 0.357   * Constants.ARC_SECONDS_TO_RADIANS;
  1133.             final double yp0Dot = 0.00395 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;

  1134.             // constants from IERS 2003, section 6.2
  1135.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1136.             final double ratio        =  0.00115;

  1137.             return new TimeVectorFunction() {

  1138.                 /** {@inheritDoc} */
  1139.                 @Override
  1140.                 public double[] value(final AbsoluteDate date) {

  1141.                     // we can't compute anything before the range covered by the annual pole file
  1142.                     if (date.compareTo(firstAnnualPoleDate) <= 0) {
  1143.                         return new double[] {
  1144.                             0.0, 0.0
  1145.                         };
  1146.                     }

  1147.                     // evaluate mean pole
  1148.                     double meanPoleX = 0;
  1149.                     double meanPoleY = 0;
  1150.                     if (date.compareTo(lastAnnualPoleDate) <= 0) {
  1151.                         // we are within the range covered by the annual pole file,
  1152.                         // we interpolate within it
  1153.                         try {
  1154.                             final HermiteInterpolator interpolator = new HermiteInterpolator();
  1155.                             annualCache.getNeighbors(date).forEach(neighbor ->
  1156.                                 interpolator.addSamplePoint(neighbor.getDate().durationFrom(date),
  1157.                                                             new double[] {
  1158.                                                                 neighbor.getX(), neighbor.getY()
  1159.                                                             }));
  1160.                             final double[] interpolated = interpolator.value(0);
  1161.                             meanPoleX = interpolated[0];
  1162.                             meanPoleY = interpolated[1];
  1163.                         } catch (TimeStampedCacheException tsce) {
  1164.                             // this should never happen
  1165.                             throw new OrekitInternalError(tsce);
  1166.                         }
  1167.                     } else {

  1168.                         // we are after the range covered by the annual pole file,
  1169.                         // we use the polynomial extension
  1170.                         final double t = date.durationFrom(
  1171.                                 eopHistory.getTimeScales().getJ2000Epoch());
  1172.                         meanPoleX = xp0 + t * xp0Dot;
  1173.                         meanPoleY = yp0 + t * yp0Dot;

  1174.                     }

  1175.                     // evaluate wobble variables
  1176.                     final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1177.                     final double m1 = correction.getXp() - meanPoleX;
  1178.                     final double m2 = meanPoleY - correction.getYp();

  1179.                     return new double[] {
  1180.                         // the following correspond to the equations published in IERS 2003 conventions,
  1181.                         // section 6.2 page 65. In the publication, the equations read:
  1182.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1183.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1184.                         // However, it seems there are sign errors in these equations, which have
  1185.                         // been fixed in IERS 2010 conventions, section 6.4 page 94. In these newer
  1186.                         // publication, the equations read:
  1187.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1188.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1189.                         // the newer equations seem more consistent with the premises as the
  1190.                         // deformation due to the centrifugal potential has the form:
  1191.                         // −Ω²r²/2 sin 2θ Re [k₂(m₁ − im₂) exp(iλ)] where k₂ is the complex
  1192.                         // number 0.3077 + 0.0036i, so the real part in the previous equation is:
  1193.                         // A[Re(k₂) m₁ + Im(k₂) m₂)] cos λ + A[Re(k₂) m₂ - Im(k₂) m₁] sin λ
  1194.                         // identifying this with ∆C₂₁ cos λ + ∆S₂₁ sin λ we get:
  1195.                         // ∆C₂₁ = A Re(k₂) [m₁ + Im(k₂)/Re(k₂) m₂)]
  1196.                         // ∆S₂₁ = A Re(k₂) [m₂ - Im(k₂)/Re(k₂) m₁)]
  1197.                         // and Im(k₂)/Re(k₂) is very close to +0.00115
  1198.                         // As the equation as written in the IERS 2003 conventions are used in
  1199.                         // legacy systems, we have reproduced this alleged error here (and fixed it in
  1200.                         // the IERS 2010 conventions below) for validation purposes. We don't recommend
  1201.                         // using the IERS 2003 conventions for solid pole tide computation other than
  1202.                         // for validation or reproducibility of legacy applications behavior.
  1203.                         // As solid pole tide is small and as the sign change is on the smallest coefficient,
  1204.                         // the effect is quite small. A test case on a propagated orbit showed a position change
  1205.                         // slightly below 0.4m after a 30 days propagation on a Low Earth Orbit
  1206.                         globalFactor * (m1 - ratio * m2),
  1207.                         globalFactor * (m2 + ratio * m1),
  1208.                     };

  1209.                 }

  1210.                 /** {@inheritDoc} */
  1211.                 @Override
  1212.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

  1213.                     final AbsoluteDate aDate = date.toAbsoluteDate();

  1214.                     // we can't compute anything before the range covered by the annual pole file
  1215.                     if (aDate.compareTo(firstAnnualPoleDate) <= 0) {
  1216.                         return MathArrays.buildArray(date.getField(), 2);
  1217.                     }

  1218.                     // evaluate mean pole
  1219.                     T meanPoleX = date.getField().getZero();
  1220.                     T meanPoleY = date.getField().getZero();
  1221.                     if (aDate.compareTo(lastAnnualPoleDate) <= 0) {
  1222.                         // we are within the range covered by the annual pole file,
  1223.                         // we interpolate within it
  1224.                         try {
  1225.                             final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
  1226.                             final T[] y = MathArrays.buildArray(date.getField(), 2);
  1227.                             final T zero = date.getField().getZero();
  1228.                             final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
  1229.                                                                                                        // for example removing derivatives
  1230.                                                                                                        // if T was DerivativeStructure
  1231.                             annualCache.getNeighbors(aDate).forEach(neighbor -> {
  1232.                                 y[0] = zero.add(neighbor.getX());
  1233.                                 y[1] = zero.add(neighbor.getY());
  1234.                                 interpolator.addSamplePoint(central.durationFrom(neighbor.getDate()).negate(), y);
  1235.                             });
  1236.                             final T[] interpolated = interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
  1237.                             meanPoleX = interpolated[0];
  1238.                             meanPoleY = interpolated[1];
  1239.                         } catch (TimeStampedCacheException tsce) {
  1240.                             // this should never happen
  1241.                             throw new OrekitInternalError(tsce);
  1242.                         }
  1243.                     } else {

  1244.                         // we are after the range covered by the annual pole file,
  1245.                         // we use the polynomial extension
  1246.                         final T t = date.durationFrom(
  1247.                                 eopHistory.getTimeScales().getJ2000Epoch());
  1248.                         meanPoleX = t.multiply(xp0Dot).add(xp0);
  1249.                         meanPoleY = t.multiply(yp0Dot).add(yp0);

  1250.                     }

  1251.                     // evaluate wobble variables
  1252.                     final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1253.                     final T m1 = correction.getXp().subtract(meanPoleX);
  1254.                     final T m2 = meanPoleY.subtract(correction.getYp());

  1255.                     final T[] a = MathArrays.buildArray(date.getField(), 2);

  1256.                     // the following correspond to the equations published in IERS 2003 conventions,
  1257.                     // section 6.2 page 65. In the publication, the equations read:
  1258.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1259.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1260.                     // However, it seems there are sign errors in these equations, which have
  1261.                     // been fixed in IERS 2010 conventions, section 6.4 page 94. In these newer
  1262.                     // publication, the equations read:
  1263.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1264.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1265.                     // the newer equations seem more consistent with the premises as the
  1266.                     // deformation due to the centrifugal potential has the form:
  1267.                     // −Ω²r²/2 sin 2θ Re [k₂(m₁ − im₂) exp(iλ)] where k₂ is the complex
  1268.                     // number 0.3077 + 0.0036i, so the real part in the previous equation is:
  1269.                     // A[Re(k₂) m₁ + Im(k₂) m₂)] cos λ + A[Re(k₂) m₂ - Im(k₂) m₁] sin λ
  1270.                     // identifying this with ∆C₂₁ cos λ + ∆S₂₁ sin λ we get:
  1271.                     // ∆C₂₁ = A Re(k₂) [m₁ + Im(k₂)/Re(k₂) m₂)]
  1272.                     // ∆S₂₁ = A Re(k₂) [m₂ - Im(k₂)/Re(k₂) m₁)]
  1273.                     // and Im(k₂)/Re(k₂) is very close to +0.00115
  1274.                     // As the equation as written in the IERS 2003 conventions are used in
  1275.                     // legacy systems, we have reproduced this alleged error here (and fixed it in
  1276.                     // the IERS 2010 conventions below) for validation purposes. We don't recommend
  1277.                     // using the IERS 2003 conventions for solid pole tide computation other than
  1278.                     // for validation or reproducibility of legacy applications behavior.
  1279.                     // As solid pole tide is small and as the sign change is on the smallest coefficient,
  1280.                     // the effect is quite small. A test case on a propagated orbit showed a position change
  1281.                     // slightly below 0.4m after a 30 days propagation on a Low Earth Orbit
  1282.                     a[0] = m1.add(m2.multiply(-ratio)).multiply(globalFactor);
  1283.                     a[1] = m2.add(m1.multiply( ratio)).multiply(globalFactor);

  1284.                     return a;

  1285.                 }

  1286.             };

  1287.         }

  1288.         /** {@inheritDoc} */
  1289.         @Override
  1290.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {

  1291.             return new TimeVectorFunction() {

  1292.                 /** {@inheritDoc} */
  1293.                 @Override
  1294.                 public double[] value(final AbsoluteDate date) {
  1295.                     // there are no model for ocean pole tide prior to conventions 2010
  1296.                     return new double[] {
  1297.                         0.0, 0.0
  1298.                     };
  1299.                 }

  1300.                 /** {@inheritDoc} */
  1301.                 @Override
  1302.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1303.                     // there are no model for ocean pole tide prior to conventions 2010
  1304.                     return MathArrays.buildArray(date.getField(), 2);
  1305.                 }

  1306.             };
  1307.         }

  1308.         /** {@inheritDoc} */
  1309.         @Override
  1310.         public double[] getNominalTidalDisplacement() {
  1311.             return new double[] {
  1312.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  1313.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  1314.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  1315.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  1316.                 // H₀
  1317.                 -0.31460
  1318.             };
  1319.         }

  1320.         /** {@inheritDoc} */
  1321.         @Override
  1322.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
  1323.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  1324.                                                                   18, 15, 16, 17, 18);
  1325.         }

  1326.         /** {@inheritDoc} */
  1327.         @Override
  1328.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  1329.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  1330.                                                                 18, 15, 16, 17, 18);
  1331.         }

  1332.     },

  1333.     /** Constant for IERS 2010 conventions. */
  1334.     IERS_2010 {

  1335.         /** Nutation arguments resources. */
  1336.         private static final String NUTATION_ARGUMENTS = IERS_BASE + "2010/nutation-arguments.txt";

  1337.         /** X series resources. */
  1338.         private static final String X_SERIES           = IERS_BASE + "2010/tab5.2a.txt";

  1339.         /** Y series resources. */
  1340.         private static final String Y_SERIES           = IERS_BASE + "2010/tab5.2b.txt";

  1341.         /** S series resources. */
  1342.         private static final String S_SERIES           = IERS_BASE + "2010/tab5.2d.txt";

  1343.         /** Psi series resources. */
  1344.         private static final String PSI_SERIES         = IERS_BASE + "2010/tab5.3a.txt";

  1345.         /** Epsilon series resources. */
  1346.         private static final String EPSILON_SERIES     = IERS_BASE + "2010/tab5.3b.txt";

  1347.         /** Greenwhich sidereal time series resources. */
  1348.         private static final String GST_SERIES         = IERS_BASE + "2010/tab5.2e.txt";

  1349.         /** Tidal correction for xp, yp series resources. */
  1350.         private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "2010/tab8.2ab.txt";

  1351.         /** Tidal correction for UT1 resources. */
  1352.         private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "2010/tab8.3ab.txt";

  1353.         /** Love numbers resources. */
  1354.         private static final String LOVE_NUMBERS = IERS_BASE + "2010/tab6.3.txt";

  1355.         /** Frequency dependence model for k₂₀. */
  1356.         private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5b.txt";

  1357.         /** Frequency dependence model for k₂₁. */
  1358.         private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5a.txt";

  1359.         /** Frequency dependence model for k₂₂. */
  1360.         private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5c.txt";

  1361.         /** Tidal displacement frequency correction for diurnal tides. */
  1362.         private static final String TIDAL_DISPLACEMENT_CORRECTION_DIURNAL = IERS_BASE + "2010/tab7.3a.txt";

  1363.         /** Tidal displacement frequency correction for zonal tides. */
  1364.         private static final String TIDAL_DISPLACEMENT_CORRECTION_ZONAL = IERS_BASE + "2010/tab7.3b.txt";

  1365.         /** {@inheritDoc} */
  1366.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
  1367.                                                                  final TimeScales timeScales) {
  1368.             try (InputStream in = getStream(NUTATION_ARGUMENTS)) {
  1369.                 return new FundamentalNutationArguments(this, timeScale,
  1370.                         in, NUTATION_ARGUMENTS, timeScales);
  1371.             } catch (IOException e) {
  1372.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1373.             }
  1374.         }

  1375.         /** {@inheritDoc} */
  1376.         @Override
  1377.         public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {

  1378.             // epsilon 0 value from chapter 5, page 56, other terms from equation 5.40 page 65
  1379.             final PolynomialNutation epsilonA =
  1380.                     new PolynomialNutation(84381.406        * Constants.ARC_SECONDS_TO_RADIANS,
  1381.                                              -46.836769     * Constants.ARC_SECONDS_TO_RADIANS,
  1382.                                               -0.0001831    * Constants.ARC_SECONDS_TO_RADIANS,
  1383.                                                0.00200340   * Constants.ARC_SECONDS_TO_RADIANS,
  1384.                                               -0.000000576  * Constants.ARC_SECONDS_TO_RADIANS,
  1385.                                               -0.0000000434 * Constants.ARC_SECONDS_TO_RADIANS);

  1386.             return new TimeScalarFunction() {

  1387.                 /** {@inheritDoc} */
  1388.                 @Override
  1389.                 public double value(final AbsoluteDate date) {
  1390.                     return epsilonA.value(evaluateTC(date, timeScales));
  1391.                 }

  1392.                 /** {@inheritDoc} */
  1393.                 @Override
  1394.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1395.                     return epsilonA.value(evaluateTC(date, timeScales));
  1396.                 }

  1397.             };

  1398.         }

  1399.         /** {@inheritDoc} */
  1400.         @Override
  1401.         public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {

  1402.             // set up nutation arguments
  1403.             final FundamentalNutationArguments arguments =
  1404.                     getNutationArguments(timeScales);

  1405.             // set up Poisson series
  1406.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1407.             final PoissonSeriesParser parser =
  1408.                     new PoissonSeriesParser(17).
  1409.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  1410.                         withFirstDelaunay(4).
  1411.                         withFirstPlanetary(9).
  1412.                         withSinCos(0, 2, microAS, 3, microAS);
  1413.             final PoissonSeries.CompiledSeries xys;
  1414.             try (InputStream xIn = getStream(X_SERIES);
  1415.                  InputStream yIn = getStream(Y_SERIES);
  1416.                  InputStream sIn = getStream(S_SERIES)) {
  1417.                 final PoissonSeries xSeries = parser.parse(xIn, X_SERIES);
  1418.                 final PoissonSeries ySeries = parser.parse(yIn, Y_SERIES);
  1419.                 final PoissonSeries sSeries = parser.parse(sIn, S_SERIES);
  1420.                 xys = PoissonSeries.compile(xSeries, ySeries, sSeries);
  1421.             } catch (IOException e) {
  1422.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1423.             }

  1424.             // create a function evaluating the series
  1425.             return new TimeVectorFunction() {

  1426.                 /** {@inheritDoc} */
  1427.                 @Override
  1428.                 public double[] value(final AbsoluteDate date) {
  1429.                     return xys.value(arguments.evaluateAll(date));
  1430.                 }

  1431.                 /** {@inheritDoc} */
  1432.                 @Override
  1433.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1434.                     return xys.value(arguments.evaluateAll(date));
  1435.                 }

  1436.             };

  1437.         }

  1438.         /** {@inheritDoc} */
  1439.         public LoveNumbers getLoveNumbers() {
  1440.             return loadLoveNumbers(LOVE_NUMBERS);
  1441.         }

  1442.         /** {@inheritDoc} */
  1443.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  1444.                                                                      final TimeScales timeScales) {

  1445.             // set up nutation arguments
  1446.             final FundamentalNutationArguments arguments = getNutationArguments(ut1, timeScales);

  1447.             // set up Poisson series
  1448.             final PoissonSeriesParser k20Parser =
  1449.                     new PoissonSeriesParser(18).
  1450.                         withOptionalColumn(1).
  1451.                         withDoodson(4, 2).
  1452.                         withFirstDelaunay(10);
  1453.             final PoissonSeriesParser k21Parser =
  1454.                     new PoissonSeriesParser(18).
  1455.                         withOptionalColumn(1).
  1456.                         withDoodson(4, 3).
  1457.                         withFirstDelaunay(10);
  1458.             final PoissonSeriesParser k22Parser =
  1459.                     new PoissonSeriesParser(16).
  1460.                         withOptionalColumn(1).
  1461.                         withDoodson(4, 2).
  1462.                         withFirstDelaunay(10);

  1463.             final double pico = 1.0e-12;
  1464.             try (InputStream k0In = getStream(K20_FREQUENCY_DEPENDENCE);
  1465.                  InputStream k21In1 = getStream(K21_FREQUENCY_DEPENDENCE);
  1466.                  InputStream k21In2 = getStream(K21_FREQUENCY_DEPENDENCE);
  1467.                  InputStream k22In1 = getStream(K22_FREQUENCY_DEPENDENCE);
  1468.                  InputStream k22In2 = getStream(K22_FREQUENCY_DEPENDENCE)) {
  1469.                 final PoissonSeries c20Series =
  1470.                         k20Parser.
  1471.                         withSinCos(0, 18, -pico, 16, pico).
  1472.                         parse(k0In, K20_FREQUENCY_DEPENDENCE);
  1473.                 final PoissonSeries c21Series =
  1474.                         k21Parser.
  1475.                         withSinCos(0, 17, pico, 18, pico).
  1476.                         parse(k21In1, K21_FREQUENCY_DEPENDENCE);
  1477.                 final PoissonSeries s21Series =
  1478.                         k21Parser.
  1479.                         withSinCos(0, 18, -pico, 17, pico).
  1480.                         parse(k21In2, K21_FREQUENCY_DEPENDENCE);
  1481.                 final PoissonSeries c22Series =
  1482.                         k22Parser.
  1483.                         withSinCos(0, -1, pico, 16, pico).
  1484.                         parse(k22In1, K22_FREQUENCY_DEPENDENCE);
  1485.                 final PoissonSeries s22Series =
  1486.                         k22Parser.
  1487.                         withSinCos(0, 16, -pico, -1, pico).
  1488.                         parse(k22In2, K22_FREQUENCY_DEPENDENCE);

  1489.                 return new TideFrequencyDependenceFunction(arguments,
  1490.                                                            c20Series,
  1491.                                                            c21Series, s21Series,
  1492.                                                            c22Series, s22Series);
  1493.             } catch (IOException e) {
  1494.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1495.             }

  1496.         }

  1497.         /** {@inheritDoc} */
  1498.         @Override
  1499.         public double getPermanentTide() {
  1500.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1501.         }

  1502.         /** Compute pole wobble variables m₁ and m₂.
  1503.          * @param date current date
  1504.          * @param eopHistory EOP history
  1505.          * @return array containing m₁ and m₂
  1506.          */
  1507.         private double[] computePoleWobble(final AbsoluteDate date, final EOPHistory eopHistory) {

  1508.             // polynomial model from IERS 2010, table 7.7
  1509.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1510.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1511.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1512.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1513.             final AbsoluteDate changeDate =
  1514.                     new AbsoluteDate(2010, 1, 1, eopHistory.getTimeScales().getTT());

  1515.             // evaluate mean pole
  1516.             final double[] xPolynomial;
  1517.             final double[] yPolynomial;
  1518.             if (date.compareTo(changeDate) <= 0) {
  1519.                 xPolynomial = new double[] {
  1520.                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
  1521.                 };
  1522.                 yPolynomial = new double[] {
  1523.                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
  1524.                 };
  1525.             } else {
  1526.                 xPolynomial = new double[] {
  1527.                     23.513 * f0, 7.6141 * f1
  1528.                 };
  1529.                 yPolynomial = new double[] {
  1530.                     358.891 * f0,  -0.6287 * f1
  1531.                 };
  1532.             }
  1533.             double meanPoleX = 0;
  1534.             double meanPoleY = 0;
  1535.             final double t = date.durationFrom(
  1536.                     eopHistory.getTimeScales().getJ2000Epoch());
  1537.             for (int i = xPolynomial.length - 1; i >= 0; --i) {
  1538.                 meanPoleX = meanPoleX * t + xPolynomial[i];
  1539.             }
  1540.             for (int i = yPolynomial.length - 1; i >= 0; --i) {
  1541.                 meanPoleY = meanPoleY * t + yPolynomial[i];
  1542.             }

  1543.             // evaluate wobble variables
  1544.             final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1545.             final double m1 = correction.getXp() - meanPoleX;
  1546.             final double m2 = meanPoleY - correction.getYp();

  1547.             return new double[] {
  1548.                 m1, m2
  1549.             };

  1550.         }

  1551.         /** Compute pole wobble variables m₁ and m₂.
  1552.          * @param date current date
  1553.          * @param <T> type of the field elements
  1554.          * @param eopHistory EOP history
  1555.          * @return array containing m₁ and m₂
  1556.          */
  1557.         private <T extends RealFieldElement<T>> T[] computePoleWobble(final FieldAbsoluteDate<T> date, final EOPHistory eopHistory) {

  1558.             // polynomial model from IERS 2010, table 7.7
  1559.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1560.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1561.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1562.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1563.             final AbsoluteDate changeDate =
  1564.                     new AbsoluteDate(2010, 1, 1, eopHistory.getTimeScales().getTT());

  1565.             // evaluate mean pole
  1566.             final double[] xPolynomial;
  1567.             final double[] yPolynomial;
  1568.             if (date.toAbsoluteDate().compareTo(changeDate) <= 0) {
  1569.                 xPolynomial = new double[] {
  1570.                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
  1571.                 };
  1572.                 yPolynomial = new double[] {
  1573.                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
  1574.                 };
  1575.             } else {
  1576.                 xPolynomial = new double[] {
  1577.                     23.513 * f0, 7.6141 * f1
  1578.                 };
  1579.                 yPolynomial = new double[] {
  1580.                     358.891 * f0,  -0.6287 * f1
  1581.                 };
  1582.             }
  1583.             T meanPoleX = date.getField().getZero();
  1584.             T meanPoleY = date.getField().getZero();
  1585.             final T t = date.durationFrom(
  1586.                     eopHistory.getTimeScales().getJ2000Epoch());
  1587.             for (int i = xPolynomial.length - 1; i >= 0; --i) {
  1588.                 meanPoleX = meanPoleX.multiply(t).add(xPolynomial[i]);
  1589.             }
  1590.             for (int i = yPolynomial.length - 1; i >= 0; --i) {
  1591.                 meanPoleY = meanPoleY.multiply(t).add(yPolynomial[i]);
  1592.             }

  1593.             // evaluate wobble variables
  1594.             final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1595.             final T[] m = MathArrays.buildArray(date.getField(), 2);
  1596.             m[0] = correction.getXp().subtract(meanPoleX);
  1597.             m[1] = meanPoleY.subtract(correction.getYp());

  1598.             return m;

  1599.         }

  1600.         /** {@inheritDoc} */
  1601.         @Override
  1602.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  1603.             // constants from IERS 2010, section 6.4
  1604.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1605.             final double ratio        =  0.00115;

  1606.             return new TimeVectorFunction() {

  1607.                 /** {@inheritDoc} */
  1608.                 @Override
  1609.                 public double[] value(final AbsoluteDate date) {

  1610.                     // evaluate wobble variables
  1611.                     final double[] wobbleM = computePoleWobble(date, eopHistory);

  1612.                     return new double[] {
  1613.                         // the following correspond to the equations published in IERS 2010 conventions,
  1614.                         // section 6.4 page 94. The equations read:
  1615.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1616.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1617.                         // These equations seem to fix what was probably a sign error in IERS 2003
  1618.                         // conventions section 6.2 page 65. In this older publication, the equations read:
  1619.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1620.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1621.                         globalFactor * (wobbleM[0] + ratio * wobbleM[1]),
  1622.                         globalFactor * (wobbleM[1] - ratio * wobbleM[0])
  1623.                     };

  1624.                 }

  1625.                 /** {@inheritDoc} */
  1626.                 @Override
  1627.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

  1628.                     // evaluate wobble variables
  1629.                     final T[] wobbleM = computePoleWobble(date, eopHistory);

  1630.                     final T[] a = MathArrays.buildArray(date.getField(), 2);

  1631.                     // the following correspond to the equations published in IERS 2010 conventions,
  1632.                     // section 6.4 page 94. The equations read:
  1633.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1634.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1635.                     // These equations seem to fix what was probably a sign error in IERS 2003
  1636.                     // conventions section 6.2 page 65. In this older publication, the equations read:
  1637.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1638.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1639.                     a[0] = wobbleM[0].add(wobbleM[1].multiply( ratio)).multiply(globalFactor);
  1640.                     a[1] = wobbleM[1].add(wobbleM[0].multiply(-ratio)).multiply(globalFactor);

  1641.                     return a;

  1642.                 }

  1643.             };

  1644.         }

  1645.         /** {@inheritDoc} */
  1646.         @Override
  1647.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {

  1648.             return new TimeVectorFunction() {

  1649.                 /** {@inheritDoc} */
  1650.                 @Override
  1651.                 public double[] value(final AbsoluteDate date) {

  1652.                     // evaluate wobble variables
  1653.                     final double[] wobbleM = computePoleWobble(date, eopHistory);

  1654.                     return new double[] {
  1655.                         // the following correspond to the equations published in IERS 2010 conventions,
  1656.                         // section 6.4 page 94 equation 6.24:
  1657.                         // ∆C₂₁ = −2.1778 × 10⁻¹⁰ (m₁ − 0.01724m₂)
  1658.                         // ∆S₂₁ = −1.7232 × 10⁻¹⁰ (m₂ − 0.03365m₁)
  1659.                         -2.1778e-10 * (wobbleM[0] - 0.01724 * wobbleM[1]) / Constants.ARC_SECONDS_TO_RADIANS,
  1660.                         -1.7232e-10 * (wobbleM[1] - 0.03365 * wobbleM[0]) / Constants.ARC_SECONDS_TO_RADIANS
  1661.                     };

  1662.                 }

  1663.                 /** {@inheritDoc} */
  1664.                 @Override
  1665.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

  1666.                     final T[] v = MathArrays.buildArray(date.getField(), 2);

  1667.                     // evaluate wobble variables
  1668.                     final T[] wobbleM = computePoleWobble(date, eopHistory);

  1669.                     // the following correspond to the equations published in IERS 2010 conventions,
  1670.                     // section 6.4 page 94 equation 6.24:
  1671.                     // ∆C₂₁ = −2.1778 × 10⁻¹⁰ (m₁ − 0.01724m₂)
  1672.                     // ∆S₂₁ = −1.7232 × 10⁻¹⁰ (m₂ − 0.03365m₁)
  1673.                     v[0] = wobbleM[0].subtract(wobbleM[1].multiply(0.01724)).multiply(-2.1778e-10 / Constants.ARC_SECONDS_TO_RADIANS);
  1674.                     v[1] = wobbleM[1].subtract(wobbleM[0].multiply(0.03365)).multiply(-1.7232e-10 / Constants.ARC_SECONDS_TO_RADIANS);

  1675.                     return v;

  1676.                 }

  1677.             };

  1678.         }

  1679.         /** {@inheritDoc} */
  1680.         @Override
  1681.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

  1682.             // set up the conventional polynomials
  1683.             // the following values are from equation 5.40 in IERS 2010 conventions
  1684.             final PolynomialNutation psiA =
  1685.                     new PolynomialNutation(   0.0,
  1686.                                            5038.481507     * Constants.ARC_SECONDS_TO_RADIANS,
  1687.                                              -1.0790069    * Constants.ARC_SECONDS_TO_RADIANS,
  1688.                                              -0.00114045   * Constants.ARC_SECONDS_TO_RADIANS,
  1689.                                               0.000132851  * Constants.ARC_SECONDS_TO_RADIANS,
  1690.                                              -0.0000000951 * Constants.ARC_SECONDS_TO_RADIANS);
  1691.             final PolynomialNutation omegaA =
  1692.                     new PolynomialNutation(getMeanObliquityFunction(timeScales)
  1693.                             .value(getNutationReferenceEpoch(timeScales)),
  1694.                                            -0.025754     * Constants.ARC_SECONDS_TO_RADIANS,
  1695.                                             0.0512623    * Constants.ARC_SECONDS_TO_RADIANS,
  1696.                                            -0.00772503   * Constants.ARC_SECONDS_TO_RADIANS,
  1697.                                            -0.000000467  * Constants.ARC_SECONDS_TO_RADIANS,
  1698.                                             0.0000003337 * Constants.ARC_SECONDS_TO_RADIANS);
  1699.             final PolynomialNutation chiA =
  1700.                     new PolynomialNutation( 0.0,
  1701.                                            10.556403     * Constants.ARC_SECONDS_TO_RADIANS,
  1702.                                            -2.3814292    * Constants.ARC_SECONDS_TO_RADIANS,
  1703.                                            -0.00121197   * Constants.ARC_SECONDS_TO_RADIANS,
  1704.                                             0.000170663  * Constants.ARC_SECONDS_TO_RADIANS,
  1705.                                            -0.0000000560 * Constants.ARC_SECONDS_TO_RADIANS);

  1706.             return new PrecessionFunction(psiA, omegaA, chiA, timeScales);

  1707.         }

  1708.         /** {@inheritDoc} */
  1709.         @Override
  1710.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  1711.             // set up nutation arguments
  1712.             final FundamentalNutationArguments arguments =
  1713.                     getNutationArguments(timeScales);

  1714.             // set up Poisson series
  1715.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1716.             final PoissonSeriesParser parser =
  1717.                     new PoissonSeriesParser(17).
  1718.                         withFirstDelaunay(4).
  1719.                         withFirstPlanetary(9).
  1720.                         withSinCos(0, 2, microAS, 3, microAS);
  1721.             final PoissonSeries.CompiledSeries psiEpsilonSeries;
  1722.             try (InputStream psiIn = getStream(PSI_SERIES);
  1723.                  InputStream epsilonIn = getStream(EPSILON_SERIES)) {
  1724.                 final PoissonSeries psiSeries     = parser.parse(psiIn, PSI_SERIES);
  1725.                 final PoissonSeries epsilonSeries = parser.parse(epsilonIn, EPSILON_SERIES);
  1726.                 psiEpsilonSeries = PoissonSeries.compile(psiSeries, epsilonSeries);
  1727.             } catch (IOException e) {
  1728.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1729.             }

  1730.             return new TimeVectorFunction() {

  1731.                 /** {@inheritDoc} */
  1732.                 @Override
  1733.                 public double[] value(final AbsoluteDate date) {
  1734.                     final BodiesElements elements = arguments.evaluateAll(date);
  1735.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  1736.                     return new double[] {
  1737.                         psiEpsilon[0], psiEpsilon[1],
  1738.                         IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  1739.                     };
  1740.                 }

  1741.                 /** {@inheritDoc} */
  1742.                 @Override
  1743.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1744.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1745.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  1746.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  1747.                     result[0] = psiEpsilon[0];
  1748.                     result[1] = psiEpsilon[1];
  1749.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  1750.                     return result;
  1751.                 }

  1752.             };

  1753.         }

  1754.         /** {@inheritDoc} */
  1755.         @Override
  1756.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  1757.                                                   final TimeScales timeScales) {

  1758.             // Earth Rotation Angle
  1759.             final StellarAngleCapitaine era =
  1760.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  1761.             // Polynomial part of the apparent sidereal time series
  1762.             // which is the opposite of Equation of Origins (EO)
  1763.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1764.             final PoissonSeriesParser parser =
  1765.                     new PoissonSeriesParser(17).
  1766.                         withFirstDelaunay(4).
  1767.                         withFirstPlanetary(9).
  1768.                         withSinCos(0, 2, microAS, 3, microAS).
  1769.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1770.             final PolynomialNutation minusEO;
  1771.             try (InputStream in = getStream(GST_SERIES)) {
  1772.                 minusEO = parser.parse(in, GST_SERIES).getPolynomial();
  1773.             } catch (IOException e) {
  1774.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1775.             }

  1776.             // create a function evaluating the series
  1777.             return new TimeScalarFunction() {

  1778.                 /** {@inheritDoc} */
  1779.                 @Override
  1780.                 public double value(final AbsoluteDate date) {
  1781.                     return era.value(date) + minusEO.value(evaluateTC(date, timeScales));
  1782.                 }

  1783.                 /** {@inheritDoc} */
  1784.                 @Override
  1785.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1786.                     return era.value(date).add(minusEO.value(evaluateTC(date, timeScales)));
  1787.                 }

  1788.             };

  1789.         }

  1790.         /** {@inheritDoc} */
  1791.         @Override
  1792.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  1793.                                                       final TimeScales timeScales) {

  1794.             // Earth Rotation Angle
  1795.             final StellarAngleCapitaine era =
  1796.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  1797.             // Polynomial part of the apparent sidereal time series
  1798.             // which is the opposite of Equation of Origins (EO)
  1799.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1800.             final PoissonSeriesParser parser =
  1801.                     new PoissonSeriesParser(17).
  1802.                         withFirstDelaunay(4).
  1803.                         withFirstPlanetary(9).
  1804.                         withSinCos(0, 2, microAS, 3, microAS).
  1805.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1806.             final PolynomialNutation minusEO;
  1807.             try (InputStream in = getStream(GST_SERIES)) {
  1808.                 minusEO = parser.parse(in, GST_SERIES).getPolynomial();
  1809.             } catch (IOException e) {
  1810.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1811.             }

  1812.             // create a function evaluating the series
  1813.             return new TimeScalarFunction() {

  1814.                 /** {@inheritDoc} */
  1815.                 @Override
  1816.                 public double value(final AbsoluteDate date) {
  1817.                     return era.getRate() + minusEO.derivative(evaluateTC(date, timeScales));
  1818.                 }

  1819.                 /** {@inheritDoc} */
  1820.                 @Override
  1821.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1822.                     return minusEO.derivative(evaluateTC(date, timeScales)).add(era.getRate());
  1823.                 }

  1824.             };

  1825.         }

  1826.         /** {@inheritDoc} */
  1827.         @Override
  1828.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  1829.                                                   final EOPHistory eopHistory,
  1830.                                                   final TimeScales timeScales) {

  1831.             // set up nutation arguments
  1832.             final FundamentalNutationArguments arguments =
  1833.                     getNutationArguments(timeScales);

  1834.             // mean obliquity function
  1835.             final TimeScalarFunction epsilon = getMeanObliquityFunction(timeScales);

  1836.             // set up Poisson series
  1837.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1838.             final PoissonSeriesParser baseParser =
  1839.                     new PoissonSeriesParser(17).
  1840.                         withFirstDelaunay(4).
  1841.                         withFirstPlanetary(9).
  1842.                         withSinCos(0, 2, microAS, 3, microAS);
  1843.             final PoissonSeriesParser gstParser  = baseParser.withPolynomialPart('t', Unit.ARC_SECONDS);
  1844.             final PoissonSeries.CompiledSeries psiGstSeries;
  1845.             try (InputStream psiIn = getStream(PSI_SERIES);
  1846.                  InputStream gstIn = getStream(GST_SERIES)) {
  1847.                 final PoissonSeries psiSeries        = baseParser.parse(psiIn, PSI_SERIES);
  1848.                 final PoissonSeries gstSeries        = gstParser.parse(gstIn, GST_SERIES);
  1849.                 psiGstSeries = PoissonSeries.compile(psiSeries, gstSeries);
  1850.             } catch (IOException e) {
  1851.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1852.             }

  1853.             // ERA function
  1854.             final TimeScalarFunction era =
  1855.                     getEarthOrientationAngleFunction(ut1, timeScales.getTAI());

  1856.             return new TimeScalarFunction() {

  1857.                 /** {@inheritDoc} */
  1858.                 @Override
  1859.                 public double value(final AbsoluteDate date) {

  1860.                     // evaluate equation of origins
  1861.                     final BodiesElements elements = arguments.evaluateAll(date);
  1862.                     final double[] angles = psiGstSeries.value(elements);
  1863.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  1864.                     final double deltaPsi = angles[0] + ddPsi;
  1865.                     final double epsilonA = epsilon.value(date);

  1866.                     // subtract equation of origin from EA
  1867.                     // (hence add the series above which have the sign included)
  1868.                     return era.value(date) + deltaPsi * FastMath.cos(epsilonA) + angles[1];

  1869.                 }

  1870.                 /** {@inheritDoc} */
  1871.                 @Override
  1872.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  1873.                     // evaluate equation of origins
  1874.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1875.                     final T[] angles = psiGstSeries.value(elements);
  1876.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  1877.                     final T deltaPsi = angles[0].add(ddPsi);
  1878.                     final T epsilonA = epsilon.value(date);

  1879.                     // subtract equation of origin from EA
  1880.                     // (hence add the series above which have the sign included)
  1881.                     return era.value(date).add(deltaPsi.multiply(epsilonA.cos())).add(angles[1]);

  1882.                 }

  1883.             };

  1884.         }

  1885.         /** {@inheritDoc} */
  1886.         @Override
  1887.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

  1888.             // set up nutation arguments
  1889.             // BEWARE! Using TT as the time scale here and not UT1 is intentional!
  1890.             // as this correction is used to compute UT1 itself, it is not surprising we cannot use UT1 yet,
  1891.             // however, using the close UTC as would seem logical make the comparison with interp.f from IERS fail
  1892.             // looking in the interp.f code, the same TT scale is used for both Delaunay and gamma argument
  1893.             final FundamentalNutationArguments arguments =
  1894.                     getNutationArguments(timeScales.getTT(), timeScales);

  1895.             // set up Poisson series
  1896.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1897.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
  1898.                     withOptionalColumn(1).
  1899.                     withGamma(2).
  1900.                     withFirstDelaunay(3);
  1901.             try (InputStream xpIn = getStream(TIDAL_CORRECTION_XP_YP_SERIES);
  1902.                  InputStream ypIn = getStream(TIDAL_CORRECTION_XP_YP_SERIES);
  1903.                  InputStream ut1In = getStream(TIDAL_CORRECTION_UT1_SERIES)) {
  1904.                 final PoissonSeries xSeries =
  1905.                         xyParser.
  1906.                         withSinCos(0, 10, microAS, 11, microAS).
  1907.                         parse(xpIn, TIDAL_CORRECTION_XP_YP_SERIES);
  1908.                 final PoissonSeries ySeries =
  1909.                         xyParser.
  1910.                         withSinCos(0, 12, microAS, 13, microAS).
  1911.                         parse(ypIn, TIDAL_CORRECTION_XP_YP_SERIES);

  1912.                 final double microS = 1.0e-6;
  1913.                 final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  1914.                         withOptionalColumn(1).
  1915.                         withGamma(2).
  1916.                         withFirstDelaunay(3).
  1917.                         withSinCos(0, 10, microS, 11, microS);
  1918.                 final PoissonSeries ut1Series =
  1919.                         ut1Parser.parse(ut1In, TIDAL_CORRECTION_UT1_SERIES);

  1920.                 return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
  1921.             } catch (IOException e) {
  1922.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1923.             }

  1924.         }

  1925.         /** {@inheritDoc} */
  1926.         @Override
  1927.         public double[] getNominalTidalDisplacement() {
  1928.             return new double[] {
  1929.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  1930.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  1931.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  1932.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  1933.                 // H₀
  1934.                 -0.31460
  1935.             };
  1936.         }

  1937.         /** {@inheritDoc} */
  1938.         @Override
  1939.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
  1940.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  1941.                                                                   18, 15, 16, 17, 18);
  1942.         }

  1943.         /** {@inheritDoc} */
  1944.         @Override
  1945.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  1946.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  1947.                                                                 18, 15, 16, 17, 18);
  1948.         }

  1949.     };

  1950.     /** IERS conventions resources base directory. */
  1951.     private static final String IERS_BASE = "/assets/org/orekit/IERS-conventions/";

  1952.     /** Get the reference epoch for fundamental nutation arguments.
  1953.      *
  1954.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  1955.      *
  1956.      * @return reference epoch for fundamental nutation arguments
  1957.      * @since 6.1
  1958.      * @see #getNutationReferenceEpoch(TimeScales)
  1959.      */
  1960.     @DefaultDataContext
  1961.     public AbsoluteDate getNutationReferenceEpoch() {
  1962.         return getNutationReferenceEpoch(DataContext.getDefault().getTimeScales());
  1963.     }

  1964.     /**
  1965.      * Get the reference epoch for fundamental nutation arguments.
  1966.      *
  1967.      * @param timeScales to use for the reference epoch.
  1968.      * @return reference epoch for fundamental nutation arguments
  1969.      * @since 10.1
  1970.      */
  1971.     public AbsoluteDate getNutationReferenceEpoch(final TimeScales timeScales) {
  1972.         // IERS 1996, IERS 2003 and IERS 2010 use the same J2000.0 reference date
  1973.         return timeScales.getJ2000Epoch();
  1974.     }

  1975.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  1976.      *
  1977.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  1978.      *
  1979.      * @param date current date
  1980.      * @return date offset in Julian centuries
  1981.      * @since 6.1
  1982.      * @see #evaluateTC(AbsoluteDate, TimeScales)
  1983.      */
  1984.     @DefaultDataContext
  1985.     public double evaluateTC(final AbsoluteDate date) {
  1986.         return evaluateTC(date, DataContext.getDefault().getTimeScales());
  1987.     }

  1988.     /**
  1989.      * Evaluate the date offset between the current date and the {@link
  1990.      * #getNutationReferenceEpoch() reference date}.
  1991.      *
  1992.      * @param date       current date
  1993.      * @param timeScales used in the evaluation.
  1994.      * @return date offset in Julian centuries
  1995.      * @since 10.1
  1996.      */
  1997.     public double evaluateTC(final AbsoluteDate date, final TimeScales timeScales) {
  1998.         return date.durationFrom(getNutationReferenceEpoch(timeScales)) /
  1999.                 Constants.JULIAN_CENTURY;
  2000.     }

  2001.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2002.      *
  2003.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2004.      *
  2005.      * @param date current date
  2006.      * @param <T> type of the field elements
  2007.      * @return date offset in Julian centuries
  2008.      * @since 9.0
  2009.      * @see #evaluateTC(FieldAbsoluteDate, TimeScales)
  2010.      */
  2011.     @DefaultDataContext
  2012.     public <T extends RealFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date) {
  2013.         return evaluateTC(date, DataContext.getDefault().getTimeScales());
  2014.     }

  2015.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2016.      * @param <T> type of the field elements
  2017.      * @param date current date
  2018.      * @param timeScales used in the evaluation.
  2019.      * @return date offset in Julian centuries
  2020.      * @since 10.1
  2021.      */
  2022.     public <T extends RealFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date,
  2023.                                                         final TimeScales timeScales) {
  2024.         return date.durationFrom(getNutationReferenceEpoch(timeScales))
  2025.                 .divide(Constants.JULIAN_CENTURY);
  2026.     }

  2027.     /**
  2028.      * Get the fundamental nutation arguments. Does not compute GMST based values: gamma,
  2029.      * gammaDot.
  2030.      *
  2031.      * @param timeScales other time scales used in the computation including TAI and TT.
  2032.      * @return fundamental nutation arguments
  2033.      * @see #getNutationArguments(TimeScale, TimeScales)
  2034.      * @since 10.1
  2035.      */
  2036.     protected FundamentalNutationArguments getNutationArguments(
  2037.             final TimeScales timeScales) {

  2038.         return getNutationArguments(null, timeScales);
  2039.     }

  2040.     /** Get the fundamental nutation arguments.
  2041.      *
  2042.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2043.      *
  2044.      * @param timeScale time scale for computing Greenwich Mean Sidereal Time
  2045.      * (typically {@link TimeScales#getUT1(IERSConventions, boolean) UT1})
  2046.      * @return fundamental nutation arguments
  2047.      * @since 6.1
  2048.      * @see #getNutationArguments(TimeScale, TimeScales)
  2049.      * @see #getNutationArguments(TimeScales)
  2050.      */
  2051.     @DefaultDataContext
  2052.     public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale) {
  2053.         return getNutationArguments(timeScale, DataContext.getDefault().getTimeScales());
  2054.     }

  2055.     /**
  2056.      * Get the fundamental nutation arguments.
  2057.      *
  2058.      * @param timeScale time scale for computing Greenwich Mean Sidereal Time (typically
  2059.      *                  {@link TimeScales#getUT1(IERSConventions, boolean) UT1})
  2060.      * @param timeScales other time scales used in the computation including TAI and TT.
  2061.      * @return fundamental nutation arguments
  2062.      * @since 10.1
  2063.      */
  2064.     public abstract FundamentalNutationArguments getNutationArguments(
  2065.             TimeScale timeScale,
  2066.             TimeScales timeScales);

  2067.     /** Get the function computing mean obliquity of the ecliptic.
  2068.      *
  2069.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2070.      *
  2071.      * @return function computing mean obliquity of the ecliptic
  2072.      * @since 6.1
  2073.      * @see #getMeanObliquityFunction(TimeScales)
  2074.      */
  2075.     @DefaultDataContext
  2076.     public TimeScalarFunction getMeanObliquityFunction() {
  2077.         return getMeanObliquityFunction(DataContext.getDefault().getTimeScales());
  2078.     }

  2079.     /**
  2080.      * Get the function computing mean obliquity of the ecliptic.
  2081.      *
  2082.      * @param timeScales used in computing the function.
  2083.      * @return function computing mean obliquity of the ecliptic
  2084.      * @since 10.1
  2085.      */
  2086.     public abstract TimeScalarFunction getMeanObliquityFunction(TimeScales timeScales);

  2087.     /** Get the function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components.
  2088.      * <p>
  2089.      * The returned function computes the two X, Y components of CIP and the S+XY/2 component of the non-rotating CIO.
  2090.      * </p>
  2091.      *
  2092.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2093.      *
  2094.      * @return function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components
  2095.      * @since 6.1
  2096.      * @see #getXYSpXY2Function(TimeScales)
  2097.      */
  2098.     @DefaultDataContext
  2099.     public TimeVectorFunction getXYSpXY2Function() {
  2100.         return getXYSpXY2Function(DataContext.getDefault().getTimeScales());
  2101.     }

  2102.     /**
  2103.      * Get the function computing the Celestial Intermediate Pole and Celestial
  2104.      * Intermediate Origin components.
  2105.      * <p>
  2106.      * The returned function computes the two X, Y components of CIP and the S+XY/2
  2107.      * component of the non-rotating CIO.
  2108.      * </p>
  2109.      *
  2110.      * @param timeScales used to define the function.
  2111.      * @return function computing the Celestial Intermediate Pole and Celestial
  2112.      * Intermediate Origin components
  2113.      * @since 10.1
  2114.      */
  2115.     public abstract TimeVectorFunction getXYSpXY2Function(TimeScales timeScales);

  2116.     /** Get the function computing the raw Earth Orientation Angle.
  2117.      *
  2118.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2119.      *
  2120.      * <p>
  2121.      * The raw angle does not contain any correction. If for example dTU1 correction
  2122.      * due to tidal effect is desired, it must be added afterward by the caller.
  2123.      * The returned value contain the angle as the value and the angular rate as
  2124.      * the first derivative.
  2125.      * </p>
  2126.      * @param ut1 UT1 time scale
  2127.      * @return function computing the rawEarth Orientation Angle, in the non-rotating origin paradigm
  2128.      * @since 6.1
  2129.      * @see #getEarthOrientationAngleFunction(TimeScale, TimeScale)
  2130.      */
  2131.     @DefaultDataContext
  2132.     public TimeScalarFunction getEarthOrientationAngleFunction(final TimeScale ut1) {
  2133.         return getEarthOrientationAngleFunction(ut1,
  2134.                 DataContext.getDefault().getTimeScales().getTAI());
  2135.     }

  2136.     /** Get the function computing the raw Earth Orientation Angle.
  2137.      * <p>
  2138.      * The raw angle does not contain any correction. If for example dTU1 correction
  2139.      * due to tidal effect is desired, it must be added afterward by the caller.
  2140.      * The returned value contain the angle as the value and the angular rate as
  2141.      * the first derivative.
  2142.      * </p>
  2143.      * @param ut1 UT1 time scale
  2144.      * @param tai TAI time scale
  2145.      * @return function computing the rawEarth Orientation Angle, in the non-rotating origin paradigm
  2146.      * @since 10.1
  2147.      */
  2148.     public TimeScalarFunction getEarthOrientationAngleFunction(final TimeScale ut1,
  2149.                                                                final TimeScale tai) {
  2150.         return new StellarAngleCapitaine(ut1, tai);
  2151.     }

  2152.     /** Get the function computing the precession angles.
  2153.      * <p>
  2154.      * The function returned computes the three precession angles
  2155.      * ψ<sub>A</sub> (around Z axis), ω<sub>A</sub> (around X axis)
  2156.      * and χ<sub>A</sub> (around Z axis). The constant angle ε₀
  2157.      * for the fourth rotation (around X axis) can be retrieved by evaluating the
  2158.      * function returned by {@link #getMeanObliquityFunction()} at {@link
  2159.      * #getNutationReferenceEpoch() nutation reference epoch}.
  2160.      * </p>
  2161.      *
  2162.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2163.      *
  2164.      * @return function computing the precession angle
  2165.      * @since 6.1
  2166.      * @see #getPrecessionFunction(TimeScales)
  2167.      */
  2168.     @DefaultDataContext
  2169.     public TimeVectorFunction getPrecessionFunction()
  2170.     {
  2171.         return getPrecessionFunction(DataContext.getDefault().getTimeScales());
  2172.     }

  2173.     /** Get the function computing the precession angles.
  2174.      * <p>
  2175.      * The function returned computes the three precession angles
  2176.      * ψ<sub>A</sub> (around Z axis), ω<sub>A</sub> (around X axis)
  2177.      * and χ<sub>A</sub> (around Z axis). The constant angle ε₀
  2178.      * for the fourth rotation (around X axis) can be retrieved by evaluating the
  2179.      * function returned by {@link #getMeanObliquityFunction()} at {@link
  2180.      * #getNutationReferenceEpoch() nutation reference epoch}.
  2181.      * </p>
  2182.      * @return function computing the precession angle
  2183.      * @since 10.1
  2184.      * @param timeScales used to define the function.
  2185.      */
  2186.     public abstract TimeVectorFunction getPrecessionFunction(TimeScales timeScales);

  2187.     /** Get the function computing the nutation angles.
  2188.      *
  2189.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2190.      *
  2191.      * <p>
  2192.      * The function returned computes the two classical angles ΔΨ and Δε,
  2193.      * and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
  2194.      * resolution C7 (the correction is forced to 0 before this date)
  2195.      * </p>
  2196.      * @return function computing the nutation in longitude ΔΨ and Δε
  2197.      * and the correction of equation of equinoxes
  2198.      * @since 6.1
  2199.      */
  2200.     @DefaultDataContext
  2201.     public TimeVectorFunction getNutationFunction() {
  2202.         return getNutationFunction(DataContext.getDefault().getTimeScales());
  2203.     }

  2204.     /** Get the function computing the nutation angles.
  2205.      * <p>
  2206.      * The function returned computes the two classical angles ΔΨ and Δε,
  2207.      * and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
  2208.      * resolution C7 (the correction is forced to 0 before this date)
  2209.      * </p>
  2210.      * @return function computing the nutation in longitude ΔΨ and Δε
  2211.      * and the correction of equation of equinoxes
  2212.      * @param timeScales used in the computation including TAI and TT.
  2213.      * @since 10.1
  2214.      */
  2215.     public abstract TimeVectorFunction getNutationFunction(TimeScales timeScales);

  2216.     /** Get the function computing Greenwich mean sidereal time, in radians.
  2217.      *
  2218.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2219.      *
  2220.      * @param ut1 UT1 time scale
  2221.      * @return function computing Greenwich mean sidereal time
  2222.      * @since 6.1
  2223.      * @see #getGMSTFunction(TimeScale, TimeScales)
  2224.      */
  2225.     @DefaultDataContext
  2226.     public TimeScalarFunction getGMSTFunction(final TimeScale ut1) {
  2227.         return getGMSTFunction(ut1, DataContext.getDefault().getTimeScales());
  2228.     }

  2229.     /**
  2230.      * Get the function computing Greenwich mean sidereal time, in radians.
  2231.      *
  2232.      * @param ut1 UT1 time scale
  2233.      * @param timeScales other time scales used in the computation including TAI and TT.
  2234.      * @return function computing Greenwich mean sidereal time
  2235.      * @since 10.1
  2236.      */
  2237.     public abstract TimeScalarFunction getGMSTFunction(TimeScale ut1,
  2238.                                                        TimeScales timeScales);

  2239.     /** Get the function computing Greenwich mean sidereal time rate, in radians per second.
  2240.      *
  2241.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2242.      *
  2243.      * @param ut1 UT1 time scale
  2244.      * @return function computing Greenwich mean sidereal time rate
  2245.      * @since 9.0
  2246.      * @see #getGMSTRateFunction(TimeScale, TimeScales)
  2247.      */
  2248.     @DefaultDataContext
  2249.     public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1) {
  2250.         return getGMSTRateFunction(ut1,
  2251.                 DataContext.getDefault().getTimeScales());
  2252.     }

  2253.     /**
  2254.      * Get the function computing Greenwich mean sidereal time rate, in radians per
  2255.      * second.
  2256.      *
  2257.      * @param ut1 UT1 time scale
  2258.      * @param timeScales other time scales used in the computation including TAI and TT.
  2259.      * @return function computing Greenwich mean sidereal time rate
  2260.      * @since 10.1
  2261.      */
  2262.     public abstract TimeScalarFunction getGMSTRateFunction(TimeScale ut1,
  2263.                                                            TimeScales timeScales);

  2264.     /**
  2265.      * Get the function computing Greenwich apparent sidereal time, in radians.
  2266.      *
  2267.      * <p>This method uses the {@link DataContext#getDefault() default data context} if
  2268.      * {@code eopHistory == null}.
  2269.      *
  2270.      * @param ut1        UT1 time scale
  2271.      * @param eopHistory EOP history. If {@code null} then no nutation correction is
  2272.      *                   applied for EOP.
  2273.      * @return function computing Greenwich apparent sidereal time
  2274.      * @since 6.1
  2275.      * @see #getGASTFunction(TimeScale, EOPHistory, TimeScales)
  2276.      */
  2277.     @DefaultDataContext
  2278.     public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  2279.                                               final EOPHistory eopHistory) {
  2280.         final TimeScales timeScales = eopHistory != null ?
  2281.                 eopHistory.getTimeScales() :
  2282.                 DataContext.getDefault().getTimeScales();
  2283.         return getGASTFunction(ut1, eopHistory, timeScales);
  2284.     }

  2285.     /**
  2286.      * Get the function computing Greenwich apparent sidereal time, in radians.
  2287.      *
  2288.      * @param ut1        UT1 time scale
  2289.      * @param eopHistory EOP history. If {@code null} then no nutation correction is
  2290.      *                   applied for EOP.
  2291.      * @param timeScales        TAI time scale.
  2292.      * @return function computing Greenwich apparent sidereal time
  2293.      * @since 10.1
  2294.      */
  2295.     public abstract TimeScalarFunction getGASTFunction(TimeScale ut1,
  2296.                                                        EOPHistory eopHistory,
  2297.                                                        TimeScales timeScales);

  2298.     /** Get the function computing tidal corrections for Earth Orientation Parameters.
  2299.      *
  2300.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2301.      *
  2302.      * @return function computing tidal corrections for Earth Orientation Parameters,
  2303.      * for xp, yp, ut1 and lod respectively
  2304.      * @since 6.1
  2305.      * @see #getEOPTidalCorrection(TimeScales)
  2306.      */
  2307.     @DefaultDataContext
  2308.     public TimeVectorFunction getEOPTidalCorrection() {
  2309.         return getEOPTidalCorrection(DataContext.getDefault().getTimeScales());
  2310.     }

  2311.     /**
  2312.      * Get the function computing tidal corrections for Earth Orientation Parameters.
  2313.      *
  2314.      * @param timeScales used in the computation. The TT and TAI scales are used.
  2315.      * @return function computing tidal corrections for Earth Orientation Parameters, for
  2316.      * xp, yp, ut1 and lod respectively
  2317.      * @since 10.1
  2318.      */
  2319.     public abstract TimeVectorFunction getEOPTidalCorrection(TimeScales timeScales);

  2320.     /** Get the Love numbers.
  2321.      * @return Love numbers
  2322.           * @since 6.1
  2323.      */
  2324.     public abstract LoveNumbers getLoveNumbers();

  2325.     /** Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2326.      *
  2327.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2328.      *
  2329.      * @param ut1 UT1 time scale
  2330.      * @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2331.      * @since 6.1
  2332.      * @see #getTideFrequencyDependenceFunction(TimeScale, TimeScales)
  2333.      */
  2334.     @DefaultDataContext
  2335.     public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1) {
  2336.         return getTideFrequencyDependenceFunction(ut1,
  2337.                 DataContext.getDefault().getTimeScales());
  2338.     }

  2339.     /**
  2340.      * Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂,
  2341.      * ΔS₂₂).
  2342.      *
  2343.      * @param ut1 UT1 time scale
  2344.      * @param timeScales other time scales used in the computation including TAI and TT.
  2345.      * @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂,
  2346.      * ΔS₂₂).
  2347.      * @since 10.1
  2348.      */
  2349.     public abstract TimeVectorFunction getTideFrequencyDependenceFunction(
  2350.             TimeScale ut1,
  2351.             TimeScales timeScales);

  2352.     /** Get the permanent tide to be <em>removed</em> from ΔC₂₀ when zero-tide potentials are used.
  2353.      * @return permanent tide to remove
  2354.      */
  2355.     public abstract double getPermanentTide();

  2356.     /** Get the function computing solid pole tide (ΔC₂₁, ΔS₂₁).
  2357.      * @param eopHistory EOP history
  2358.      * @return model for solid pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2359.           * @since 6.1
  2360.      */
  2361.     public abstract TimeVectorFunction getSolidPoleTide(EOPHistory eopHistory);

  2362.     /** Get the function computing ocean pole tide (ΔC₂₁, ΔS₂₁).
  2363.      * @param eopHistory EOP history
  2364.      * @return model for ocean pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2365.           * @since 6.1
  2366.      */
  2367.     public abstract TimeVectorFunction getOceanPoleTide(EOPHistory eopHistory);

  2368.     /** Get the nominal values of the displacement numbers.
  2369.      * @return an array containing h⁽⁰⁾, h⁽²⁾, h₃, hI diurnal, hI semi-diurnal,
  2370.      * l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾, l₃, lI diurnal, lI semi-diurnal,
  2371.      * H₀ permanent deformation amplitude
  2372.      * @since 9.1
  2373.      */
  2374.     public abstract double[] getNominalTidalDisplacement();

  2375.     /** Get the correction function for tidal displacement for diurnal tides.
  2376.      * <ul>
  2377.      *  <li>f[0]: radial correction, longitude cosine part</li>
  2378.      *  <li>f[1]: radial correction, longitude sine part</li>
  2379.      *  <li>f[2]: North correction, longitude cosine part</li>
  2380.      *  <li>f[3]: North correction, longitude sine part</li>
  2381.      *  <li>f[4]: East correction, longitude cosine part</li>
  2382.      *  <li>f[5]: East correction, longitude sine part</li>
  2383.      * </ul>
  2384.      * @return correction function for tidal displacement
  2385.      * @since 9.1
  2386.      */
  2387.     public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal();

  2388.     /** Get the correction function for tidal displacement for diurnal tides.
  2389.      * <ul>
  2390.      *  <li>f[0]: radial correction, longitude cosine part</li>
  2391.      *  <li>f[1]: radial correction, longitude sine part</li>
  2392.      *  <li>f[2]: North correction, longitude cosine part</li>
  2393.      *  <li>f[3]: North correction, longitude sine part</li>
  2394.      *  <li>f[4]: East correction, longitude cosine part</li>
  2395.      *  <li>f[5]: East correction, longitude sine part</li>
  2396.      * </ul>
  2397.      * @param tableName name for the diurnal tides table
  2398.      * @param cols total number of columns of the diurnal tides table
  2399.      * @param rIp column holding ∆Rf(ip) in the diurnal tides table, counting from 1
  2400.      * @param rOp column holding ∆Rf(op) in the diurnal tides table, counting from 1
  2401.      * @param tIp column holding ∆Tf(ip) in the diurnal tides table, counting from 1
  2402.      * @param tOp column holding ∆Tf(op) in the diurnal tides table, counting from 1
  2403.      * @return correction function for tidal displacement for diurnal tides
  2404.           * @since 9.1
  2405.      */
  2406.     protected static CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal(final String tableName, final int cols,
  2407.                                                                                    final int rIp, final int rOp,
  2408.                                                                                    final int tIp, final int tOp) {

  2409.         try (InputStream in1 = getStream(tableName);
  2410.              InputStream in2 = getStream(tableName);
  2411.              InputStream in3 = getStream(tableName);
  2412.              InputStream in4 = getStream(tableName);
  2413.              InputStream in5 = getStream(tableName);
  2414.              InputStream in6 = getStream(tableName)) {
  2415.             // radial component, missing the sin 2φ factor; this corresponds to:
  2416.             //  - equation 15a in IERS conventions 1996, chapter 7
  2417.             //  - equation 16a in IERS conventions 2003, chapter 7
  2418.             //  - equation 7.12a in IERS conventions 2010, chapter 7
  2419.             final PoissonSeries drCos = new PoissonSeriesParser(cols).
  2420.                                         withOptionalColumn(1).
  2421.                                         withDoodson(4, 3).
  2422.                                         withFirstDelaunay(10).
  2423.                                         withSinCos(0, rIp, +1.0e-3, rOp, +1.0e-3).
  2424.                                         parse(in1, tableName);
  2425.             final PoissonSeries drSin = new PoissonSeriesParser(cols).
  2426.                                         withOptionalColumn(1).
  2427.                                         withDoodson(4, 3).
  2428.                                         withFirstDelaunay(10).
  2429.                                         withSinCos(0, rOp, -1.0e-3, rIp, +1.0e-3).
  2430.                                         parse(in2, tableName);

  2431.             // North component, missing the cos 2φ factor; this corresponds to:
  2432.             //  - equation 15b in IERS conventions 1996, chapter 7
  2433.             //  - equation 16b in IERS conventions 2003, chapter 7
  2434.             //  - equation 7.12b in IERS conventions 2010, chapter 7
  2435.             final PoissonSeries dnCos = new PoissonSeriesParser(cols).
  2436.                                         withOptionalColumn(1).
  2437.                                         withDoodson(4, 3).
  2438.                                         withFirstDelaunay(10).
  2439.                                         withSinCos(0, tIp, +1.0e-3, tOp, +1.0e-3).
  2440.                                         parse(in3, tableName);
  2441.             final PoissonSeries dnSin = new PoissonSeriesParser(cols).
  2442.                                         withOptionalColumn(1).
  2443.                                         withDoodson(4, 3).
  2444.                                         withFirstDelaunay(10).
  2445.                                         withSinCos(0, tOp, -1.0e-3, tIp, +1.0e-3).
  2446.                                         parse(in4, tableName);

  2447.             // East component, missing the sin φ factor; this corresponds to:
  2448.             //  - equation 15b in IERS conventions 1996, chapter 7
  2449.             //  - equation 16b in IERS conventions 2003, chapter 7
  2450.             //  - equation 7.12b in IERS conventions 2010, chapter 7
  2451.             final PoissonSeries deCos = new PoissonSeriesParser(cols).
  2452.                                         withOptionalColumn(1).
  2453.                                         withDoodson(4, 3).
  2454.                                         withFirstDelaunay(10).
  2455.                                         withSinCos(0, tOp, -1.0e-3, tIp, +1.0e-3).
  2456.                                         parse(in5, tableName);
  2457.             final PoissonSeries deSin = new PoissonSeriesParser(cols).
  2458.                                         withOptionalColumn(1).
  2459.                                         withDoodson(4, 3).
  2460.                                         withFirstDelaunay(10).
  2461.                                         withSinCos(0, tIp, -1.0e-3, tOp, -1.0e-3).
  2462.                                         parse(in6, tableName);

  2463.             return PoissonSeries.compile(drCos, drSin,
  2464.                                          dnCos, dnSin,
  2465.                                          deCos, deSin);
  2466.         } catch (IOException e) {
  2467.             throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2468.         }

  2469.     }

  2470.     /** Get the correction function for tidal displacement for zonal tides.
  2471.      * <ul>
  2472.      *  <li>f[0]: radial correction</li>
  2473.      *  <li>f[1]: North correction</li>
  2474.      * </ul>
  2475.      * @return correction function for tidal displacement
  2476.      * @since 9.1
  2477.      */
  2478.     public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionZonal();

  2479.     /** Get the correction function for tidal displacement for zonal tides.
  2480.      * <ul>
  2481.      *  <li>f[0]: radial correction</li>
  2482.      *  <li>f[1]: North correction</li>
  2483.      * </ul>
  2484.      * @param tableName name for the zonal tides table
  2485.      * @param cols total number of columns of the table
  2486.      * @param rIp column holding ∆Rf(ip) in the table, counting from 1
  2487.      * @param rOp column holding ∆Rf(op) in the table, counting from 1
  2488.      * @param tIp column holding ∆Tf(ip) in the table, counting from 1
  2489.      * @param tOp column holding ∆Tf(op) in the table, counting from 1
  2490.      * @return correction function for tidal displacement for zonal tides
  2491.           * @since 9.1
  2492.      */
  2493.     protected static CompiledSeries getTidalDisplacementFrequencyCorrectionZonal(final String tableName, final int cols,
  2494.                                                                                  final int rIp, final int rOp,
  2495.                                                                                  final int tIp, final int tOp) {

  2496.         try (InputStream in1 = getStream(tableName);
  2497.              InputStream in2 = getStream(tableName)) {
  2498.             // radial component, missing the 3⁄2 sin² φ - 1⁄2 factor; this corresponds to:
  2499.             //  - equation 16a in IERS conventions 1996, chapter 7
  2500.             //  - equation 17a in IERS conventions 2003, chapter 7
  2501.             //  - equation 7.13a in IERS conventions 2010, chapter 7
  2502.             final PoissonSeries dr = new PoissonSeriesParser(cols).
  2503.                                      withOptionalColumn(1).
  2504.                                      withDoodson(4, 3).
  2505.                                      withFirstDelaunay(10).
  2506.                                      withSinCos(0, rOp, +1.0e-3, rIp, +1.0e-3).
  2507.                                      parse(in1, tableName);

  2508.             // North component, missing the sin 2φ factor; this corresponds to:
  2509.             //  - equation 16b in IERS conventions 1996, chapter 7
  2510.             //  - equation 17b in IERS conventions 2003, chapter 7
  2511.             //  - equation 7.13b in IERS conventions 2010, chapter 7
  2512.             final PoissonSeries dn = new PoissonSeriesParser(cols).
  2513.                                      withOptionalColumn(1).
  2514.                                      withDoodson(4, 3).
  2515.                                      withFirstDelaunay(10).
  2516.                                      withSinCos(0, tOp, +1.0e-3, tIp, +1.0e-3).
  2517.                                      parse(in2, tableName);

  2518.             return PoissonSeries.compile(dr, dn);
  2519.         } catch (IOException e) {
  2520.             throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2521.         }

  2522.     }

  2523.     /** Interface for functions converting nutation corrections between
  2524.      * δΔψ/δΔε to δX/δY.
  2525.      * <ul>
  2526.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  2527.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  2528.      * </ul>
  2529.      * @since 6.1
  2530.      */
  2531.     public interface NutationCorrectionConverter {

  2532.         /** Convert nutation corrections.
  2533.          * @param date current date
  2534.          * @param ddPsi δΔψ part of the nutation correction
  2535.          * @param ddEpsilon δΔε part of the nutation correction
  2536.          * @return array containing δX and δY
  2537.          */
  2538.         double[] toNonRotating(AbsoluteDate date, double ddPsi, double ddEpsilon);

  2539.         /** Convert nutation corrections.
  2540.          * @param date current date
  2541.          * @param dX δX part of the nutation correction
  2542.          * @param dY δY part of the nutation correction
  2543.          * @return array containing δΔψ and δΔε
  2544.          */
  2545.         double[] toEquinox(AbsoluteDate date, double dX, double dY);

  2546.     }

  2547.     /** Create a function converting nutation corrections between
  2548.      * δX/δY and δΔψ/δΔε.
  2549.      * <ul>
  2550.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  2551.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  2552.      * </ul>
  2553.      *
  2554.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2555.      *
  2556.      * @return a new converter
  2557.      * @since 6.1
  2558.      * @see #getNutationCorrectionConverter(TimeScales)
  2559.      */
  2560.     @DefaultDataContext
  2561.     public NutationCorrectionConverter getNutationCorrectionConverter() {
  2562.         return getNutationCorrectionConverter(DataContext.getDefault().getTimeScales());
  2563.     }

  2564.     /** Create a function converting nutation corrections between
  2565.      * δX/δY and δΔψ/δΔε.
  2566.      * <ul>
  2567.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  2568.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  2569.      * </ul>
  2570.      * @return a new converter
  2571.      * @since 10.1
  2572.      * @param timeScales used to define the conversion.
  2573.      */
  2574.     public NutationCorrectionConverter getNutationCorrectionConverter(
  2575.             final TimeScales timeScales) {

  2576.         // get models parameters
  2577.         final TimeVectorFunction precessionFunction = getPrecessionFunction(timeScales);
  2578.         final TimeScalarFunction epsilonAFunction = getMeanObliquityFunction(timeScales);
  2579.         final AbsoluteDate date0 = getNutationReferenceEpoch(timeScales);
  2580.         final double cosE0 = FastMath.cos(epsilonAFunction.value(date0));

  2581.         return new NutationCorrectionConverter() {

  2582.             /** {@inheritDoc} */
  2583.             @Override
  2584.             public double[] toNonRotating(final AbsoluteDate date,
  2585.                                           final double ddPsi, final double ddEpsilon) {
  2586.                 // compute precession angles psiA, omegaA and chiA
  2587.                 final double[] angles = precessionFunction.value(date);

  2588.                 // conversion coefficients
  2589.                 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
  2590.                 final double c     = angles[0] * cosE0 - angles[2];

  2591.                 // convert nutation corrections (equation 23/IERS-2003 or 5.25/IERS-2010)
  2592.                 return new double[] {
  2593.                     sinEA * ddPsi + c * ddEpsilon,
  2594.                     ddEpsilon - c * sinEA * ddPsi
  2595.                 };

  2596.             }

  2597.             /** {@inheritDoc} */
  2598.             @Override
  2599.             public double[] toEquinox(final AbsoluteDate date,
  2600.                                       final double dX, final double dY) {
  2601.                 // compute precession angles psiA, omegaA and chiA
  2602.                 final double[] angles   = precessionFunction.value(date);

  2603.                 // conversion coefficients
  2604.                 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
  2605.                 final double c     = angles[0] * cosE0 - angles[2];
  2606.                 final double opC2  = 1 + c * c;

  2607.                 // convert nutation corrections (inverse of equation 23/IERS-2003 or 5.25/IERS-2010)
  2608.                 return new double[] {
  2609.                     (dX - c * dY) / (sinEA * opC2),
  2610.                     (dY + c * dX) / opC2
  2611.                 };
  2612.             }

  2613.         };

  2614.     }

  2615.     /** Load the Love numbers.
  2616.      * @param nameLove name of the Love number resource
  2617.      * @return Love numbers
  2618.      */
  2619.     protected LoveNumbers loadLoveNumbers(final String nameLove) {
  2620.         try {

  2621.             // allocate the three triangular arrays for real, imaginary and time-dependent numbers
  2622.             final double[][] real      = new double[4][];
  2623.             final double[][] imaginary = new double[4][];
  2624.             final double[][] plus      = new double[4][];
  2625.             for (int i = 0; i < real.length; ++i) {
  2626.                 real[i]      = new double[i + 1];
  2627.                 imaginary[i] = new double[i + 1];
  2628.                 plus[i]      = new double[i + 1];
  2629.             }

  2630.             try (InputStream stream = IERSConventions.class.getResourceAsStream(nameLove)) {

  2631.                 if (stream == null) {
  2632.                     // this should never happen with files embedded within Orekit
  2633.                     throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, nameLove);
  2634.                 }

  2635.                 // setup the reader
  2636.                 try (BufferedReader reader = new BufferedReader(new InputStreamReader(stream, StandardCharsets.UTF_8))) {

  2637.                     String line = reader.readLine();
  2638.                     int lineNumber = 1;

  2639.                     // look for the Love numbers
  2640.                     while (line != null) {

  2641.                         line = line.trim();
  2642.                         if (!(line.isEmpty() || line.startsWith("#"))) {
  2643.                             final String[] fields = line.split("\\p{Space}+");
  2644.                             if (fields.length != 5) {
  2645.                                 // this should never happen with files embedded within Orekit
  2646.                                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  2647.                                                           lineNumber, nameLove, line);
  2648.                             }
  2649.                             final int n = Integer.parseInt(fields[0]);
  2650.                             final int m = Integer.parseInt(fields[1]);
  2651.                             if ((n < 2) || (n > 3) || (m < 0) || (m > n)) {
  2652.                                 // this should never happen with files embedded within Orekit
  2653.                                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  2654.                                                           lineNumber, nameLove, line);

  2655.                             }
  2656.                             real[n][m]      = Double.parseDouble(fields[2]);
  2657.                             imaginary[n][m] = Double.parseDouble(fields[3]);
  2658.                             plus[n][m]      = Double.parseDouble(fields[4]);
  2659.                         }

  2660.                         // next line
  2661.                         lineNumber++;
  2662.                         line = reader.readLine();

  2663.                     }
  2664.                 }
  2665.             }

  2666.             return new LoveNumbers(real, imaginary, plus);

  2667.         } catch (IOException ioe) {
  2668.             // this should never happen with files embedded within Orekit
  2669.             throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, nameLove);
  2670.         }
  2671.     }

  2672.     /** Get a data stream.
  2673.      * @param name file name of the resource stream
  2674.      * @return stream
  2675.      */
  2676.     private static InputStream getStream(final String name) {
  2677.         return IERSConventions.class.getResourceAsStream(name);
  2678.     }

  2679.     /** Correction to equation of equinoxes.
  2680.      * <p>IAU 1994 resolution C7 added two terms to the equation of equinoxes
  2681.      * taking effect since 1997-02-27 for continuity
  2682.      * </p>
  2683.      */
  2684.     private static class IAU1994ResolutionC7 {

  2685.         /** First Moon correction term for the Equation of the Equinoxes. */
  2686.         private static final double EQE1 =     0.00264  * Constants.ARC_SECONDS_TO_RADIANS;

  2687.         /** Second Moon correction term for the Equation of the Equinoxes. */
  2688.         private static final double EQE2 =     0.000063 * Constants.ARC_SECONDS_TO_RADIANS;


  2689.         /** Evaluate the correction.
  2690.          * @param arguments Delaunay for nutation
  2691.          * @param tai TAI time scale.
  2692.          * @return correction value (0 before 1997-02-27)
  2693.          */
  2694.         public static double value(final DelaunayArguments arguments,
  2695.                                    final TimeScale tai) {
  2696.             /* Start date for applying Moon corrections to the equation of the equinoxes.
  2697.              * This date corresponds to 1997-02-27T00:00:00 UTC, hence the 30s offset from TAI.
  2698.              */
  2699.             final AbsoluteDate modelStart = new AbsoluteDate(1997, 2, 27, 0, 0, 30, tai);
  2700.             if (arguments.getDate().compareTo(modelStart) >= 0) {

  2701.                 // IAU 1994 resolution C7 added two terms to the equation of equinoxes
  2702.                 // taking effect since 1997-02-27 for continuity

  2703.                 // Mean longitude of the ascending node of the Moon
  2704.                 final double om = arguments.getOmega();

  2705.                 // add the two correction terms
  2706.                 return EQE1 * FastMath.sin(om) + EQE2 * FastMath.sin(om + om);

  2707.             } else {
  2708.                 return 0.0;
  2709.             }
  2710.         }

  2711.         /** Evaluate the correction.
  2712.          * @param arguments Delaunay for nutation
  2713.          * @param tai TAI time scale.
  2714.          * @param <T> type of the field elements
  2715.          * @return correction value (0 before 1997-02-27)
  2716.          */
  2717.         public static <T extends RealFieldElement<T>> T value(
  2718.                 final FieldDelaunayArguments<T> arguments,
  2719.                 final TimeScale tai) {
  2720.             /* Start date for applying Moon corrections to the equation of the equinoxes.
  2721.              * This date corresponds to 1997-02-27T00:00:00 UTC, hence the 30s offset from TAI.
  2722.              */
  2723.             final AbsoluteDate modelStart = new AbsoluteDate(1997, 2, 27, 0, 0, 30, tai);
  2724.             if (arguments.getDate().toAbsoluteDate().compareTo(modelStart) >= 0) {

  2725.                 // IAU 1994 resolution C7 added two terms to the equation of equinoxes
  2726.                 // taking effect since 1997-02-27 for continuity

  2727.                 // Mean longitude of the ascending node of the Moon
  2728.                 final T om = arguments.getOmega();

  2729.                 // add the two correction terms
  2730.                 return om.sin().multiply(EQE1).add(om.add(om).sin().multiply(EQE2));

  2731.             } else {
  2732.                 return arguments.getDate().getField().getZero();
  2733.             }
  2734.         }

  2735.     };

  2736.     /** Stellar angle model.
  2737.      * <p>
  2738.      * The stellar angle computed here has been defined in the paper "A non-rotating origin on the
  2739.      * instantaneous equator: Definition, properties and use", N. Capitaine, Guinot B., and Souchay J.,
  2740.      * Celestial Mechanics, Volume 39, Issue 3, pp 283-307. It has been proposed as a conventional
  2741.      * conventional relationship between UT1 and Earth rotation in the paper "Definition of the Celestial
  2742.      * Ephemeris origin and of UT1 in the International Celestial Reference Frame", Capitaine, N.,
  2743.      * Guinot, B., and McCarthy, D. D., 2000, Astronomy and Astrophysics, 355(1), pp. 398–405.
  2744.      * </p>
  2745.      * <p>
  2746.      * It is presented simply as stellar angle in IERS conventions 1996 but as since been adopted as
  2747.      * the conventional relationship defining UT1 from ICRF and is called Earth Rotation Angle in
  2748.      * IERS conventions 2003 and 2010.
  2749.      * </p>
  2750.      */
  2751.     private static class StellarAngleCapitaine implements TimeScalarFunction {

  2752.         /** Constant term of Capitaine's Earth Rotation Angle model. */
  2753.         private static final double ERA_0   = MathUtils.TWO_PI * 0.7790572732640;

  2754.         /** Rate term of Capitaine's Earth Rotation Angle model.
  2755.          * (radians per day, main part) */
  2756.         private static final double ERA_1A  = MathUtils.TWO_PI / Constants.JULIAN_DAY;

  2757.         /** Rate term of Capitaine's Earth Rotation Angle model.
  2758.          * (radians per day, fractional part) */
  2759.         private static final double ERA_1B  = ERA_1A * 0.00273781191135448;

  2760.         /** UT1 time scale. */
  2761.         private final TimeScale ut1;

  2762.         /** Reference date of Capitaine's Earth Rotation Angle model. */
  2763.         private final AbsoluteDate referenceDate;

  2764.         /** Simple constructor.
  2765.          * @param ut1 UT1 time scale
  2766.          * @param tai TAI time scale
  2767.          */
  2768.         StellarAngleCapitaine(final TimeScale ut1, final TimeScale tai) {
  2769.             this.ut1 = ut1;
  2770.             referenceDate = new AbsoluteDate(
  2771.                     DateComponents.J2000_EPOCH,
  2772.                     TimeComponents.H12,
  2773.                     tai);
  2774.         }

  2775.         /** Get the rotation rate.
  2776.          * @return rotation rate
  2777.          */
  2778.         public double getRate() {
  2779.             return ERA_1A + ERA_1B;
  2780.         }

  2781.         /** {@inheritDoc} */
  2782.         @Override
  2783.         public double value(final AbsoluteDate date) {

  2784.             // split the date offset as a full number of days plus a smaller part
  2785.             final int secondsInDay = 86400;
  2786.             final double dt  = date.durationFrom(referenceDate);
  2787.             final long days  = ((long) dt) / secondsInDay;
  2788.             final double dtA = secondsInDay * days;
  2789.             final double dtB = (dt - dtA) + ut1.offsetFromTAI(date);

  2790.             return ERA_0 + ERA_1A * dtB + ERA_1B * (dtA + dtB);

  2791.         }

  2792.         /** {@inheritDoc} */
  2793.         @Override
  2794.         public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  2795.             // split the date offset as a full number of days plus a smaller part
  2796.             final int secondsInDay = 86400;
  2797.             final T dt  = date.durationFrom(referenceDate);
  2798.             final long days  = ((long) dt.getReal()) / secondsInDay;
  2799.             final double dtA = secondsInDay * days;
  2800.             final T dtB = dt.subtract(dtA).add(ut1.offsetFromTAI(date.toAbsoluteDate()));

  2801.             return dtB.add(dtA).multiply(ERA_1B).add(dtB.multiply(ERA_1A)).add(ERA_0);

  2802.         }

  2803.     }

  2804.     /** Mean pole. */
  2805.     private static class MeanPole implements TimeStamped, Serializable {

  2806.         /** Serializable UID. */
  2807.         private static final long serialVersionUID = 20131028l;

  2808.         /** Date. */
  2809.         private final AbsoluteDate date;

  2810.         /** X coordinate. */
  2811.         private double x;

  2812.         /** Y coordinate. */
  2813.         private double y;

  2814.         /** Simple constructor.
  2815.          * @param date date
  2816.          * @param x x coordinate
  2817.          * @param y y coordinate
  2818.          */
  2819.         MeanPole(final AbsoluteDate date, final double x, final double y) {
  2820.             this.date = date;
  2821.             this.x    = x;
  2822.             this.y    = y;
  2823.         }

  2824.         /** {@inheritDoc} */
  2825.         @Override
  2826.         public AbsoluteDate getDate() {
  2827.             return date;
  2828.         }

  2829.         /** Get x coordinate.
  2830.          * @return x coordinate
  2831.          */
  2832.         public double getX() {
  2833.             return x;
  2834.         }

  2835.         /** Get y coordinate.
  2836.          * @return y coordinate
  2837.          */
  2838.         public double getY() {
  2839.             return y;
  2840.         }

  2841.     }

  2842.     /** Local class for precession function. */
  2843.     private class PrecessionFunction implements TimeVectorFunction {

  2844.         /** Polynomial nutation for psiA. */
  2845.         private final PolynomialNutation psiA;

  2846.         /** Polynomial nutation for omegaA. */
  2847.         private final PolynomialNutation omegaA;

  2848.         /** Polynomial nutation for chiA. */
  2849.         private final PolynomialNutation chiA;

  2850.         /** Time scales to use. */
  2851.         private final TimeScales timeScales;

  2852.         /** Simple constructor.
  2853.          * @param psiA polynomial nutation for psiA
  2854.          * @param omegaA polynomial nutation for omegaA
  2855.          * @param chiA polynomial nutation for chiA
  2856.          * @param timeScales used in the computation.
  2857.          */
  2858.         PrecessionFunction(final PolynomialNutation psiA,
  2859.                            final PolynomialNutation omegaA,
  2860.                            final PolynomialNutation chiA,
  2861.                            final TimeScales timeScales) {
  2862.             this.psiA   = psiA;
  2863.             this.omegaA = omegaA;
  2864.             this.chiA   = chiA;
  2865.             this.timeScales = timeScales;
  2866.         }


  2867.         /** {@inheritDoc} */
  2868.         @Override
  2869.         public double[] value(final AbsoluteDate date) {
  2870.             final double tc = evaluateTC(date, timeScales);
  2871.             return new double[] {
  2872.                 psiA.value(tc), omegaA.value(tc), chiA.value(tc)
  2873.             };
  2874.         }

  2875.         /** {@inheritDoc} */
  2876.         @Override
  2877.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  2878.             final T[] a = MathArrays.buildArray(date.getField(), 3);
  2879.             final T tc = evaluateTC(date, timeScales);
  2880.             a[0] = psiA.value(tc);
  2881.             a[1] = omegaA.value(tc);
  2882.             a[2] = chiA.value(tc);
  2883.             return a;
  2884.         }

  2885.     }

  2886.     /** Local class for tides frequency function. */
  2887.     private static class TideFrequencyDependenceFunction implements TimeVectorFunction {

  2888.         /** Nutation arguments. */
  2889.         private final FundamentalNutationArguments arguments;

  2890.         /** Correction series. */
  2891.         private final PoissonSeries.CompiledSeries kSeries;

  2892.         /** Simple constructor.
  2893.          * @param arguments nutation arguments
  2894.          * @param c20Series correction series for the C20 term
  2895.          * @param c21Series correction series for the C21 term
  2896.          * @param s21Series correction series for the S21 term
  2897.          * @param c22Series correction series for the C22 term
  2898.          * @param s22Series correction series for the S22 term
  2899.          */
  2900.         TideFrequencyDependenceFunction(final FundamentalNutationArguments arguments,
  2901.                                         final PoissonSeries c20Series,
  2902.                                         final PoissonSeries c21Series,
  2903.                                         final PoissonSeries s21Series,
  2904.                                         final PoissonSeries c22Series,
  2905.                                         final PoissonSeries s22Series) {
  2906.             this.arguments = arguments;
  2907.             this.kSeries   = PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);
  2908.         }


  2909.         /** {@inheritDoc} */
  2910.         @Override
  2911.         public double[] value(final AbsoluteDate date) {
  2912.             return kSeries.value(arguments.evaluateAll(date));
  2913.         }

  2914.         /** {@inheritDoc} */
  2915.         @Override
  2916.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  2917.             return kSeries.value(arguments.evaluateAll(date));
  2918.         }

  2919.     }

  2920.     /** Local class for EOP tidal corrections. */
  2921.     private static class EOPTidalCorrection implements TimeVectorFunction {

  2922.         /** Nutation arguments. */
  2923.         private final FundamentalNutationArguments arguments;

  2924.         /** Correction series. */
  2925.         private final PoissonSeries.CompiledSeries correctionSeries;

  2926.         /** Simple constructor.
  2927.          * @param arguments nutation arguments
  2928.          * @param xSeries correction series for the x coordinate
  2929.          * @param ySeries correction series for the y coordinate
  2930.          * @param ut1Series correction series for the UT1
  2931.          */
  2932.         EOPTidalCorrection(final FundamentalNutationArguments arguments,
  2933.                            final PoissonSeries xSeries,
  2934.                            final PoissonSeries ySeries,
  2935.                            final PoissonSeries ut1Series) {
  2936.             this.arguments        = arguments;
  2937.             this.correctionSeries = PoissonSeries.compile(xSeries, ySeries, ut1Series);
  2938.         }

  2939.         /** {@inheritDoc} */
  2940.         @Override
  2941.         public double[] value(final AbsoluteDate date) {
  2942.             final BodiesElements elements = arguments.evaluateAll(date);
  2943.             final double[] correction    = correctionSeries.value(elements);
  2944.             final double[] correctionDot = correctionSeries.derivative(elements);
  2945.             return new double[] {
  2946.                 correction[0],
  2947.                 correction[1],
  2948.                 correction[2],
  2949.                 correctionDot[2] * (-Constants.JULIAN_DAY)
  2950.             };
  2951.         }

  2952.         /** {@inheritDoc} */
  2953.         @Override
  2954.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

  2955.             final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  2956.             final T[] correction    = correctionSeries.value(elements);
  2957.             final T[] correctionDot = correctionSeries.derivative(elements);
  2958.             final T[] a = MathArrays.buildArray(date.getField(), 4);
  2959.             a[0] = correction[0];
  2960.             a[1] = correction[1];
  2961.             a[2] = correction[2];
  2962.             a[3] = correctionDot[2].multiply(-Constants.JULIAN_DAY);
  2963.             return a;
  2964.         }

  2965.     }

  2966. }