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;

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

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


  66. /** Supported IERS conventions.
  67.  * @since 6.0
  68.  * @author Luc Maisonobe
  69.  */
  70. public enum IERSConventions {

  72.     /** Constant for IERS 1996 conventions. */
  73.     IERS_1996 {

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

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

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

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

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

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

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

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

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

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

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

  108.         /** {@inheritDoc} */
  109.         @Override
  110.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
  111.                                                                  final TimeScales timeScales) {
  112.             try (InputStream in = getStream(NUTATION_ARGUMENTS)) {
  113.                 return new FundamentalNutationArguments(this, timeScale,
  114.                         in, NUTATION_ARGUMENTS, timeScales);
  115.             } catch (IOException e) {
  116.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  117.             }
  118.         }

  120.         /** {@inheritDoc} */
  121.         @Override
  122.         public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {

  124.             // value from chapter 5, page 22
  125.             final PolynomialNutation epsilonA =
  126.                     new PolynomialNutation(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  127.                                              -46.8150   * Constants.ARC_SECONDS_TO_RADIANS,
  128.                                               -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  129.                                                0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  131.             return new TimeScalarFunction() {

  133.                 /** {@inheritDoc} */
  134.                 @Override
  135.                 public double value(final AbsoluteDate date) {
  136.                     return epsilonA.value(evaluateTC(date, timeScales));
  137.                 }

  139.                 /** {@inheritDoc} */
  140.                 @Override
  141.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  142.                     return epsilonA.value(evaluateTC(date, timeScales));
  143.                 }

  145.             };

  147.         }

  149.         /** {@inheritDoc} */
  150.         @Override
  151.         public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {

  153.             // set up nutation arguments
  154.             final FundamentalNutationArguments arguments =
  155.                     getNutationArguments(timeScales);

  157.             // X = 2004.3109″t - 0.42665″t² - 0.198656″t³ + 0.0000140″t⁴
  158.             //     + 0.00006″t² cos Ω + sin ε0 { Σ [(Ai + Ai' t) sin(ARGUMENT) + Ai'' t cos(ARGUMENT)]}
  159.             //     + 0.00204″t² sin Ω + 0.00016″t² sin 2(F - D + Ω),
  160.             final PolynomialNutation xPolynomial =
  161.                     new PolynomialNutation(0,
  162.                                            2004.3109 * Constants.ARC_SECONDS_TO_RADIANS,
  163.                                            -0.42665  * Constants.ARC_SECONDS_TO_RADIANS,
  164.                                            -0.198656 * Constants.ARC_SECONDS_TO_RADIANS,
  165.                                            0.0000140 * Constants.ARC_SECONDS_TO_RADIANS);

  167.             final double fXCosOm    = 0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
  168.             final double fXSinOm    = 0.00204 * Constants.ARC_SECONDS_TO_RADIANS;
  169.             final double fXSin2FDOm = 0.00016 * Constants.ARC_SECONDS_TO_RADIANS;
  170.             final double sinEps0   = FastMath.sin(getMeanObliquityFunction(timeScales)
  171.                     .value(getNutationReferenceEpoch(timeScales)));

  173.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  174.             final PoissonSeriesParser baseParser =
  175.                     new PoissonSeriesParser(12).withFirstDelaunay(1);

  177.             final PoissonSeriesParser xParser =
  178.                     baseParser.
  179.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  180.                     withSinCos(1, 8, deciMilliAS,  9, deciMilliAS);
  181.             final PoissonSeries xSum;
  182.             try (InputStream in = getStream(X_Y_SERIES)) {
  183.                 xSum = xParser.parse(in, X_Y_SERIES);
  184.             } catch (IOException e) {
  185.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  186.             }

  188.             // Y = -0.00013″ - 22.40992″t² + 0.001836″t³ + 0.0011130″t⁴
  189.             //     + Σ [(Bi + Bi' t) cos(ARGUMENT) + Bi'' t sin(ARGUMENT)]
  190.             //    - 0.00231″t² cos Ω − 0.00014″t² cos 2(F - D + Ω)
  191.             final PolynomialNutation yPolynomial =
  192.                     new PolynomialNutation(-0.00013  * Constants.ARC_SECONDS_TO_RADIANS,
  193.                                            0.0,
  194.                                            -22.40992 * Constants.ARC_SECONDS_TO_RADIANS,
  195.                                            0.001836  * Constants.ARC_SECONDS_TO_RADIANS,
  196.                                            0.0011130 * Constants.ARC_SECONDS_TO_RADIANS);

  198.             final double fYCosOm    = -0.00231 * Constants.ARC_SECONDS_TO_RADIANS;
  199.             final double fYCos2FDOm = -0.00014 * Constants.ARC_SECONDS_TO_RADIANS;

  201.             final PoissonSeriesParser yParser =
  202.                     baseParser.
  203.                     withSinCos(0, -1, deciMilliAS, 10, deciMilliAS).
  204.                     withSinCos(1, 12, deciMilliAS, 11, deciMilliAS);
  205.             final PoissonSeries ySum;
  206.             try (InputStream in = getStream(X_Y_SERIES)) {
  207.                 ySum = yParser.parse(in, X_Y_SERIES);
  208.             } catch (IOException e) {
  209.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  210.             }

  212.             final PoissonSeries.CompiledSeries xySum =
  213.                     PoissonSeries.compile(xSum, ySum);

  215.             // s = -XY/2 + 0.00385″t - 0.07259″t³ - 0.00264″ sin Ω - 0.00006″ sin 2Ω
  216.             //     + 0.00074″t² sin Ω + 0.00006″t² sin 2(F - D + Ω)
  217.             final double fST          =  0.00385 * Constants.ARC_SECONDS_TO_RADIANS;
  218.             final double fST3         = -0.07259 * Constants.ARC_SECONDS_TO_RADIANS;
  219.             final double fSSinOm      = -0.00264 * Constants.ARC_SECONDS_TO_RADIANS;
  220.             final double fSSin2Om     = -0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
  221.             final double fST2SinOm    =  0.00074 * Constants.ARC_SECONDS_TO_RADIANS;
  222.             final double fST2Sin2FDOm =  0.00006 * Constants.ARC_SECONDS_TO_RADIANS;

  224.             return new TimeVectorFunction() {

  226.                 /** {@inheritDoc} */
  227.                 @Override
  228.                 public double[] value(final AbsoluteDate date) {

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

  233.                     final double omega     = elements.getOmega();
  234.                     final double f         = elements.getF();
  235.                     final double d         = elements.getD();
  236.                     final double t         = elements.getTC();

  238.                     final double cosOmega  = FastMath.cos(omega);
  239.                     final double sinOmega  = FastMath.sin(omega);
  240.                     final double sin2Omega = FastMath.sin(2 * omega);
  241.                     final double cos2FDOm  = FastMath.cos(2 * (f - d + omega));
  242.                     final double sin2FDOm  = FastMath.sin(2 * (f - d + omega));

  244.                     final double x = xPolynomial.value(t) + sinEps0 * xy[0] +
  245.                             t * t * (fXCosOm * cosOmega + fXSinOm * sinOmega + fXSin2FDOm * cos2FDOm);
  246.                     final double y = yPolynomial.value(t) + xy[1] +
  247.                             t * t * (fYCosOm * cosOmega + fYCos2FDOm * cos2FDOm);
  248.                     final double sPxy2 = fSSinOm * sinOmega + fSSin2Om * sin2Omega +
  249.                             t * (fST + t * (fST2SinOm * sinOmega + fST2Sin2FDOm * sin2FDOm + t * fST3));

  251.                     return new double[] {
  252.                         x, y, sPxy2
  253.                     };

  255.                 }

  257.                 /** {@inheritDoc} */
  258.                 @Override
  259.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

  264.                     final T omega     = elements.getOmega();
  265.                     final T f         = elements.getF();
  266.                     final T d         = elements.getD();
  267.                     final T t         = elements.getTC();
  268.                     final T t2        = t.multiply(t);

  270.                     final T cosOmega  = omega.cos();
  271.                     final T sinOmega  = omega.sin();
  272.                     final T sin2Omega = omega.multiply(2).sin();
  273.                     final T fMDPO2 = f.subtract(d).add(omega).multiply(2);
  274.                     final T cos2FDOm  = fMDPO2.cos();
  275.                     final T sin2FDOm  = fMDPO2.sin();

  277.                     final T x = xPolynomial.value(t).
  278.                                 add(xy[0].multiply(sinEps0)).
  279.                                 add(t2.multiply(cosOmega.multiply(fXCosOm).add(sinOmega.multiply(fXSinOm)).add(cos2FDOm.multiply(fXSin2FDOm))));
  280.                     final T y = yPolynomial.value(t).
  281.                                 add(xy[1]).
  282.                                 add(t2.multiply(cosOmega.multiply(fYCosOm).add(cos2FDOm.multiply(fYCos2FDOm))));
  283.                     final T sPxy2 = sinOmega.multiply(fSSinOm).
  284.                                     add(sin2Omega.multiply(fSSin2Om)).
  285.                                     add(t.multiply(fST3).add(sinOmega.multiply(fST2SinOm)).add(sin2FDOm.multiply(fST2Sin2FDOm)).multiply(t).add(fST).multiply(t));

  287.                     final T[] a = MathArrays.buildArray(date.getField(), 3);
  288.                     a[0] = x;
  289.                     a[1] = y;
  290.                     a[2] = sPxy2;
  291.                     return a;

  293.                 }

  295.             };

  297.         }

  299.         /** {@inheritDoc} */
  300.         @Override
  301.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

  303.             // set up the conventional polynomials
  304.             // the following values are from Lieske et al. paper:
  305.             // Expressions for the precession quantities based upon the IAU(1976) system of astronomical constants
  306.             // http://articles.adsabs.harvard.edu/full/1977A%26A....58....1L
  307.             // also available as equation 30 in IERS 2003 conventions
  308.             final PolynomialNutation psiA =
  309.                     new PolynomialNutation(   0.0,
  310.                                            5038.7784   * Constants.ARC_SECONDS_TO_RADIANS,
  311.                                              -1.07259  * Constants.ARC_SECONDS_TO_RADIANS,
  312.                                              -0.001147 * Constants.ARC_SECONDS_TO_RADIANS);
  313.             final PolynomialNutation omegaA =
  314.                     new PolynomialNutation(getMeanObliquityFunction(timeScales)
  315.                             .value(getNutationReferenceEpoch(timeScales)),
  316.                                             0.0,
  317.                                             0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
  318.                                            -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
  319.             final PolynomialNutation chiA =
  320.                     new PolynomialNutation( 0.0,
  321.                                            10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
  322.                                            -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
  323.                                            -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);

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

  327.         }

  329.         /** {@inheritDoc} */
  330.         @Override
  331.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  333.             // set up nutation arguments
  334.             final FundamentalNutationArguments arguments =
  335.                     getNutationArguments(timeScales);

  337.             // set up Poisson series
  338.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  339.             final PoissonSeriesParser baseParser =
  340.                     new PoissonSeriesParser(10).withFirstDelaunay(1);

  342.             final PoissonSeriesParser psiParser =
  343.                     baseParser.
  344.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  345.                     withSinCos(1, 8, deciMilliAS, -1, deciMilliAS);
  346.             final PoissonSeries psiSeries;
  347.             try (InputStream in = getStream(PSI_EPSILON_SERIES)) {
  348.                 psiSeries = psiParser.parse(in, PSI_EPSILON_SERIES);
  349.             } catch (IOException e) {
  350.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  351.             }

  353.             final PoissonSeriesParser epsilonParser =
  354.                     baseParser.
  355.                     withSinCos(0, -1, deciMilliAS, 9, deciMilliAS).
  356.                     withSinCos(1, -1, deciMilliAS, 10, deciMilliAS);
  357.             final PoissonSeries epsilonSeries;
  358.             try (InputStream in = getStream(PSI_EPSILON_SERIES)) {
  359.                 epsilonSeries = epsilonParser.parse(in, PSI_EPSILON_SERIES);
  360.             } catch (IOException e) {
  361.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  362.             }

  364.             final PoissonSeries.CompiledSeries psiEpsilonSeries =
  365.                     PoissonSeries.compile(psiSeries, epsilonSeries);

  367.             return new TimeVectorFunction() {

  369.                 /** {@inheritDoc} */
  370.                 @Override
  371.                 public double[] value(final AbsoluteDate date) {
  372.                     final BodiesElements elements = arguments.evaluateAll(date);
  373.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  374.                     return new double[] {
  375.                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  376.                     };
  377.                 }

  379.                 /** {@inheritDoc} */
  380.                 @Override
  381.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  382.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  383.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  384.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  385.                     result[0] = psiEpsilon[0];
  386.                     result[1] = psiEpsilon[1];
  387.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  388.                     return result;
  389.                 }

  391.             };

  393.         }

  395.         /** {@inheritDoc} */
  396.         @Override
  397.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  398.                                                   final TimeScales timeScales) {

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

  403.             // constants from IERS 1996 page 21
  404.             // the underlying model is IAU 1982 GMST-UT1
  405.             final AbsoluteDate gmstReference = new AbsoluteDate(
  406.                     DateComponents.J2000_EPOCH, TimeComponents.H12, timeScales.getTAI());
  407.             final double gmst0 = 24110.54841;
  408.             final double gmst1 = 8640184.812866;
  409.             final double gmst2 = 0.093104;
  410.             final double gmst3 = -6.2e-6;

  412.             return new TimeScalarFunction() {

  414.                 /** {@inheritDoc} */
  415.                 @Override
  416.                 public double value(final AbsoluteDate date) {

  418.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  419.                     final double dtai = date.durationFrom(gmstReference);
  420.                     final double tut1 = dtai + ut1.offsetFromTAI(date);
  421.                     final double tt   = tut1 / Constants.JULIAN_CENTURY;

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

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

  430.                 }

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

  436.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  437.                     final T dtai = date.durationFrom(gmstReference);
  438.                     final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
  439.                     final T tt   = tut1.divide(Constants.JULIAN_CENTURY);

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

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

  448.                 }

  450.             };

  452.         }

  454.         /** {@inheritDoc} */
  455.         @Override
  456.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  457.                                                       final TimeScales timeScales) {

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

  462.             // constants from IERS 1996 page 21
  463.             // the underlying model is IAU 1982 GMST-UT1
  464.             final AbsoluteDate gmstReference = new AbsoluteDate(
  465.                     DateComponents.J2000_EPOCH, TimeComponents.H12, timeScales.getTAI());
  466.             final double gmst1 = 8640184.812866;
  467.             final double gmst2 = 0.093104;
  468.             final double gmst3 = -6.2e-6;

  470.             return new TimeScalarFunction() {

  472.                 /** {@inheritDoc} */
  473.                 @Override
  474.                 public double value(final AbsoluteDate date) {

  476.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  477.                     final double dtai = date.durationFrom(gmstReference);
  478.                     final double tut1 = dtai + ut1.offsetFromTAI(date);
  479.                     final double tt   = tut1 / Constants.JULIAN_CENTURY;

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

  484.                 }

  486.                 /** {@inheritDoc} */
  487.                 @Override
  488.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  490.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  491.                     final T dtai = date.durationFrom(gmstReference);
  492.                     final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
  493.                     final T tt   = tut1.divide(Constants.JULIAN_CENTURY);

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

  498.                 }

  500.             };

  502.         }

  504.         /** {@inheritDoc} */
  505.         @Override
  506.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  507.                                                   final EOPHistory eopHistory,
  508.                                                   final TimeScales timeScales) {

  510.             // obliquity
  511.             final TimeScalarFunction epsilonA = getMeanObliquityFunction(timeScales);

  513.             // GMST function
  514.             final TimeScalarFunction gmst = getGMSTFunction(ut1, timeScales);

  516.             // nutation function
  517.             final TimeVectorFunction nutation = getNutationFunction(timeScales);

  519.             return new TimeScalarFunction() {

  521.                 /** {@inheritDoc} */
  522.                 @Override
  523.                 public double value(final AbsoluteDate date) {

  525.                     // compute equation of equinoxes
  526.                     final double[] angles = nutation.value(date);
  527.                     double deltaPsi = angles[0];
  528.                     if (eopHistory != null) {
  529.                         deltaPsi += eopHistory.getEquinoxNutationCorrection(date)[0];
  530.                     }
  531.                     final double eqe = deltaPsi  * FastMath.cos(epsilonA.value(date)) + angles[2];

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

  536.                 }

  538.                 /** {@inheritDoc} */
  539.                 @Override
  540.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  542.                     // compute equation of equinoxes
  543.                     final T[] angles = nutation.value(date);
  544.                     T deltaPsi = angles[0];
  545.                     if (eopHistory != null) {
  546.                         deltaPsi = deltaPsi.add(eopHistory.getEquinoxNutationCorrection(date)[0]);
  547.                     }
  548.                     final T eqe = deltaPsi.multiply(epsilonA.value(date).cos()).add(angles[2]);

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

  553.                 }

  555.             };

  557.         }

  559.         /** {@inheritDoc} */
  560.         @Override
  561.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

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

  571.             // set up Poisson series
  572.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  573.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(17).
  574.                     withOptionalColumn(1).
  575.                     withGamma(7).
  576.                     withFirstDelaunay(2);
  577.             final PoissonSeries xSeries;
  578.             try (InputStream in = getStream(TIDAL_CORRECTION_XP_YP_SERIES)) {
  579.                 xSeries = xyParser.
  580.                         withSinCos(0, 14, milliAS, 15, milliAS).
  581.                         parse(in, TIDAL_CORRECTION_XP_YP_SERIES);
  582.             } catch (IOException e) {
  583.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  584.             }
  585.             final PoissonSeries ySeries;
  586.             try (InputStream in = getStream(TIDAL_CORRECTION_XP_YP_SERIES)) {
  587.                 ySeries = xyParser.
  588.                         withSinCos(0, 16, milliAS, 17, milliAS).
  589.                         parse(in, TIDAL_CORRECTION_XP_YP_SERIES);
  590.             } catch (IOException e) {
  591.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  592.             }

  594.             final double deciMilliS = 1.0e-4;
  595.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(17).
  596.                     withOptionalColumn(1).
  597.                     withGamma(7).
  598.                     withFirstDelaunay(2).
  599.                     withSinCos(0, 16, deciMilliS, 17, deciMilliS);
  600.             final PoissonSeries ut1Series;
  601.             try (InputStream in = getStream(TIDAL_CORRECTION_UT1_SERIES)) {
  602.                 ut1Series = ut1Parser.parse(in, TIDAL_CORRECTION_UT1_SERIES);
  603.             } catch (IOException e) {
  604.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  605.             }

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

  609.         }

  611.         /** {@inheritDoc} */
  612.         public LoveNumbers getLoveNumbers() {
  613.             return loadLoveNumbers(LOVE_NUMBERS);
  614.         }

  616.         /** {@inheritDoc} */
  617.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  618.                                                                      final TimeScales timeScales) {

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

  623.             // set up Poisson series
  624.             final PoissonSeriesParser k20Parser =
  625.                     new PoissonSeriesParser(18).
  626.                         withOptionalColumn(1).
  627.                         withDoodson(4, 2).
  628.                         withFirstDelaunay(10);
  629.             final PoissonSeriesParser k21Parser =
  630.                     new PoissonSeriesParser(18).
  631.                         withOptionalColumn(1).
  632.                         withDoodson(4, 3).
  633.                         withFirstDelaunay(10);
  634.             final PoissonSeriesParser k22Parser =
  635.                     new PoissonSeriesParser(16).
  636.                         withOptionalColumn(1).
  637.                         withDoodson(4, 2).
  638.                         withFirstDelaunay(10);

  640.             final double pico = 1.0e-12;
  641.             try (InputStream k20 = getStream(K20_FREQUENCY_DEPENDENCE);
  642.                  InputStream k21In1 = getStream(K21_FREQUENCY_DEPENDENCE);
  643.                  InputStream k21In2 = getStream(K21_FREQUENCY_DEPENDENCE);
  644.                  InputStream k22In1 = getStream(K22_FREQUENCY_DEPENDENCE);
  645.                  InputStream k22In2 = getStream(K22_FREQUENCY_DEPENDENCE)) {
  646.                 final PoissonSeries c20Series =
  647.                         k20Parser.
  648.                         withSinCos(0, 18, -pico, 16, pico).
  649.                         parse(k20, K20_FREQUENCY_DEPENDENCE);
  650.                 final PoissonSeries c21Series =
  651.                         k21Parser.
  652.                         withSinCos(0, 17, pico, 18, pico).
  653.                         parse(k21In1, K21_FREQUENCY_DEPENDENCE);
  654.                 final PoissonSeries s21Series =
  655.                         k21Parser.
  656.                         withSinCos(0, 18, -pico, 17, pico).
  657.                         parse(k21In2, K21_FREQUENCY_DEPENDENCE);
  658.                 final PoissonSeries c22Series =
  659.                         k22Parser.
  660.                         withSinCos(0, -1, pico, 16, pico).
  661.                         parse(k22In1, K22_FREQUENCY_DEPENDENCE);
  662.                 final PoissonSeries s22Series =
  663.                         k22Parser.
  664.                         withSinCos(0, 16, -pico, -1, pico).
  665.                         parse(k22In2, K22_FREQUENCY_DEPENDENCE);

  667.                 return new TideFrequencyDependenceFunction(arguments,
  668.                                                            c20Series,
  669.                                                            c21Series, s21Series,
  670.                                                            c22Series, s22Series);
  671.             } catch (IOException e) {
  672.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  673.             }

  675.         }

  677.         /** {@inheritDoc} */
  678.         @Override
  679.         public double getPermanentTide() {
  680.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  681.         }

  683.         /** {@inheritDoc} */
  684.         @Override
  685.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  687.             // constants from IERS 1996 page 47
  688.             final double globalFactor = -1.348e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  689.             final double coupling     =  0.00112;

  691.             return new TimeVectorFunction() {

  693.                 /** {@inheritDoc} */
  694.                 @Override
  695.                 public double[] value(final AbsoluteDate date) {
  696.                     final PoleCorrection pole = eopHistory.getPoleCorrection(date);
  697.                     return new double[] {
  698.                         globalFactor * (pole.getXp() + coupling * pole.getYp()),
  699.                         globalFactor * (coupling * pole.getXp() - pole.getYp()),
  700.                     };
  701.                 }

  703.                 /** {@inheritDoc} */
  704.                 @Override
  705.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  706.                     final FieldPoleCorrection<T> pole = eopHistory.getPoleCorrection(date);
  707.                     final T[] a = MathArrays.buildArray(date.getField(), 2);
  708.                     a[0] = pole.getXp().add(pole.getYp().multiply(coupling)).multiply(globalFactor);
  709.                     a[1] = pole.getXp().multiply(coupling).subtract(pole.getYp()).multiply(globalFactor);
  710.                     return a;
  711.                 }

  713.             };
  714.         }

  716.         /** {@inheritDoc} */
  717.         @Override
  718.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {

  720.             return new TimeVectorFunction() {

  722.                 /** {@inheritDoc} */
  723.                 @Override
  724.                 public double[] value(final AbsoluteDate date) {
  725.                     // there are no model for ocean pole tide prior to conventions 2010
  726.                     return new double[] {
  727.                         0.0, 0.0
  728.                     };
  729.                 }

  731.                 /** {@inheritDoc} */
  732.                 @Override
  733.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  734.                     // there are no model for ocean pole tide prior to conventions 2010
  735.                     return MathArrays.buildArray(date.getField(), 2);
  736.                 }

  738.             };
  739.         }

  741.         /** {@inheritDoc} */
  742.         @Override
  743.         public double[] getNominalTidalDisplacement() {

  745.             //  // elastic Earth values
  746.             //  return new double[] {
  747.             //      // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  748.             //      0.6026,                                  -0.0006, 0.292, -0.0025, -0.0022,
  749.             //      // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  750.             //      0.0831,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  751.             //      // H₀
  752.             //      -0.31460
  753.             //  };

  755.             // anelastic Earth values
  756.             return new double[] {
  757.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  758.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  759.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  760.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  761.                 // H₀
  762.                 -0.31460
  763.             };

  765.         }

  767.         /** {@inheritDoc} */
  768.         @Override
  769.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
  770.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  771.                                                                   18, 17, -1, 18, -1);
  772.         }

  774.         /** {@inheritDoc} */
  775.         @Override
  776.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  777.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  778.                                                                 20, 17, 19, 18, 20);
  779.         }

  781.     },

  783.     /** Constant for IERS 2003 conventions. */
  784.     IERS_2003 {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  834.         /** {@inheritDoc} */
  835.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
  836.                                                                  final TimeScales timeScales) {
  837.             try (InputStream in = getStream(NUTATION_ARGUMENTS)) {
  838.                 return new FundamentalNutationArguments(this, timeScale,
  839.                         in, NUTATION_ARGUMENTS, timeScales);
  840.             } catch (IOException e) {
  841.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  842.             }
  843.         }

  845.         /** {@inheritDoc} */
  846.         @Override
  847.         public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {

  849.             // epsilon 0 value from chapter 5, page 41, other terms from equation 32 page 45
  850.             final PolynomialNutation epsilonA =
  851.                     new PolynomialNutation(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  852.                                              -46.84024  * Constants.ARC_SECONDS_TO_RADIANS,
  853.                                               -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  854.                                                0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  856.             return new TimeScalarFunction() {

  858.                 /** {@inheritDoc} */
  859.                 @Override
  860.                 public double value(final AbsoluteDate date) {
  861.                     return epsilonA.value(evaluateTC(date, timeScales));
  862.                 }

  864.                 /** {@inheritDoc} */
  865.                 @Override
  866.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  867.                     return epsilonA.value(evaluateTC(date, timeScales));
  868.                 }

  870.             };

  872.         }

  874.         /** {@inheritDoc} */
  875.         @Override
  876.         public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {

  878.             // set up nutation arguments
  879.             final FundamentalNutationArguments arguments =
  880.                     getNutationArguments(timeScales);

  882.             // set up Poisson series
  883.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  884.             final PoissonSeriesParser parser =
  885.                     new PoissonSeriesParser(17).
  886.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  887.                         withFirstDelaunay(4).
  888.                         withFirstPlanetary(9).
  889.                         withSinCos(0, 2, microAS, 3, microAS);

  891.             final PoissonSeries.CompiledSeries xys;
  892.             try (InputStream xIn = getStream(X_SERIES);
  893.                  InputStream yIn = getStream(Y_SERIES);
  894.                  InputStream sIn = getStream(S_SERIES)) {
  895.                 final PoissonSeries xSeries = parser.parse(xIn, X_SERIES);
  896.                 final PoissonSeries ySeries = parser.parse(yIn, Y_SERIES);
  897.                 final PoissonSeries sSeries = parser.parse(sIn, S_SERIES);
  898.                 xys = PoissonSeries.compile(xSeries, ySeries, sSeries);
  899.             } catch (IOException e) {
  900.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  901.             }

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

  906.                 /** {@inheritDoc} */
  907.                 @Override
  908.                 public double[] value(final AbsoluteDate date) {
  909.                     return xys.value(arguments.evaluateAll(date));
  910.                 }

  912.                 /** {@inheritDoc} */
  913.                 @Override
  914.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  915.                     return xys.value(arguments.evaluateAll(date));
  916.                 }

  918.             };

  920.         }


  923.         /** {@inheritDoc} */
  924.         @Override
  925.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

  927.             // set up the conventional polynomials
  928.             // the following values are from equation 32 in IERS 2003 conventions
  929.             final PolynomialNutation psiA =
  930.                     new PolynomialNutation(    0.0,
  931.                                             5038.47875   * Constants.ARC_SECONDS_TO_RADIANS,
  932.                                               -1.07259   * Constants.ARC_SECONDS_TO_RADIANS,
  933.                                               -0.001147  * Constants.ARC_SECONDS_TO_RADIANS);
  934.             final PolynomialNutation omegaA =
  935.                     new PolynomialNutation(getMeanObliquityFunction(timeScales)
  936.                             .value(getNutationReferenceEpoch(timeScales)),
  937.                                            -0.02524   * Constants.ARC_SECONDS_TO_RADIANS,
  938.                                             0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
  939.                                            -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
  940.             final PolynomialNutation chiA =
  941.                     new PolynomialNutation( 0.0,
  942.                                            10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
  943.                                            -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
  944.                                            -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);

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

  948.         }

  950.         /** {@inheritDoc} */
  951.         @Override
  952.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  954.             // set up nutation arguments
  955.             final FundamentalNutationArguments arguments =
  956.                     getNutationArguments(timeScales);

  958.             // set up Poisson series
  959.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  960.             final PoissonSeriesParser luniSolarParser =
  961.                     new PoissonSeriesParser(14).withFirstDelaunay(1);
  962.             final PoissonSeriesParser luniSolarPsiParser =
  963.                     luniSolarParser.
  964.                     withSinCos(0, 7, milliAS, 11, milliAS).
  965.                     withSinCos(1, 8, milliAS, 12, milliAS);
  966.             final PoissonSeries psiLuniSolarSeries;
  967.             try (InputStream in = getStream(LUNI_SOLAR_SERIES)) {
  968.                 psiLuniSolarSeries = luniSolarPsiParser.parse(in, LUNI_SOLAR_SERIES);
  969.             } catch (IOException e) {
  970.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  971.             }
  972.             final PoissonSeriesParser luniSolarEpsilonParser =
  973.                     luniSolarParser.
  974.                     withSinCos(0, 13, milliAS, 9, milliAS).
  975.                     withSinCos(1, 14, milliAS, 10, milliAS);
  976.             final PoissonSeries epsilonLuniSolarSeries;
  977.             try (InputStream in = getStream(LUNI_SOLAR_SERIES)) {
  978.                 epsilonLuniSolarSeries = luniSolarEpsilonParser.parse(in, LUNI_SOLAR_SERIES);
  979.             } catch (IOException e) {
  980.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  981.             }

  983.             final PoissonSeriesParser planetaryParser =
  984.                     new PoissonSeriesParser(21).
  985.                         withFirstDelaunay(2).
  986.                         withFirstPlanetary(7);
  987.             final PoissonSeriesParser planetaryPsiParser =
  988.                     planetaryParser.withSinCos(0, 17, milliAS, 18, milliAS);
  989.             final PoissonSeries psiPlanetarySeries;
  990.             try (InputStream in = getStream(PLANETARY_SERIES)) {
  991.                 psiPlanetarySeries = planetaryPsiParser.parse(in, PLANETARY_SERIES);
  992.             } catch (IOException e) {
  993.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  994.             }
  995.             final PoissonSeriesParser planetaryEpsilonParser =
  996.                     planetaryParser.withSinCos(0, 19, milliAS, 20, milliAS);
  997.             final PoissonSeries epsilonPlanetarySeries;
  998.             try (InputStream in = getStream(PLANETARY_SERIES)) {
  999.                 epsilonPlanetarySeries = planetaryEpsilonParser.parse(in, PLANETARY_SERIES);
  1000.             } catch (IOException e) {
  1001.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1002.             }

  1004.             final PoissonSeries.CompiledSeries luniSolarSeries =
  1005.                     PoissonSeries.compile(psiLuniSolarSeries, epsilonLuniSolarSeries);
  1006.             final PoissonSeries.CompiledSeries planetarySeries =
  1007.                     PoissonSeries.compile(psiPlanetarySeries, epsilonPlanetarySeries);

  1009.             return new TimeVectorFunction() {

  1011.                 /** {@inheritDoc} */
  1012.                 @Override
  1013.                 public double[] value(final AbsoluteDate date) {
  1014.                     final BodiesElements elements = arguments.evaluateAll(date);
  1015.                     final double[] luniSolar = luniSolarSeries.value(elements);
  1016.                     final double[] planetary = planetarySeries.value(elements);
  1017.                     return new double[] {
  1018.                         luniSolar[0] + planetary[0], luniSolar[1] + planetary[1],
  1019.                         IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  1020.                     };
  1021.                 }

  1023.                 /** {@inheritDoc} */
  1024.                 @Override
  1025.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1026.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1027.                     final T[] luniSolar = luniSolarSeries.value(elements);
  1028.                     final T[] planetary = planetarySeries.value(elements);
  1029.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  1030.                     result[0] = luniSolar[0].add(planetary[0]);
  1031.                     result[1] = luniSolar[1].add(planetary[1]);
  1032.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  1033.                     return result;
  1034.                 }

  1036.             };

  1038.         }

  1040.         /** {@inheritDoc} */
  1041.         @Override
  1042.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  1043.                                                   final TimeScales timeScales) {

  1045.             // Earth Rotation Angle
  1046.             final StellarAngleCapitaine era =
  1047.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  1049.             // Polynomial part of the apparent sidereal time series
  1050.             // which is the opposite of Equation of Origins (EO)
  1051.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1052.             final PoissonSeriesParser parser =
  1053.                     new PoissonSeriesParser(17).
  1054.                         withFirstDelaunay(4).
  1055.                         withFirstPlanetary(9).
  1056.                         withSinCos(0, 2, microAS, 3, microAS).
  1057.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1058.             final PolynomialNutation minusEO;
  1059.             try (InputStream in = getStream(GST_SERIES)) {
  1060.                 minusEO = parser.parse(in, GST_SERIES).getPolynomial();
  1061.             } catch (IOException e) {
  1062.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1063.             }

  1065.             // create a function evaluating the series
  1066.             return new TimeScalarFunction() {

  1068.                 /** {@inheritDoc} */
  1069.                 @Override
  1070.                 public double value(final AbsoluteDate date) {
  1071.                     return era.value(date) + minusEO.value(evaluateTC(date, timeScales));
  1072.                 }

  1074.                 /** {@inheritDoc} */
  1075.                 @Override
  1076.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1077.                     return era.value(date).add(minusEO.value(evaluateTC(date, timeScales)));
  1078.                 }

  1080.             };

  1082.         }

  1084.         /** {@inheritDoc} */
  1085.         @Override
  1086.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  1087.                                                       final TimeScales timeScales) {

  1089.             // Earth Rotation Angle
  1090.             final StellarAngleCapitaine era =
  1091.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  1093.             // Polynomial part of the apparent sidereal time series
  1094.             // which is the opposite of Equation of Origins (EO)
  1095.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1096.             final PoissonSeriesParser parser =
  1097.                     new PoissonSeriesParser(17).
  1098.                         withFirstDelaunay(4).
  1099.                         withFirstPlanetary(9).
  1100.                         withSinCos(0, 2, microAS, 3, microAS).
  1101.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1102.             final PolynomialNutation minusEO;
  1103.             try (InputStream in = getStream(GST_SERIES)) {
  1104.                 minusEO = parser.parse(in, GST_SERIES).getPolynomial();
  1105.             } catch (IOException e) {
  1106.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1107.             }

  1109.             // create a function evaluating the series
  1110.             return new TimeScalarFunction() {

  1112.                 /** {@inheritDoc} */
  1113.                 @Override
  1114.                 public double value(final AbsoluteDate date) {
  1115.                     return era.getRate() + minusEO.derivative(evaluateTC(date, timeScales));
  1116.                 }

  1118.                 /** {@inheritDoc} */
  1119.                 @Override
  1120.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1121.                     return minusEO.derivative(evaluateTC(date, timeScales)).add(era.getRate());
  1122.                 }

  1124.             };

  1126.         }

  1128.         /** {@inheritDoc} */
  1129.         @Override
  1130.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  1131.                                                   final EOPHistory eopHistory,
  1132.                                                   final TimeScales timeScales) {

  1134.             // set up nutation arguments
  1135.             final FundamentalNutationArguments arguments =
  1136.                     getNutationArguments(timeScales);

  1138.             // mean obliquity function
  1139.             final TimeScalarFunction epsilon = getMeanObliquityFunction(timeScales);

  1141.             // set up Poisson series
  1142.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  1143.             final PoissonSeriesParser luniSolarPsiParser =
  1144.                     new PoissonSeriesParser(14).
  1145.                         withFirstDelaunay(1).
  1146.                         withSinCos(0, 7, milliAS, 11, milliAS).
  1147.                         withSinCos(1, 8, milliAS, 12, milliAS);
  1148.             final PoissonSeries psiLuniSolarSeries;
  1149.             try (InputStream in = getStream(LUNI_SOLAR_SERIES)) {
  1150.                 psiLuniSolarSeries = luniSolarPsiParser.parse(in, LUNI_SOLAR_SERIES);
  1151.             } catch (IOException e) {
  1152.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1153.             }

  1155.             final PoissonSeriesParser planetaryPsiParser =
  1156.                     new PoissonSeriesParser(21).
  1157.                         withFirstDelaunay(2).
  1158.                         withFirstPlanetary(7).
  1159.                         withSinCos(0, 17, milliAS, 18, milliAS);
  1160.             final PoissonSeries psiPlanetarySeries;
  1161.             try (InputStream in = getStream(PLANETARY_SERIES)) {
  1162.                 psiPlanetarySeries = planetaryPsiParser.parse(in, PLANETARY_SERIES);
  1163.             } catch (IOException e) {
  1164.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1165.             }


  1168.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1169.             final PoissonSeriesParser gstParser =
  1170.                     new PoissonSeriesParser(17).
  1171.                         withFirstDelaunay(4).
  1172.                         withFirstPlanetary(9).
  1173.                         withSinCos(0, 2, microAS, 3, microAS).
  1174.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1175.             final PoissonSeries gstSeries;
  1176.             try (InputStream in = getStream(GST_SERIES)) {
  1177.                 gstSeries = gstParser.parse(in, GST_SERIES);
  1178.             } catch (IOException e) {
  1179.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1180.             }
  1181.             final PoissonSeries.CompiledSeries psiGstSeries =
  1182.                     PoissonSeries.compile(psiLuniSolarSeries, psiPlanetarySeries, gstSeries);

  1184.             // ERA function
  1185.             final TimeScalarFunction era =
  1186.                     getEarthOrientationAngleFunction(ut1, timeScales.getTAI());

  1188.             return new TimeScalarFunction() {

  1190.                 /** {@inheritDoc} */
  1191.                 @Override
  1192.                 public double value(final AbsoluteDate date) {

  1194.                     // evaluate equation of origins
  1195.                     final BodiesElements elements = arguments.evaluateAll(date);
  1196.                     final double[] angles = psiGstSeries.value(elements);
  1197.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  1198.                     final double deltaPsi = angles[0] + angles[1] + ddPsi;
  1199.                     final double epsilonA = epsilon.value(date);

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

  1205.                 }

  1207.                 /** {@inheritDoc} */
  1208.                 @Override
  1209.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  1211.                     // evaluate equation of origins
  1212.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1213.                     final T[] angles = psiGstSeries.value(elements);
  1214.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  1215.                     final T deltaPsi = angles[0].add(angles[1]).add(ddPsi);
  1216.                     final T epsilonA = epsilon.value(date);

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

  1222.                 }

  1224.             };

  1226.         }

  1228.         /** {@inheritDoc} */
  1229.         @Override
  1230.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

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

  1240.             // set up Poisson series
  1241.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1242.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
  1243.                     withOptionalColumn(1).
  1244.                     withGamma(2).
  1245.                     withFirstDelaunay(3);
  1246.             try (InputStream tidalIn1 = getStream(TIDAL_CORRECTION_XP_YP_SERIES);
  1247.                  InputStream tidalIn2 = getStream(TIDAL_CORRECTION_XP_YP_SERIES);
  1248.                  InputStream tidalUt1 = getStream(TIDAL_CORRECTION_UT1_SERIES)) {
  1249.                 final PoissonSeries xSeries =
  1250.                         xyParser.
  1251.                         withSinCos(0, 10, microAS, 11, microAS).
  1252.                         parse(tidalIn1, TIDAL_CORRECTION_XP_YP_SERIES);
  1253.                 final PoissonSeries ySeries =
  1254.                         xyParser.
  1255.                         withSinCos(0, 12, microAS, 13, microAS).
  1256.                         parse(tidalIn2, TIDAL_CORRECTION_XP_YP_SERIES);

  1258.                 final double microS = 1.0e-6;
  1259.                 final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  1260.                         withOptionalColumn(1).
  1261.                         withGamma(2).
  1262.                         withFirstDelaunay(3).
  1263.                         withSinCos(0, 10, microS, 11, microS);
  1264.                 final PoissonSeries ut1Series =
  1265.                         ut1Parser.parse(tidalUt1, TIDAL_CORRECTION_UT1_SERIES);

  1267.                 return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
  1268.             } catch (IOException e) {
  1269.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1270.             }

  1272.         }

  1274.         /** {@inheritDoc} */
  1275.         public LoveNumbers getLoveNumbers() {
  1276.             return loadLoveNumbers(LOVE_NUMBERS);
  1277.         }

  1279.         /** {@inheritDoc} */
  1280.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  1281.                                                                      final TimeScales timeScales) {

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

  1286.             // set up Poisson series
  1287.             final PoissonSeriesParser k20Parser =
  1288.                     new PoissonSeriesParser(18).
  1289.                         withOptionalColumn(1).
  1290.                         withDoodson(4, 2).
  1291.                         withFirstDelaunay(10);
  1292.             final PoissonSeriesParser k21Parser =
  1293.                     new PoissonSeriesParser(18).
  1294.                         withOptionalColumn(1).
  1295.                         withDoodson(4, 3).
  1296.                         withFirstDelaunay(10);
  1297.             final PoissonSeriesParser k22Parser =
  1298.                     new PoissonSeriesParser(16).
  1299.                         withOptionalColumn(1).
  1300.                         withDoodson(4, 2).
  1301.                         withFirstDelaunay(10);

  1303.             final double pico = 1.0e-12;
  1304.             try (InputStream k20In = getStream(K20_FREQUENCY_DEPENDENCE);
  1305.                  InputStream k21In1 = getStream(K21_FREQUENCY_DEPENDENCE);
  1306.                  InputStream k21In2 = getStream(K21_FREQUENCY_DEPENDENCE);
  1307.                  InputStream k22In1 = getStream(K22_FREQUENCY_DEPENDENCE);
  1308.                  InputStream k22In2 = getStream(K22_FREQUENCY_DEPENDENCE)) {
  1309.                 final PoissonSeries c20Series =
  1310.                         k20Parser.
  1311.                         withSinCos(0, 18, -pico, 16, pico).
  1312.                         parse(k20In, K20_FREQUENCY_DEPENDENCE);
  1313.                 final PoissonSeries c21Series =
  1314.                         k21Parser.
  1315.                         withSinCos(0, 17, pico, 18, pico).
  1316.                         parse(k21In1, K21_FREQUENCY_DEPENDENCE);
  1317.                 final PoissonSeries s21Series =
  1318.                         k21Parser.
  1319.                         withSinCos(0, 18, -pico, 17, pico).
  1320.                         parse(k21In2, K21_FREQUENCY_DEPENDENCE);
  1321.                 final PoissonSeries c22Series =
  1322.                         k22Parser.
  1323.                         withSinCos(0, -1, pico, 16, pico).
  1324.                         parse(k22In1, K22_FREQUENCY_DEPENDENCE);
  1325.                 final PoissonSeries s22Series =
  1326.                         k22Parser.
  1327.                         withSinCos(0, 16, -pico, -1, pico).
  1328.                         parse(k22In2, K22_FREQUENCY_DEPENDENCE);

  1330.                 return new TideFrequencyDependenceFunction(arguments,
  1331.                                                            c20Series,
  1332.                                                            c21Series, s21Series,
  1333.                                                            c22Series, s22Series);
  1334.             } catch (IOException e) {
  1335.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1336.             }

  1338.         }

  1340.         /** {@inheritDoc} */
  1341.         @Override
  1342.         public double getPermanentTide() {
  1343.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1344.         }

  1346.         /** {@inheritDoc} */
  1347.         @Override
  1348.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  1350.             // annual pole from ftp://tai.bipm.org/iers/conv2003/chapter7/annual.pole
  1351.             final TimeScale utc = eopHistory.getTimeScales().getUTC();
  1352.             final SimpleTimeStampedTableParser.RowConverter<MeanPole> converter =
  1353.                 new SimpleTimeStampedTableParser.RowConverter<MeanPole>() {
  1354.                     /** {@inheritDoc} */
  1355.                     @Override
  1356.                     public MeanPole convert(final double[] rawFields) {
  1357.                         return new MeanPole(new AbsoluteDate((int) rawFields[0], 1, 1, utc),
  1358.                                             rawFields[1] * Constants.ARC_SECONDS_TO_RADIANS,
  1359.                                             rawFields[2] * Constants.ARC_SECONDS_TO_RADIANS);
  1360.                     }
  1361.                 };
  1362.             final SimpleTimeStampedTableParser<MeanPole> parser =
  1363.                     new SimpleTimeStampedTableParser<MeanPole>(3, converter);
  1364.             final List<MeanPole> annualPoleList;
  1365.             try (InputStream in = getStream(ANNUAL_POLE)) {
  1366.                 annualPoleList = parser.parse(in, ANNUAL_POLE);
  1367.             } catch (IOException e) {
  1368.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1369.             }
  1370.             final AbsoluteDate firstAnnualPoleDate = annualPoleList.get(0).getDate();
  1371.             final AbsoluteDate lastAnnualPoleDate  = annualPoleList.get(annualPoleList.size() - 1).getDate();
  1372.             final ImmutableTimeStampedCache<MeanPole> annualCache =
  1373.                     new ImmutableTimeStampedCache<MeanPole>(2, annualPoleList);

  1375.             // polynomial extension from IERS 2003, section 7.1.4, equations 23a and 23b
  1376.             final double xp0    = 0.054   * Constants.ARC_SECONDS_TO_RADIANS;
  1377.             final double xp0Dot = 0.00083 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
  1378.             final double yp0    = 0.357   * Constants.ARC_SECONDS_TO_RADIANS;
  1379.             final double yp0Dot = 0.00395 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;

  1381.             // constants from IERS 2003, section 6.2
  1382.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1383.             final double ratio        =  0.00115;

  1385.             return new TimeVectorFunction() {

  1387.                 /** {@inheritDoc} */
  1388.                 @Override
  1389.                 public double[] value(final AbsoluteDate date) {

  1391.                     // we can't compute anything before the range covered by the annual pole file
  1392.                     if (date.compareTo(firstAnnualPoleDate) <= 0) {
  1393.                         return new double[] {
  1394.                             0.0, 0.0
  1395.                         };
  1396.                     }

  1398.                     // evaluate mean pole
  1399.                     double meanPoleX = 0;
  1400.                     double meanPoleY = 0;
  1401.                     if (date.compareTo(lastAnnualPoleDate) <= 0) {
  1402.                         // we are within the range covered by the annual pole file,
  1403.                         // we interpolate within it
  1404.                         try {
  1405.                             final HermiteInterpolator interpolator = new HermiteInterpolator();
  1406.                             annualCache.getNeighbors(date).forEach(neighbor ->
  1407.                                 interpolator.addSamplePoint(neighbor.getDate().durationFrom(date),
  1408.                                                             new double[] {
  1409.                                                                 neighbor.getX(), neighbor.getY()
  1410.                                                             }));
  1411.                             final double[] interpolated = interpolator.value(0);
  1412.                             meanPoleX = interpolated[0];
  1413.                             meanPoleY = interpolated[1];
  1414.                         } catch (TimeStampedCacheException tsce) {
  1415.                             // this should never happen
  1416.                             throw new OrekitInternalError(tsce);
  1417.                         }
  1418.                     } else {

  1420.                         // we are after the range covered by the annual pole file,
  1421.                         // we use the polynomial extension
  1422.                         final double t = date.durationFrom(
  1423.                                 eopHistory.getTimeScales().getJ2000Epoch());
  1424.                         meanPoleX = xp0 + t * xp0Dot;
  1425.                         meanPoleY = yp0 + t * yp0Dot;

  1427.                     }

  1429.                     // evaluate wobble variables
  1430.                     final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1431.                     final double m1 = correction.getXp() - meanPoleX;
  1432.                     final double m2 = meanPoleY - correction.getYp();

  1434.                     return new double[] {
  1435.                         // the following correspond to the equations published in IERS 2003 conventions,
  1436.                         // section 6.2 page 65. In the publication, the equations read:
  1437.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1438.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1439.                         // However, it seems there are sign errors in these equations, which have
  1440.                         // been fixed in IERS 2010 conventions, section 6.4 page 94. In these newer
  1441.                         // publication, the equations read:
  1442.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1443.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1444.                         // the newer equations seem more consistent with the premises as the
  1445.                         // deformation due to the centrifugal potential has the form:
  1446.                         // −Ω²r²/2 sin 2θ Re [k₂(m₁ − im₂) exp(iλ)] where k₂ is the complex
  1447.                         // number 0.3077 + 0.0036i, so the real part in the previous equation is:
  1448.                         // A[Re(k₂) m₁ + Im(k₂) m₂)] cos λ + A[Re(k₂) m₂ - Im(k₂) m₁] sin λ
  1449.                         // identifying this with ∆C₂₁ cos λ + ∆S₂₁ sin λ we get:
  1450.                         // ∆C₂₁ = A Re(k₂) [m₁ + Im(k₂)/Re(k₂) m₂)]
  1451.                         // ∆S₂₁ = A Re(k₂) [m₂ - Im(k₂)/Re(k₂) m₁)]
  1452.                         // and Im(k₂)/Re(k₂) is very close to +0.00115
  1453.                         // As the equation as written in the IERS 2003 conventions are used in
  1454.                         // legacy systems, we have reproduced this alleged error here (and fixed it in
  1455.                         // the IERS 2010 conventions below) for validation purposes. We don't recommend
  1456.                         // using the IERS 2003 conventions for solid pole tide computation other than
  1457.                         // for validation or reproducibility of legacy applications behavior.
  1458.                         // As solid pole tide is small and as the sign change is on the smallest coefficient,
  1459.                         // the effect is quite small. A test case on a propagated orbit showed a position change
  1460.                         // slightly below 0.4m after a 30 days propagation on a Low Earth Orbit
  1461.                         globalFactor * (m1 - ratio * m2),
  1462.                         globalFactor * (m2 + ratio * m1),
  1463.                     };

  1465.                 }

  1467.                 /** {@inheritDoc} */
  1468.                 @Override
  1469.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

  1478.                     // evaluate mean pole
  1479.                     T meanPoleX = date.getField().getZero();
  1480.                     T meanPoleY = date.getField().getZero();
  1481.                     if (aDate.compareTo(lastAnnualPoleDate) <= 0) {
  1482.                         // we are within the range covered by the annual pole file,
  1483.                         // we interpolate within it
  1484.                         try {
  1485.                             final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
  1486.                             final T[] y = MathArrays.buildArray(date.getField(), 2);
  1487.                             final T zero = date.getField().getZero();
  1488.                             final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
  1489.                                                                                                        // for example removing derivatives
  1490.                                                                                                        // if T was DerivativeStructure
  1491.                             annualCache.getNeighbors(aDate).forEach(neighbor -> {
  1492.                                 y[0] = zero.add(neighbor.getX());
  1493.                                 y[1] = zero.add(neighbor.getY());
  1494.                                 interpolator.addSamplePoint(central.durationFrom(neighbor.getDate()).negate(), y);
  1495.                             });
  1496.                             final T[] interpolated = interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
  1497.                             meanPoleX = interpolated[0];
  1498.                             meanPoleY = interpolated[1];
  1499.                         } catch (TimeStampedCacheException tsce) {
  1500.                             // this should never happen
  1501.                             throw new OrekitInternalError(tsce);
  1502.                         }
  1503.                     } else {

  1505.                         // we are after the range covered by the annual pole file,
  1506.                         // we use the polynomial extension
  1507.                         final T t = date.durationFrom(
  1508.                                 eopHistory.getTimeScales().getJ2000Epoch());
  1509.                         meanPoleX = t.multiply(xp0Dot).add(xp0);
  1510.                         meanPoleY = t.multiply(yp0Dot).add(yp0);

  1512.                     }

  1514.                     // evaluate wobble variables
  1515.                     final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1516.                     final T m1 = correction.getXp().subtract(meanPoleX);
  1517.                     final T m2 = meanPoleY.subtract(correction.getYp());

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

  1521.                     // the following correspond to the equations published in IERS 2003 conventions,
  1522.                     // section 6.2 page 65. In the publication, the equations read:
  1523.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1524.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1525.                     // However, it seems there are sign errors in these equations, which have
  1526.                     // been fixed in IERS 2010 conventions, section 6.4 page 94. In these newer
  1527.                     // publication, the equations read:
  1528.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1529.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1530.                     // the newer equations seem more consistent with the premises as the
  1531.                     // deformation due to the centrifugal potential has the form:
  1532.                     // −Ω²r²/2 sin 2θ Re [k₂(m₁ − im₂) exp(iλ)] where k₂ is the complex
  1533.                     // number 0.3077 + 0.0036i, so the real part in the previous equation is:
  1534.                     // A[Re(k₂) m₁ + Im(k₂) m₂)] cos λ + A[Re(k₂) m₂ - Im(k₂) m₁] sin λ
  1535.                     // identifying this with ∆C₂₁ cos λ + ∆S₂₁ sin λ we get:
  1536.                     // ∆C₂₁ = A Re(k₂) [m₁ + Im(k₂)/Re(k₂) m₂)]
  1537.                     // ∆S₂₁ = A Re(k₂) [m₂ - Im(k₂)/Re(k₂) m₁)]
  1538.                     // and Im(k₂)/Re(k₂) is very close to +0.00115
  1539.                     // As the equation as written in the IERS 2003 conventions are used in
  1540.                     // legacy systems, we have reproduced this alleged error here (and fixed it in
  1541.                     // the IERS 2010 conventions below) for validation purposes. We don't recommend
  1542.                     // using the IERS 2003 conventions for solid pole tide computation other than
  1543.                     // for validation or reproducibility of legacy applications behavior.
  1544.                     // As solid pole tide is small and as the sign change is on the smallest coefficient,
  1545.                     // the effect is quite small. A test case on a propagated orbit showed a position change
  1546.                     // slightly below 0.4m after a 30 days propagation on a Low Earth Orbit
  1547.                     a[0] = m1.add(m2.multiply(-ratio)).multiply(globalFactor);
  1548.                     a[1] = m2.add(m1.multiply( ratio)).multiply(globalFactor);

  1550.                     return a;

  1552.                 }

  1554.             };

  1556.         }

  1558.         /** {@inheritDoc} */
  1559.         @Override
  1560.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {

  1562.             return new TimeVectorFunction() {

  1564.                 /** {@inheritDoc} */
  1565.                 @Override
  1566.                 public double[] value(final AbsoluteDate date) {
  1567.                     // there are no model for ocean pole tide prior to conventions 2010
  1568.                     return new double[] {
  1569.                         0.0, 0.0
  1570.                     };
  1571.                 }

  1573.                 /** {@inheritDoc} */
  1574.                 @Override
  1575.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1576.                     // there are no model for ocean pole tide prior to conventions 2010
  1577.                     return MathArrays.buildArray(date.getField(), 2);
  1578.                 }

  1580.             };
  1581.         }

  1583.         /** {@inheritDoc} */
  1584.         @Override
  1585.         public double[] getNominalTidalDisplacement() {
  1586.             return new double[] {
  1587.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  1588.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  1589.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  1590.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  1591.                 // H₀
  1592.                 -0.31460
  1593.             };
  1594.         }

  1596.         /** {@inheritDoc} */
  1597.         @Override
  1598.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
  1599.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  1600.                                                                   18, 15, 16, 17, 18);
  1601.         }

  1603.         /** {@inheritDoc} */
  1604.         @Override
  1605.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  1606.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  1607.                                                                 18, 15, 16, 17, 18);
  1608.         }

  1610.     },

  1612.     /** Constant for IERS 2010 conventions. */
  1613.     IERS_2010 {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  1660.         /** {@inheritDoc} */
  1661.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
  1662.                                                                  final TimeScales timeScales) {
  1663.             try (InputStream in = getStream(NUTATION_ARGUMENTS)) {
  1664.                 return new FundamentalNutationArguments(this, timeScale,
  1665.                         in, NUTATION_ARGUMENTS, timeScales);
  1666.             } catch (IOException e) {
  1667.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1668.             }
  1669.         }

  1671.         /** {@inheritDoc} */
  1672.         @Override
  1673.         public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {

  1675.             // epsilon 0 value from chapter 5, page 56, other terms from equation 5.40 page 65
  1676.             final PolynomialNutation epsilonA =
  1677.                     new PolynomialNutation(84381.406        * Constants.ARC_SECONDS_TO_RADIANS,
  1678.                                              -46.836769     * Constants.ARC_SECONDS_TO_RADIANS,
  1679.                                               -0.0001831    * Constants.ARC_SECONDS_TO_RADIANS,
  1680.                                                0.00200340   * Constants.ARC_SECONDS_TO_RADIANS,
  1681.                                               -0.000000576  * Constants.ARC_SECONDS_TO_RADIANS,
  1682.                                               -0.0000000434 * Constants.ARC_SECONDS_TO_RADIANS);

  1684.             return new TimeScalarFunction() {

  1686.                 /** {@inheritDoc} */
  1687.                 @Override
  1688.                 public double value(final AbsoluteDate date) {
  1689.                     return epsilonA.value(evaluateTC(date, timeScales));
  1690.                 }

  1692.                 /** {@inheritDoc} */
  1693.                 @Override
  1694.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1695.                     return epsilonA.value(evaluateTC(date, timeScales));
  1696.                 }

  1698.             };

  1700.         }

  1702.         /** {@inheritDoc} */
  1703.         @Override
  1704.         public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {

  1706.             // set up nutation arguments
  1707.             final FundamentalNutationArguments arguments =
  1708.                     getNutationArguments(timeScales);

  1710.             // set up Poisson series
  1711.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1712.             final PoissonSeriesParser parser =
  1713.                     new PoissonSeriesParser(17).
  1714.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  1715.                         withFirstDelaunay(4).
  1716.                         withFirstPlanetary(9).
  1717.                         withSinCos(0, 2, microAS, 3, microAS);
  1718.             final PoissonSeries.CompiledSeries xys;
  1719.             try (InputStream xIn = getStream(X_SERIES);
  1720.                  InputStream yIn = getStream(Y_SERIES);
  1721.                  InputStream sIn = getStream(S_SERIES)) {
  1722.                 final PoissonSeries xSeries = parser.parse(xIn, X_SERIES);
  1723.                 final PoissonSeries ySeries = parser.parse(yIn, Y_SERIES);
  1724.                 final PoissonSeries sSeries = parser.parse(sIn, S_SERIES);
  1725.                 xys = PoissonSeries.compile(xSeries, ySeries, sSeries);
  1726.             } catch (IOException e) {
  1727.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1728.             }

  1730.             // create a function evaluating the series
  1731.             return new TimeVectorFunction() {

  1733.                 /** {@inheritDoc} */
  1734.                 @Override
  1735.                 public double[] value(final AbsoluteDate date) {
  1736.                     return xys.value(arguments.evaluateAll(date));
  1737.                 }

  1739.                 /** {@inheritDoc} */
  1740.                 @Override
  1741.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1742.                     return xys.value(arguments.evaluateAll(date));
  1743.                 }

  1745.             };

  1747.         }

  1749.         /** {@inheritDoc} */
  1750.         public LoveNumbers getLoveNumbers() {
  1751.             return loadLoveNumbers(LOVE_NUMBERS);
  1752.         }

  1754.         /** {@inheritDoc} */
  1755.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  1756.                                                                      final TimeScales timeScales) {

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

  1761.             // set up Poisson series
  1762.             final PoissonSeriesParser k20Parser =
  1763.                     new PoissonSeriesParser(18).
  1764.                         withOptionalColumn(1).
  1765.                         withDoodson(4, 2).
  1766.                         withFirstDelaunay(10);
  1767.             final PoissonSeriesParser k21Parser =
  1768.                     new PoissonSeriesParser(18).
  1769.                         withOptionalColumn(1).
  1770.                         withDoodson(4, 3).
  1771.                         withFirstDelaunay(10);
  1772.             final PoissonSeriesParser k22Parser =
  1773.                     new PoissonSeriesParser(16).
  1774.                         withOptionalColumn(1).
  1775.                         withDoodson(4, 2).
  1776.                         withFirstDelaunay(10);

  1778.             final double pico = 1.0e-12;
  1779.             try (InputStream k0In = getStream(K20_FREQUENCY_DEPENDENCE);
  1780.                  InputStream k21In1 = getStream(K21_FREQUENCY_DEPENDENCE);
  1781.                  InputStream k21In2 = getStream(K21_FREQUENCY_DEPENDENCE);
  1782.                  InputStream k22In1 = getStream(K22_FREQUENCY_DEPENDENCE);
  1783.                  InputStream k22In2 = getStream(K22_FREQUENCY_DEPENDENCE)) {
  1784.                 final PoissonSeries c20Series =
  1785.                         k20Parser.
  1786.                         withSinCos(0, 18, -pico, 16, pico).
  1787.                         parse(k0In, K20_FREQUENCY_DEPENDENCE);
  1788.                 final PoissonSeries c21Series =
  1789.                         k21Parser.
  1790.                         withSinCos(0, 17, pico, 18, pico).
  1791.                         parse(k21In1, K21_FREQUENCY_DEPENDENCE);
  1792.                 final PoissonSeries s21Series =
  1793.                         k21Parser.
  1794.                         withSinCos(0, 18, -pico, 17, pico).
  1795.                         parse(k21In2, K21_FREQUENCY_DEPENDENCE);
  1796.                 final PoissonSeries c22Series =
  1797.                         k22Parser.
  1798.                         withSinCos(0, -1, pico, 16, pico).
  1799.                         parse(k22In1, K22_FREQUENCY_DEPENDENCE);
  1800.                 final PoissonSeries s22Series =
  1801.                         k22Parser.
  1802.                         withSinCos(0, 16, -pico, -1, pico).
  1803.                         parse(k22In2, K22_FREQUENCY_DEPENDENCE);

  1805.                 return new TideFrequencyDependenceFunction(arguments,
  1806.                                                            c20Series,
  1807.                                                            c21Series, s21Series,
  1808.                                                            c22Series, s22Series);
  1809.             } catch (IOException e) {
  1810.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1811.             }

  1813.         }

  1815.         /** {@inheritDoc} */
  1816.         @Override
  1817.         public double getPermanentTide() {
  1818.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1819.         }

  1821.         /** Compute pole wobble variables m₁ and m₂.
  1822.          * @param date current date
  1823.          * @param eopHistory EOP history
  1824.          * @return array containing m₁ and m₂
  1825.          */
  1826.         private double[] computePoleWobble(final AbsoluteDate date, final EOPHistory eopHistory) {

  1828.             // polynomial model from IERS 2010, table 7.7
  1829.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1830.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1831.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1832.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1833.             final AbsoluteDate changeDate =
  1834.                     new AbsoluteDate(2010, 1, 1, eopHistory.getTimeScales().getTT());

  1836.             // evaluate mean pole
  1837.             final double[] xPolynomial;
  1838.             final double[] yPolynomial;
  1839.             if (date.compareTo(changeDate) <= 0) {
  1840.                 xPolynomial = new double[] {
  1841.                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
  1842.                 };
  1843.                 yPolynomial = new double[] {
  1844.                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
  1845.                 };
  1846.             } else {
  1847.                 xPolynomial = new double[] {
  1848.                     23.513 * f0, 7.6141 * f1
  1849.                 };
  1850.                 yPolynomial = new double[] {
  1851.                     358.891 * f0,  -0.6287 * f1
  1852.                 };
  1853.             }
  1854.             double meanPoleX = 0;
  1855.             double meanPoleY = 0;
  1856.             final double t = date.durationFrom(
  1857.                     eopHistory.getTimeScales().getJ2000Epoch());
  1858.             for (int i = xPolynomial.length - 1; i >= 0; --i) {
  1859.                 meanPoleX = meanPoleX * t + xPolynomial[i];
  1860.             }
  1861.             for (int i = yPolynomial.length - 1; i >= 0; --i) {
  1862.                 meanPoleY = meanPoleY * t + yPolynomial[i];
  1863.             }

  1865.             // evaluate wobble variables
  1866.             final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1867.             final double m1 = correction.getXp() - meanPoleX;
  1868.             final double m2 = meanPoleY - correction.getYp();

  1870.             return new double[] {
  1871.                 m1, m2
  1872.             };

  1874.         }

  1876.         /** Compute pole wobble variables m₁ and m₂.
  1877.          * @param date current date
  1878.          * @param <T> type of the field elements
  1879.          * @param eopHistory EOP history
  1880.          * @return array containing m₁ and m₂
  1881.          */
  1882.         private <T extends RealFieldElement<T>> T[] computePoleWobble(final FieldAbsoluteDate<T> date, final EOPHistory eopHistory) {

  1884.             // polynomial model from IERS 2010, table 7.7
  1885.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1886.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1887.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1888.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1889.             final AbsoluteDate changeDate =
  1890.                     new AbsoluteDate(2010, 1, 1, eopHistory.getTimeScales().getTT());

  1892.             // evaluate mean pole
  1893.             final double[] xPolynomial;
  1894.             final double[] yPolynomial;
  1895.             if (date.toAbsoluteDate().compareTo(changeDate) <= 0) {
  1896.                 xPolynomial = new double[] {
  1897.                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
  1898.                 };
  1899.                 yPolynomial = new double[] {
  1900.                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
  1901.                 };
  1902.             } else {
  1903.                 xPolynomial = new double[] {
  1904.                     23.513 * f0, 7.6141 * f1
  1905.                 };
  1906.                 yPolynomial = new double[] {
  1907.                     358.891 * f0,  -0.6287 * f1
  1908.                 };
  1909.             }
  1910.             T meanPoleX = date.getField().getZero();
  1911.             T meanPoleY = date.getField().getZero();
  1912.             final T t = date.durationFrom(
  1913.                     eopHistory.getTimeScales().getJ2000Epoch());
  1914.             for (int i = xPolynomial.length - 1; i >= 0; --i) {
  1915.                 meanPoleX = meanPoleX.multiply(t).add(xPolynomial[i]);
  1916.             }
  1917.             for (int i = yPolynomial.length - 1; i >= 0; --i) {
  1918.                 meanPoleY = meanPoleY.multiply(t).add(yPolynomial[i]);
  1919.             }

  1921.             // evaluate wobble variables
  1922.             final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1923.             final T[] m = MathArrays.buildArray(date.getField(), 2);
  1924.             m[0] = correction.getXp().subtract(meanPoleX);
  1925.             m[1] = meanPoleY.subtract(correction.getYp());

  1927.             return m;

  1929.         }

  1931.         /** {@inheritDoc} */
  1932.         @Override
  1933.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  1935.             // constants from IERS 2010, section 6.4
  1936.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1937.             final double ratio        =  0.00115;

  1939.             return new TimeVectorFunction() {

  1941.                 /** {@inheritDoc} */
  1942.                 @Override
  1943.                 public double[] value(final AbsoluteDate date) {

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

  1948.                     return new double[] {
  1949.                         // the following correspond to the equations published in IERS 2010 conventions,
  1950.                         // section 6.4 page 94. The equations read:
  1951.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1952.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1953.                         // These equations seem to fix what was probably a sign error in IERS 2003
  1954.                         // conventions section 6.2 page 65. In this older publication, the equations read:
  1955.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1956.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1957.                         globalFactor * (wobbleM[0] + ratio * wobbleM[1]),
  1958.                         globalFactor * (wobbleM[1] - ratio * wobbleM[0])
  1959.                     };

  1961.                 }

  1963.                 /** {@inheritDoc} */
  1964.                 @Override
  1965.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

  1972.                     // the following correspond to the equations published in IERS 2010 conventions,
  1973.                     // section 6.4 page 94. The equations read:
  1974.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1975.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1976.                     // These equations seem to fix what was probably a sign error in IERS 2003
  1977.                     // conventions section 6.2 page 65. In this older publication, the equations read:
  1978.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1979.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1980.                     a[0] = wobbleM[0].add(wobbleM[1].multiply( ratio)).multiply(globalFactor);
  1981.                     a[1] = wobbleM[1].add(wobbleM[0].multiply(-ratio)).multiply(globalFactor);

  1983.                     return a;

  1985.                 }

  1987.             };

  1989.         }

  1991.         /** {@inheritDoc} */
  1992.         @Override
  1993.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {

  1995.             return new TimeVectorFunction() {

  1997.                 /** {@inheritDoc} */
  1998.                 @Override
  1999.                 public double[] value(final AbsoluteDate date) {

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

  2004.                     return new double[] {
  2005.                         // the following correspond to the equations published in IERS 2010 conventions,
  2006.                         // section 6.4 page 94 equation 6.24:
  2007.                         // ∆C₂₁ = −2.1778 × 10⁻¹⁰ (m₁ − 0.01724m₂)
  2008.                         // ∆S₂₁ = −1.7232 × 10⁻¹⁰ (m₂ − 0.03365m₁)
  2009.                         -2.1778e-10 * (wobbleM[0] - 0.01724 * wobbleM[1]) / Constants.ARC_SECONDS_TO_RADIANS,
  2010.                         -1.7232e-10 * (wobbleM[1] - 0.03365 * wobbleM[0]) / Constants.ARC_SECONDS_TO_RADIANS
  2011.                     };

  2013.                 }

  2015.                 /** {@inheritDoc} */
  2016.                 @Override
  2017.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

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

  2031.                     return v;

  2033.                 }

  2035.             };

  2037.         }

  2039.         /** {@inheritDoc} */
  2040.         @Override
  2041.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

  2043.             // set up the conventional polynomials
  2044.             // the following values are from equation 5.40 in IERS 2010 conventions
  2045.             final PolynomialNutation psiA =
  2046.                     new PolynomialNutation(   0.0,
  2047.                                            5038.481507     * Constants.ARC_SECONDS_TO_RADIANS,
  2048.                                              -1.0790069    * Constants.ARC_SECONDS_TO_RADIANS,
  2049.                                              -0.00114045   * Constants.ARC_SECONDS_TO_RADIANS,
  2050.                                               0.000132851  * Constants.ARC_SECONDS_TO_RADIANS,
  2051.                                              -0.0000000951 * Constants.ARC_SECONDS_TO_RADIANS);
  2052.             final PolynomialNutation omegaA =
  2053.                     new PolynomialNutation(getMeanObliquityFunction(timeScales)
  2054.                             .value(getNutationReferenceEpoch(timeScales)),
  2055.                                            -0.025754     * Constants.ARC_SECONDS_TO_RADIANS,
  2056.                                             0.0512623    * Constants.ARC_SECONDS_TO_RADIANS,
  2057.                                            -0.00772503   * Constants.ARC_SECONDS_TO_RADIANS,
  2058.                                            -0.000000467  * Constants.ARC_SECONDS_TO_RADIANS,
  2059.                                             0.0000003337 * Constants.ARC_SECONDS_TO_RADIANS);
  2060.             final PolynomialNutation chiA =
  2061.                     new PolynomialNutation( 0.0,
  2062.                                            10.556403     * Constants.ARC_SECONDS_TO_RADIANS,
  2063.                                            -2.3814292    * Constants.ARC_SECONDS_TO_RADIANS,
  2064.                                            -0.00121197   * Constants.ARC_SECONDS_TO_RADIANS,
  2065.                                             0.000170663  * Constants.ARC_SECONDS_TO_RADIANS,
  2066.                                            -0.0000000560 * Constants.ARC_SECONDS_TO_RADIANS);

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

  2070.         }

  2072.         /** {@inheritDoc} */
  2073.         @Override
  2074.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  2076.             // set up nutation arguments
  2077.             final FundamentalNutationArguments arguments =
  2078.                     getNutationArguments(timeScales);

  2080.             // set up Poisson series
  2081.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2082.             final PoissonSeriesParser parser =
  2083.                     new PoissonSeriesParser(17).
  2084.                         withFirstDelaunay(4).
  2085.                         withFirstPlanetary(9).
  2086.                         withSinCos(0, 2, microAS, 3, microAS);
  2087.             final PoissonSeries.CompiledSeries psiEpsilonSeries;
  2088.             try (InputStream psiIn = getStream(PSI_SERIES);
  2089.                  InputStream epsilonIn = getStream(EPSILON_SERIES)) {
  2090.                 final PoissonSeries psiSeries     = parser.parse(psiIn, PSI_SERIES);
  2091.                 final PoissonSeries epsilonSeries = parser.parse(epsilonIn, EPSILON_SERIES);
  2092.                 psiEpsilonSeries = PoissonSeries.compile(psiSeries, epsilonSeries);
  2093.             } catch (IOException e) {
  2094.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2095.             }

  2097.             return new TimeVectorFunction() {

  2099.                 /** {@inheritDoc} */
  2100.                 @Override
  2101.                 public double[] value(final AbsoluteDate date) {
  2102.                     final BodiesElements elements = arguments.evaluateAll(date);
  2103.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  2104.                     return new double[] {
  2105.                         psiEpsilon[0], psiEpsilon[1],
  2106.                         IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  2107.                     };
  2108.                 }

  2110.                 /** {@inheritDoc} */
  2111.                 @Override
  2112.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  2113.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  2114.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  2115.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  2116.                     result[0] = psiEpsilon[0];
  2117.                     result[1] = psiEpsilon[1];
  2118.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  2119.                     return result;
  2120.                 }

  2122.             };

  2124.         }

  2126.         /** {@inheritDoc} */
  2127.         @Override
  2128.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  2129.                                                   final TimeScales timeScales) {

  2131.             // Earth Rotation Angle
  2132.             final StellarAngleCapitaine era =
  2133.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  2135.             // Polynomial part of the apparent sidereal time series
  2136.             // which is the opposite of Equation of Origins (EO)
  2137.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2138.             final PoissonSeriesParser parser =
  2139.                     new PoissonSeriesParser(17).
  2140.                         withFirstDelaunay(4).
  2141.                         withFirstPlanetary(9).
  2142.                         withSinCos(0, 2, microAS, 3, microAS).
  2143.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  2144.             final PolynomialNutation minusEO;
  2145.             try (InputStream in = getStream(GST_SERIES)) {
  2146.                 minusEO = parser.parse(in, GST_SERIES).getPolynomial();
  2147.             } catch (IOException e) {
  2148.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2149.             }

  2151.             // create a function evaluating the series
  2152.             return new TimeScalarFunction() {

  2154.                 /** {@inheritDoc} */
  2155.                 @Override
  2156.                 public double value(final AbsoluteDate date) {
  2157.                     return era.value(date) + minusEO.value(evaluateTC(date, timeScales));
  2158.                 }

  2160.                 /** {@inheritDoc} */
  2161.                 @Override
  2162.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  2163.                     return era.value(date).add(minusEO.value(evaluateTC(date, timeScales)));
  2164.                 }

  2166.             };

  2168.         }

  2170.         /** {@inheritDoc} */
  2171.         @Override
  2172.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  2173.                                                       final TimeScales timeScales) {

  2175.             // Earth Rotation Angle
  2176.             final StellarAngleCapitaine era =
  2177.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  2179.             // Polynomial part of the apparent sidereal time series
  2180.             // which is the opposite of Equation of Origins (EO)
  2181.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2182.             final PoissonSeriesParser parser =
  2183.                     new PoissonSeriesParser(17).
  2184.                         withFirstDelaunay(4).
  2185.                         withFirstPlanetary(9).
  2186.                         withSinCos(0, 2, microAS, 3, microAS).
  2187.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  2188.             final PolynomialNutation minusEO;
  2189.             try (InputStream in = getStream(GST_SERIES)) {
  2190.                 minusEO = parser.parse(in, GST_SERIES).getPolynomial();
  2191.             } catch (IOException e) {
  2192.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2193.             }

  2195.             // create a function evaluating the series
  2196.             return new TimeScalarFunction() {

  2198.                 /** {@inheritDoc} */
  2199.                 @Override
  2200.                 public double value(final AbsoluteDate date) {
  2201.                     return era.getRate() + minusEO.derivative(evaluateTC(date, timeScales));
  2202.                 }

  2204.                 /** {@inheritDoc} */
  2205.                 @Override
  2206.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  2207.                     return minusEO.derivative(evaluateTC(date, timeScales)).add(era.getRate());
  2208.                 }

  2210.             };

  2212.         }

  2214.         /** {@inheritDoc} */
  2215.         @Override
  2216.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  2217.                                                   final EOPHistory eopHistory,
  2218.                                                   final TimeScales timeScales) {

  2220.             // set up nutation arguments
  2221.             final FundamentalNutationArguments arguments =
  2222.                     getNutationArguments(timeScales);

  2224.             // mean obliquity function
  2225.             final TimeScalarFunction epsilon = getMeanObliquityFunction(timeScales);

  2227.             // set up Poisson series
  2228.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2229.             final PoissonSeriesParser baseParser =
  2230.                     new PoissonSeriesParser(17).
  2231.                         withFirstDelaunay(4).
  2232.                         withFirstPlanetary(9).
  2233.                         withSinCos(0, 2, microAS, 3, microAS);
  2234.             final PoissonSeriesParser gstParser  = baseParser.withPolynomialPart('t', Unit.ARC_SECONDS);
  2235.             final PoissonSeries.CompiledSeries psiGstSeries;
  2236.             try (InputStream psiIn = getStream(PSI_SERIES);
  2237.                  InputStream gstIn = getStream(GST_SERIES)) {
  2238.                 final PoissonSeries psiSeries        = baseParser.parse(psiIn, PSI_SERIES);
  2239.                 final PoissonSeries gstSeries        = gstParser.parse(gstIn, GST_SERIES);
  2240.                 psiGstSeries = PoissonSeries.compile(psiSeries, gstSeries);
  2241.             } catch (IOException e) {
  2242.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2243.             }

  2245.             // ERA function
  2246.             final TimeScalarFunction era =
  2247.                     getEarthOrientationAngleFunction(ut1, timeScales.getTAI());

  2249.             return new TimeScalarFunction() {

  2251.                 /** {@inheritDoc} */
  2252.                 @Override
  2253.                 public double value(final AbsoluteDate date) {

  2255.                     // evaluate equation of origins
  2256.                     final BodiesElements elements = arguments.evaluateAll(date);
  2257.                     final double[] angles = psiGstSeries.value(elements);
  2258.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  2259.                     final double deltaPsi = angles[0] + ddPsi;
  2260.                     final double epsilonA = epsilon.value(date);

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

  2266.                 }

  2268.                 /** {@inheritDoc} */
  2269.                 @Override
  2270.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  2272.                     // evaluate equation of origins
  2273.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  2274.                     final T[] angles = psiGstSeries.value(elements);
  2275.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  2276.                     final T deltaPsi = angles[0].add(ddPsi);
  2277.                     final T epsilonA = epsilon.value(date);

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

  2283.                 }

  2285.             };

  2287.         }

  2289.         /** {@inheritDoc} */
  2290.         @Override
  2291.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

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

  2301.             // set up Poisson series
  2302.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2303.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
  2304.                     withOptionalColumn(1).
  2305.                     withGamma(2).
  2306.                     withFirstDelaunay(3);
  2307.             try (InputStream xpIn = getStream(TIDAL_CORRECTION_XP_YP_SERIES);
  2308.                  InputStream ypIn = getStream(TIDAL_CORRECTION_XP_YP_SERIES);
  2309.                  InputStream ut1In = getStream(TIDAL_CORRECTION_UT1_SERIES)) {
  2310.                 final PoissonSeries xSeries =
  2311.                         xyParser.
  2312.                         withSinCos(0, 10, microAS, 11, microAS).
  2313.                         parse(xpIn, TIDAL_CORRECTION_XP_YP_SERIES);
  2314.                 final PoissonSeries ySeries =
  2315.                         xyParser.
  2316.                         withSinCos(0, 12, microAS, 13, microAS).
  2317.                         parse(ypIn, TIDAL_CORRECTION_XP_YP_SERIES);

  2319.                 final double microS = 1.0e-6;
  2320.                 final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  2321.                         withOptionalColumn(1).
  2322.                         withGamma(2).
  2323.                         withFirstDelaunay(3).
  2324.                         withSinCos(0, 10, microS, 11, microS);
  2325.                 final PoissonSeries ut1Series =
  2326.                         ut1Parser.parse(ut1In, TIDAL_CORRECTION_UT1_SERIES);

  2328.                 return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
  2329.             } catch (IOException e) {
  2330.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2331.             }

  2333.         }

  2335.         /** {@inheritDoc} */
  2336.         @Override
  2337.         public double[] getNominalTidalDisplacement() {
  2338.             return new double[] {
  2339.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  2340.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  2341.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  2342.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  2343.                 // H₀
  2344.                 -0.31460
  2345.             };
  2346.         }

  2348.         /** {@inheritDoc} */
  2349.         @Override
  2350.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
  2351.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  2352.                                                                   18, 15, 16, 17, 18);
  2353.         }

  2355.         /** {@inheritDoc} */
  2356.         @Override
  2357.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  2358.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  2359.                                                                 18, 15, 16, 17, 18);
  2360.         }

  2362.     };

  2364.     /** Pattern for delimiting regular expressions. */
  2365.     private static final Pattern SEPARATOR = Pattern.compile("\\p{Space}+");

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

  2370.     /** Get the reference epoch for fundamental nutation arguments.
  2371.      *
  2372.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2373.      *
  2374.      * @return reference epoch for fundamental nutation arguments
  2375.      * @since 6.1
  2376.      * @see #getNutationReferenceEpoch(TimeScales)
  2377.      */
  2378.     @DefaultDataContext
  2379.     public AbsoluteDate getNutationReferenceEpoch() {
  2380.         return getNutationReferenceEpoch(DataContext.getDefault().getTimeScales());
  2381.     }

  2383.     /**
  2384.      * Get the reference epoch for fundamental nutation arguments.
  2385.      *
  2386.      * @param timeScales to use for the reference epoch.
  2387.      * @return reference epoch for fundamental nutation arguments
  2388.      * @since 10.1
  2389.      */
  2390.     public AbsoluteDate getNutationReferenceEpoch(final TimeScales timeScales) {
  2391.         // IERS 1996, IERS 2003 and IERS 2010 use the same J2000.0 reference date
  2392.         return timeScales.getJ2000Epoch();
  2393.     }

  2395.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2396.      *
  2397.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2398.      *
  2399.      * @param date current date
  2400.      * @return date offset in Julian centuries
  2401.      * @since 6.1
  2402.      * @see #evaluateTC(AbsoluteDate, TimeScales)
  2403.      */
  2404.     @DefaultDataContext
  2405.     public double evaluateTC(final AbsoluteDate date) {
  2406.         return evaluateTC(date, DataContext.getDefault().getTimeScales());
  2407.     }

  2409.     /**
  2410.      * Evaluate the date offset between the current date and the {@link
  2411.      * #getNutationReferenceEpoch() reference date}.
  2412.      *
  2413.      * @param date       current date
  2414.      * @param timeScales used in the evaluation.
  2415.      * @return date offset in Julian centuries
  2416.      * @since 10.1
  2417.      */
  2418.     public double evaluateTC(final AbsoluteDate date, final TimeScales timeScales) {
  2419.         return date.durationFrom(getNutationReferenceEpoch(timeScales)) /
  2420.                 Constants.JULIAN_CENTURY;
  2421.     }

  2423.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2424.      *
  2425.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2426.      *
  2427.      * @param date current date
  2428.      * @param <T> type of the field elements
  2429.      * @return date offset in Julian centuries
  2430.      * @since 9.0
  2431.      * @see #evaluateTC(FieldAbsoluteDate, TimeScales)
  2432.      */
  2433.     @DefaultDataContext
  2434.     public <T extends RealFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date) {
  2435.         return evaluateTC(date, DataContext.getDefault().getTimeScales());
  2436.     }

  2438.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2439.      * @param <T> type of the field elements
  2440.      * @param date current date
  2441.      * @param timeScales used in the evaluation.
  2442.      * @return date offset in Julian centuries
  2443.      * @since 10.1
  2444.      */
  2445.     public <T extends RealFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date,
  2446.                                                         final TimeScales timeScales) {
  2447.         return date.durationFrom(getNutationReferenceEpoch(timeScales))
  2448.                 .divide(Constants.JULIAN_CENTURY);
  2449.     }

  2451.     /**
  2452.      * Get the fundamental nutation arguments. Does not compute GMST based values: gamma,
  2453.      * gammaDot.
  2454.      *
  2455.      * @param timeScales other time scales used in the computation including TAI and TT.
  2456.      * @return fundamental nutation arguments
  2457.      * @see #getNutationArguments(TimeScale, TimeScales)
  2458.      * @since 10.1
  2459.      */
  2460.     protected FundamentalNutationArguments getNutationArguments(
  2461.             final TimeScales timeScales) {

  2463.         return getNutationArguments(null, timeScales);
  2464.     }

  2466.     /** Get the fundamental nutation arguments.
  2467.      *
  2468.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2469.      *
  2470.      * @param timeScale time scale for computing Greenwich Mean Sidereal Time
  2471.      * (typically {@link TimeScales#getUT1(IERSConventions, boolean) UT1})
  2472.      * @return fundamental nutation arguments
  2473.      * @since 6.1
  2474.      * @see #getNutationArguments(TimeScale, TimeScales)
  2475.      * @see #getNutationArguments(TimeScales)
  2476.      */
  2477.     @DefaultDataContext
  2478.     public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale) {
  2479.         return getNutationArguments(timeScale, DataContext.getDefault().getTimeScales());
  2480.     }

  2482.     /**
  2483.      * Get the fundamental nutation arguments.
  2484.      *
  2485.      * @param timeScale time scale for computing Greenwich Mean Sidereal Time (typically
  2486.      *                  {@link TimeScales#getUT1(IERSConventions, boolean) UT1})
  2487.      * @param timeScales other time scales used in the computation including TAI and TT.
  2488.      * @return fundamental nutation arguments
  2489.      * @since 10.1
  2490.      */
  2491.     public abstract FundamentalNutationArguments getNutationArguments(
  2492.             TimeScale timeScale,
  2493.             TimeScales timeScales);

  2495.     /** Get the function computing mean obliquity of the ecliptic.
  2496.      *
  2497.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2498.      *
  2499.      * @return function computing mean obliquity of the ecliptic
  2500.      * @since 6.1
  2501.      * @see #getMeanObliquityFunction(TimeScales)
  2502.      */
  2503.     @DefaultDataContext
  2504.     public TimeScalarFunction getMeanObliquityFunction() {
  2505.         return getMeanObliquityFunction(DataContext.getDefault().getTimeScales());
  2506.     }

  2508.     /**
  2509.      * Get the function computing mean obliquity of the ecliptic.
  2510.      *
  2511.      * @param timeScales used in computing the function.
  2512.      * @return function computing mean obliquity of the ecliptic
  2513.      * @since 10.1
  2514.      */
  2515.     public abstract TimeScalarFunction getMeanObliquityFunction(TimeScales timeScales);

  2517.     /** Get the function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components.
  2518.      * <p>
  2519.      * The returned function computes the two X, Y components of CIP and the S+XY/2 component of the non-rotating CIO.
  2520.      * </p>
  2521.      *
  2522.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2523.      *
  2524.      * @return function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components
  2525.      * @since 6.1
  2526.      * @see #getXYSpXY2Function(TimeScales)
  2527.      */
  2528.     @DefaultDataContext
  2529.     public TimeVectorFunction getXYSpXY2Function() {
  2530.         return getXYSpXY2Function(DataContext.getDefault().getTimeScales());
  2531.     }

  2533.     /**
  2534.      * Get the function computing the Celestial Intermediate Pole and Celestial
  2535.      * Intermediate Origin components.
  2536.      * <p>
  2537.      * The returned function computes the two X, Y components of CIP and the S+XY/2
  2538.      * component of the non-rotating CIO.
  2539.      * </p>
  2540.      *
  2541.      * @param timeScales used to define the function.
  2542.      * @return function computing the Celestial Intermediate Pole and Celestial
  2543.      * Intermediate Origin components
  2544.      * @since 10.1
  2545.      */
  2546.     public abstract TimeVectorFunction getXYSpXY2Function(TimeScales timeScales);

  2548.     /** Get the function computing the raw Earth Orientation Angle.
  2549.      *
  2550.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2551.      *
  2552.      * <p>
  2553.      * The raw angle does not contain any correction. If for example dTU1 correction
  2554.      * due to tidal effect is desired, it must be added afterward by the caller.
  2555.      * The returned value contain the angle as the value and the angular rate as
  2556.      * the first derivative.
  2557.      * </p>
  2558.      * @param ut1 UT1 time scale
  2559.      * @return function computing the rawEarth Orientation Angle, in the non-rotating origin paradigm
  2560.      * @since 6.1
  2561.      * @see #getEarthOrientationAngleFunction(TimeScale, TimeScale)
  2562.      */
  2563.     @DefaultDataContext
  2564.     public TimeScalarFunction getEarthOrientationAngleFunction(final TimeScale ut1) {
  2565.         return getEarthOrientationAngleFunction(ut1,
  2566.                 DataContext.getDefault().getTimeScales().getTAI());
  2567.     }

  2569.     /** Get the function computing the raw Earth Orientation Angle.
  2570.      * <p>
  2571.      * The raw angle does not contain any correction. If for example dTU1 correction
  2572.      * due to tidal effect is desired, it must be added afterward by the caller.
  2573.      * The returned value contain the angle as the value and the angular rate as
  2574.      * the first derivative.
  2575.      * </p>
  2576.      * @param ut1 UT1 time scale
  2577.      * @param tai TAI time scale
  2578.      * @return function computing the rawEarth Orientation Angle, in the non-rotating origin paradigm
  2579.      * @since 10.1
  2580.      */
  2581.     public TimeScalarFunction getEarthOrientationAngleFunction(final TimeScale ut1,
  2582.                                                                final TimeScale tai) {
  2583.         return new StellarAngleCapitaine(ut1, tai);
  2584.     }

  2586.     /** Get the function computing the precession angles.
  2587.      * <p>
  2588.      * The function returned computes the three precession angles
  2589.      * ψ<sub>A</sub> (around Z axis), ω<sub>A</sub> (around X axis)
  2590.      * and χ<sub>A</sub> (around Z axis). The constant angle ε₀
  2591.      * for the fourth rotation (around X axis) can be retrieved by evaluating the
  2592.      * function returned by {@link #getMeanObliquityFunction()} at {@link
  2593.      * #getNutationReferenceEpoch() nutation reference epoch}.
  2594.      * </p>
  2595.      *
  2596.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2597.      *
  2598.      * @return function computing the precession angle
  2599.      * @since 6.1
  2600.      * @see #getPrecessionFunction(TimeScales)
  2601.      */
  2602.     @DefaultDataContext
  2603.     public TimeVectorFunction getPrecessionFunction()
  2604.     {
  2605.         return getPrecessionFunction(DataContext.getDefault().getTimeScales());
  2606.     }

  2608.     /** Get the function computing the precession angles.
  2609.      * <p>
  2610.      * The function returned computes the three precession angles
  2611.      * ψ<sub>A</sub> (around Z axis), ω<sub>A</sub> (around X axis)
  2612.      * and χ<sub>A</sub> (around Z axis). The constant angle ε₀
  2613.      * for the fourth rotation (around X axis) can be retrieved by evaluating the
  2614.      * function returned by {@link #getMeanObliquityFunction()} at {@link
  2615.      * #getNutationReferenceEpoch() nutation reference epoch}.
  2616.      * </p>
  2617.      * @return function computing the precession angle
  2618.      * @since 10.1
  2619.      * @param timeScales used to define the function.
  2620.      */
  2621.     public abstract TimeVectorFunction getPrecessionFunction(TimeScales timeScales);

  2623.     /** Get the function computing the nutation angles.
  2624.      *
  2625.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2626.      *
  2627.      * <p>
  2628.      * The function returned computes the two classical angles ΔΨ and Δε,
  2629.      * and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
  2630.      * resolution C7 (the correction is forced to 0 before this date)
  2631.      * </p>
  2632.      * @return function computing the nutation in longitude ΔΨ and Δε
  2633.      * and the correction of equation of equinoxes
  2634.      * @since 6.1
  2635.      */
  2636.     @DefaultDataContext
  2637.     public TimeVectorFunction getNutationFunction() {
  2638.         return getNutationFunction(DataContext.getDefault().getTimeScales());
  2639.     }

  2641.     /** Get the function computing the nutation angles.
  2642.      * <p>
  2643.      * The function returned computes the two classical angles ΔΨ and Δε,
  2644.      * and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
  2645.      * resolution C7 (the correction is forced to 0 before this date)
  2646.      * </p>
  2647.      * @return function computing the nutation in longitude ΔΨ and Δε
  2648.      * and the correction of equation of equinoxes
  2649.      * @param timeScales used in the computation including TAI and TT.
  2650.      * @since 10.1
  2651.      */
  2652.     public abstract TimeVectorFunction getNutationFunction(TimeScales timeScales);

  2654.     /** Get the function computing Greenwich mean sidereal time, in radians.
  2655.      *
  2656.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2657.      *
  2658.      * @param ut1 UT1 time scale
  2659.      * @return function computing Greenwich mean sidereal time
  2660.      * @since 6.1
  2661.      * @see #getGMSTFunction(TimeScale, TimeScales)
  2662.      */
  2663.     @DefaultDataContext
  2664.     public TimeScalarFunction getGMSTFunction(final TimeScale ut1) {
  2665.         return getGMSTFunction(ut1, DataContext.getDefault().getTimeScales());
  2666.     }

  2668.     /**
  2669.      * Get the function computing Greenwich mean sidereal time, in radians.
  2670.      *
  2671.      * @param ut1 UT1 time scale
  2672.      * @param timeScales other time scales used in the computation including TAI and TT.
  2673.      * @return function computing Greenwich mean sidereal time
  2674.      * @since 10.1
  2675.      */
  2676.     public abstract TimeScalarFunction getGMSTFunction(TimeScale ut1,
  2677.                                                        TimeScales timeScales);

  2679.     /** Get the function computing Greenwich mean sidereal time rate, in radians per second.
  2680.      *
  2681.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2682.      *
  2683.      * @param ut1 UT1 time scale
  2684.      * @return function computing Greenwich mean sidereal time rate
  2685.      * @since 9.0
  2686.      * @see #getGMSTRateFunction(TimeScale, TimeScales)
  2687.      */
  2688.     @DefaultDataContext
  2689.     public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1) {
  2690.         return getGMSTRateFunction(ut1,
  2691.                 DataContext.getDefault().getTimeScales());
  2692.     }

  2694.     /**
  2695.      * Get the function computing Greenwich mean sidereal time rate, in radians per
  2696.      * second.
  2697.      *
  2698.      * @param ut1 UT1 time scale
  2699.      * @param timeScales other time scales used in the computation including TAI and TT.
  2700.      * @return function computing Greenwich mean sidereal time rate
  2701.      * @since 10.1
  2702.      */
  2703.     public abstract TimeScalarFunction getGMSTRateFunction(TimeScale ut1,
  2704.                                                            TimeScales timeScales);

  2706.     /**
  2707.      * Get the function computing Greenwich apparent sidereal time, in radians.
  2708.      *
  2709.      * <p>This method uses the {@link DataContext#getDefault() default data context} if
  2710.      * {@code eopHistory == null}.
  2711.      *
  2712.      * @param ut1        UT1 time scale
  2713.      * @param eopHistory EOP history. If {@code null} then no nutation correction is
  2714.      *                   applied for EOP.
  2715.      * @return function computing Greenwich apparent sidereal time
  2716.      * @since 6.1
  2717.      * @see #getGASTFunction(TimeScale, EOPHistory, TimeScales)
  2718.      */
  2719.     @DefaultDataContext
  2720.     public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  2721.                                               final EOPHistory eopHistory) {
  2722.         final TimeScales timeScales = eopHistory != null ?
  2723.                 eopHistory.getTimeScales() :
  2724.                 DataContext.getDefault().getTimeScales();
  2725.         return getGASTFunction(ut1, eopHistory, timeScales);
  2726.     }

  2728.     /**
  2729.      * Get the function computing Greenwich apparent sidereal time, in radians.
  2730.      *
  2731.      * @param ut1        UT1 time scale
  2732.      * @param eopHistory EOP history. If {@code null} then no nutation correction is
  2733.      *                   applied for EOP.
  2734.      * @param timeScales        TAI time scale.
  2735.      * @return function computing Greenwich apparent sidereal time
  2736.      * @since 10.1
  2737.      */
  2738.     public abstract TimeScalarFunction getGASTFunction(TimeScale ut1,
  2739.                                                        EOPHistory eopHistory,
  2740.                                                        TimeScales timeScales);

  2742.     /** Get the function computing tidal corrections for Earth Orientation Parameters.
  2743.      *
  2744.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2745.      *
  2746.      * @return function computing tidal corrections for Earth Orientation Parameters,
  2747.      * for xp, yp, ut1 and lod respectively
  2748.      * @since 6.1
  2749.      * @see #getEOPTidalCorrection(TimeScales)
  2750.      */
  2751.     @DefaultDataContext
  2752.     public TimeVectorFunction getEOPTidalCorrection() {
  2753.         return getEOPTidalCorrection(DataContext.getDefault().getTimeScales());
  2754.     }

  2756.     /**
  2757.      * Get the function computing tidal corrections for Earth Orientation Parameters.
  2758.      *
  2759.      * @param timeScales used in the computation. The TT and TAI scales are used.
  2760.      * @return function computing tidal corrections for Earth Orientation Parameters, for
  2761.      * xp, yp, ut1 and lod respectively
  2762.      * @since 10.1
  2763.      */
  2764.     public abstract TimeVectorFunction getEOPTidalCorrection(TimeScales timeScales);

  2766.     /** Get the Love numbers.
  2767.      * @return Love numbers
  2768.           * @since 6.1
  2769.      */
  2770.     public abstract LoveNumbers getLoveNumbers();

  2772.     /** Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2773.      *
  2774.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2775.      *
  2776.      * @param ut1 UT1 time scale
  2777.      * @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2778.      * @since 6.1
  2779.      * @see #getTideFrequencyDependenceFunction(TimeScale, TimeScales)
  2780.      */
  2781.     @DefaultDataContext
  2782.     public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1) {
  2783.         return getTideFrequencyDependenceFunction(ut1,
  2784.                 DataContext.getDefault().getTimeScales());
  2785.     }

  2787.     /**
  2788.      * Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂,
  2789.      * ΔS₂₂).
  2790.      *
  2791.      * @param ut1 UT1 time scale
  2792.      * @param timeScales other time scales used in the computation including TAI and TT.
  2793.      * @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂,
  2794.      * ΔS₂₂).
  2795.      * @since 10.1
  2796.      */
  2797.     public abstract TimeVectorFunction getTideFrequencyDependenceFunction(
  2798.             TimeScale ut1,
  2799.             TimeScales timeScales);

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

  2806.     /** Get the function computing solid pole tide (ΔC₂₁, ΔS₂₁).
  2807.      * @param eopHistory EOP history
  2808.      * @return model for solid pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2809.           * @since 6.1
  2810.      */
  2811.     public abstract TimeVectorFunction getSolidPoleTide(EOPHistory eopHistory);

  2813.     /** Get the function computing ocean pole tide (ΔC₂₁, ΔS₂₁).
  2814.      * @param eopHistory EOP history
  2815.      * @return model for ocean pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2816.           * @since 6.1
  2817.      */
  2818.     public abstract TimeVectorFunction getOceanPoleTide(EOPHistory eopHistory);

  2820.     /** Get the nominal values of the displacement numbers.
  2821.      * @return an array containing h⁽⁰⁾, h⁽²⁾, h₃, hI diurnal, hI semi-diurnal,
  2822.      * l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾, l₃, lI diurnal, lI semi-diurnal,
  2823.      * H₀ permanent deformation amplitude
  2824.      * @since 9.1
  2825.      */
  2826.     public abstract double[] getNominalTidalDisplacement();

  2828.     /** Get the correction function for tidal displacement for diurnal tides.
  2829.      * <ul>
  2830.      *  <li>f[0]: radial correction, longitude cosine part</li>
  2831.      *  <li>f[1]: radial correction, longitude sine part</li>
  2832.      *  <li>f[2]: North correction, longitude cosine part</li>
  2833.      *  <li>f[3]: North correction, longitude sine part</li>
  2834.      *  <li>f[4]: East correction, longitude cosine part</li>
  2835.      *  <li>f[5]: East correction, longitude sine part</li>
  2836.      * </ul>
  2837.      * @return correction function for tidal displacement
  2838.      * @since 9.1
  2839.      */
  2840.     public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal();

  2842.     /** Get the correction function for tidal displacement for diurnal tides.
  2843.      * <ul>
  2844.      *  <li>f[0]: radial correction, longitude cosine part</li>
  2845.      *  <li>f[1]: radial correction, longitude sine part</li>
  2846.      *  <li>f[2]: North correction, longitude cosine part</li>
  2847.      *  <li>f[3]: North correction, longitude sine part</li>
  2848.      *  <li>f[4]: East correction, longitude cosine part</li>
  2849.      *  <li>f[5]: East correction, longitude sine part</li>
  2850.      * </ul>
  2851.      * @param tableName name for the diurnal tides table
  2852.      * @param cols total number of columns of the diurnal tides table
  2853.      * @param rIp column holding ∆Rf(ip) in the diurnal tides table, counting from 1
  2854.      * @param rOp column holding ∆Rf(op) in the diurnal tides table, counting from 1
  2855.      * @param tIp column holding ∆Tf(ip) in the diurnal tides table, counting from 1
  2856.      * @param tOp column holding ∆Tf(op) in the diurnal tides table, counting from 1
  2857.      * @return correction function for tidal displacement for diurnal tides
  2858.           * @since 9.1
  2859.      */
  2860.     protected static CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal(final String tableName, final int cols,
  2861.                                                                                    final int rIp, final int rOp,
  2862.                                                                                    final int tIp, final int tOp) {

  2864.         try (InputStream in1 = getStream(tableName);
  2865.              InputStream in2 = getStream(tableName);
  2866.              InputStream in3 = getStream(tableName);
  2867.              InputStream in4 = getStream(tableName);
  2868.              InputStream in5 = getStream(tableName);
  2869.              InputStream in6 = getStream(tableName)) {
  2870.             // radial component, missing the sin 2φ factor; this corresponds to:
  2871.             //  - equation 15a in IERS conventions 1996, chapter 7
  2872.             //  - equation 16a in IERS conventions 2003, chapter 7
  2873.             //  - equation 7.12a in IERS conventions 2010, chapter 7
  2874.             final PoissonSeries drCos = new PoissonSeriesParser(cols).
  2875.                                         withOptionalColumn(1).
  2876.                                         withDoodson(4, 3).
  2877.                                         withFirstDelaunay(10).
  2878.                                         withSinCos(0, rIp, +1.0e-3, rOp, +1.0e-3).
  2879.                                         parse(in1, tableName);
  2880.             final PoissonSeries drSin = new PoissonSeriesParser(cols).
  2881.                                         withOptionalColumn(1).
  2882.                                         withDoodson(4, 3).
  2883.                                         withFirstDelaunay(10).
  2884.                                         withSinCos(0, rOp, -1.0e-3, rIp, +1.0e-3).
  2885.                                         parse(in2, tableName);

  2887.             // North component, missing the cos 2φ factor; this corresponds to:
  2888.             //  - equation 15b in IERS conventions 1996, chapter 7
  2889.             //  - equation 16b in IERS conventions 2003, chapter 7
  2890.             //  - equation 7.12b in IERS conventions 2010, chapter 7
  2891.             final PoissonSeries dnCos = new PoissonSeriesParser(cols).
  2892.                                         withOptionalColumn(1).
  2893.                                         withDoodson(4, 3).
  2894.                                         withFirstDelaunay(10).
  2895.                                         withSinCos(0, tIp, +1.0e-3, tOp, +1.0e-3).
  2896.                                         parse(in3, tableName);
  2897.             final PoissonSeries dnSin = new PoissonSeriesParser(cols).
  2898.                                         withOptionalColumn(1).
  2899.                                         withDoodson(4, 3).
  2900.                                         withFirstDelaunay(10).
  2901.                                         withSinCos(0, tOp, -1.0e-3, tIp, +1.0e-3).
  2902.                                         parse(in4, tableName);

  2904.             // East component, missing the sin φ factor; this corresponds to:
  2905.             //  - equation 15b in IERS conventions 1996, chapter 7
  2906.             //  - equation 16b in IERS conventions 2003, chapter 7
  2907.             //  - equation 7.12b in IERS conventions 2010, chapter 7
  2908.             final PoissonSeries deCos = new PoissonSeriesParser(cols).
  2909.                                         withOptionalColumn(1).
  2910.                                         withDoodson(4, 3).
  2911.                                         withFirstDelaunay(10).
  2912.                                         withSinCos(0, tOp, -1.0e-3, tIp, +1.0e-3).
  2913.                                         parse(in5, tableName);
  2914.             final PoissonSeries deSin = new PoissonSeriesParser(cols).
  2915.                                         withOptionalColumn(1).
  2916.                                         withDoodson(4, 3).
  2917.                                         withFirstDelaunay(10).
  2918.                                         withSinCos(0, tIp, -1.0e-3, tOp, -1.0e-3).
  2919.                                         parse(in6, tableName);

  2921.             return PoissonSeries.compile(drCos, drSin,
  2922.                                          dnCos, dnSin,
  2923.                                          deCos, deSin);
  2924.         } catch (IOException e) {
  2925.             throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2926.         }

  2928.     }

  2930.     /** Get the correction function for tidal displacement for zonal tides.
  2931.      * <ul>
  2932.      *  <li>f[0]: radial correction</li>
  2933.      *  <li>f[1]: North correction</li>
  2934.      * </ul>
  2935.      * @return correction function for tidal displacement
  2936.      * @since 9.1
  2937.      */
  2938.     public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionZonal();

  2940.     /** Get the correction function for tidal displacement for zonal tides.
  2941.      * <ul>
  2942.      *  <li>f[0]: radial correction</li>
  2943.      *  <li>f[1]: North correction</li>
  2944.      * </ul>
  2945.      * @param tableName name for the zonal tides table
  2946.      * @param cols total number of columns of the table
  2947.      * @param rIp column holding ∆Rf(ip) in the table, counting from 1
  2948.      * @param rOp column holding ∆Rf(op) in the table, counting from 1
  2949.      * @param tIp column holding ∆Tf(ip) in the table, counting from 1
  2950.      * @param tOp column holding ∆Tf(op) in the table, counting from 1
  2951.      * @return correction function for tidal displacement for zonal tides
  2952.           * @since 9.1
  2953.      */
  2954.     protected static CompiledSeries getTidalDisplacementFrequencyCorrectionZonal(final String tableName, final int cols,
  2955.                                                                                  final int rIp, final int rOp,
  2956.                                                                                  final int tIp, final int tOp) {

  2958.         try (InputStream in1 = getStream(tableName);
  2959.              InputStream in2 = getStream(tableName)) {
  2960.             // radial component, missing the 3⁄2 sin² φ - 1⁄2 factor; this corresponds to:
  2961.             //  - equation 16a in IERS conventions 1996, chapter 7
  2962.             //  - equation 17a in IERS conventions 2003, chapter 7
  2963.             //  - equation 7.13a in IERS conventions 2010, chapter 7
  2964.             final PoissonSeries dr = new PoissonSeriesParser(cols).
  2965.                                      withOptionalColumn(1).
  2966.                                      withDoodson(4, 3).
  2967.                                      withFirstDelaunay(10).
  2968.                                      withSinCos(0, rOp, +1.0e-3, rIp, +1.0e-3).
  2969.                                      parse(in1, tableName);

  2971.             // North component, missing the sin 2φ factor; this corresponds to:
  2972.             //  - equation 16b in IERS conventions 1996, chapter 7
  2973.             //  - equation 17b in IERS conventions 2003, chapter 7
  2974.             //  - equation 7.13b in IERS conventions 2010, chapter 7
  2975.             final PoissonSeries dn = new PoissonSeriesParser(cols).
  2976.                                      withOptionalColumn(1).
  2977.                                      withDoodson(4, 3).
  2978.                                      withFirstDelaunay(10).
  2979.                                      withSinCos(0, tOp, +1.0e-3, tIp, +1.0e-3).
  2980.                                      parse(in2, tableName);

  2982.             return PoissonSeries.compile(dr, dn);
  2983.         } catch (IOException e) {
  2984.             throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2985.         }

  2987.     }

  2989.     /** Interface for functions converting nutation corrections between
  2990.      * δΔψ/δΔε to δX/δY.
  2991.      * <ul>
  2992.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  2993.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  2994.      * </ul>
  2995.      * @since 6.1
  2996.      */
  2997.     public interface NutationCorrectionConverter {

  2999.         /** Convert nutation corrections.
  3000.          * @param date current date
  3001.          * @param ddPsi δΔψ part of the nutation correction
  3002.          * @param ddEpsilon δΔε part of the nutation correction
  3003.          * @return array containing δX and δY
  3004.          */
  3005.         double[] toNonRotating(AbsoluteDate date, double ddPsi, double ddEpsilon);

  3007.         /** Convert nutation corrections.
  3008.          * @param date current date
  3009.          * @param dX δX part of the nutation correction
  3010.          * @param dY δY part of the nutation correction
  3011.          * @return array containing δΔψ and δΔε
  3012.          */
  3013.         double[] toEquinox(AbsoluteDate date, double dX, double dY);

  3015.     }

  3017.     /** Create a function converting nutation corrections between
  3018.      * δX/δY and δΔψ/δΔε.
  3019.      * <ul>
  3020.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  3021.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  3022.      * </ul>
  3023.      *
  3024.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  3025.      *
  3026.      * @return a new converter
  3027.      * @since 6.1
  3028.      * @see #getNutationCorrectionConverter(TimeScales)
  3029.      */
  3030.     @DefaultDataContext
  3031.     public NutationCorrectionConverter getNutationCorrectionConverter() {
  3032.         return getNutationCorrectionConverter(DataContext.getDefault().getTimeScales());
  3033.     }

  3035.     /** Create a function converting nutation corrections between
  3036.      * δX/δY and δΔψ/δΔε.
  3037.      * <ul>
  3038.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  3039.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  3040.      * </ul>
  3041.      * @return a new converter
  3042.      * @since 10.1
  3043.      * @param timeScales used to define the conversion.
  3044.      */
  3045.     public NutationCorrectionConverter getNutationCorrectionConverter(
  3046.             final TimeScales timeScales) {

  3048.         // get models parameters
  3049.         final TimeVectorFunction precessionFunction = getPrecessionFunction(timeScales);
  3050.         final TimeScalarFunction epsilonAFunction = getMeanObliquityFunction(timeScales);
  3051.         final AbsoluteDate date0 = getNutationReferenceEpoch(timeScales);
  3052.         final double cosE0 = FastMath.cos(epsilonAFunction.value(date0));

  3054.         return new NutationCorrectionConverter() {

  3056.             /** {@inheritDoc} */
  3057.             @Override
  3058.             public double[] toNonRotating(final AbsoluteDate date,
  3059.                                           final double ddPsi, final double ddEpsilon) {
  3060.                 // compute precession angles psiA, omegaA and chiA
  3061.                 final double[] angles = precessionFunction.value(date);

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

  3067.                 // convert nutation corrections (equation 23/IERS-2003 or 5.25/IERS-2010)
  3068.                 return new double[] {
  3069.                     sinEA * ddPsi + c * ddEpsilon,
  3070.                     ddEpsilon - c * sinEA * ddPsi
  3071.                 };

  3073.             }

  3075.             /** {@inheritDoc} */
  3076.             @Override
  3077.             public double[] toEquinox(final AbsoluteDate date,
  3078.                                       final double dX, final double dY) {
  3079.                 // compute precession angles psiA, omegaA and chiA
  3080.                 final double[] angles   = precessionFunction.value(date);

  3082.                 // conversion coefficients
  3083.                 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
  3084.                 final double c     = angles[0] * cosE0 - angles[2];
  3085.                 final double opC2  = 1 + c * c;

  3087.                 // convert nutation corrections (inverse of equation 23/IERS-2003 or 5.25/IERS-2010)
  3088.                 return new double[] {
  3089.                     (dX - c * dY) / (sinEA * opC2),
  3090.                     (dY + c * dX) / opC2
  3091.                 };
  3092.             }

  3094.         };

  3096.     }

  3098.     /** Load the Love numbers.
  3099.      * @param nameLove name of the Love number resource
  3100.      * @return Love numbers
  3101.      */
  3102.     protected LoveNumbers loadLoveNumbers(final String nameLove) {
  3103.         try {

  3105.             // allocate the three triangular arrays for real, imaginary and time-dependent numbers
  3106.             final double[][] real      = new double[4][];
  3107.             final double[][] imaginary = new double[4][];
  3108.             final double[][] plus      = new double[4][];
  3109.             for (int i = 0; i < real.length; ++i) {
  3110.                 real[i]      = new double[i + 1];
  3111.                 imaginary[i] = new double[i + 1];
  3112.                 plus[i]      = new double[i + 1];
  3113.             }

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

  3117.                 if (stream == null) {
  3118.                     // this should never happen with files embedded within Orekit
  3119.                     throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, nameLove);
  3120.                 }

  3122.                 int lineNumber = 1;
  3123.                 String line = null;
  3124.                 // setup the reader
  3125.                 try (BufferedReader reader = new BufferedReader(new InputStreamReader(stream, StandardCharsets.UTF_8))) {

  3127.                     line = reader.readLine();

  3129.                     // look for the Love numbers
  3130.                     while (line != null) {

  3132.                         line = line.trim();
  3133.                         if (!(line.isEmpty() || line.startsWith("#"))) {
  3134.                             final String[] fields = SEPARATOR.split(line);
  3135.                             if (fields.length != 5) {
  3136.                                 // this should never happen with files embedded within Orekit
  3137.                                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  3138.                                                           lineNumber, nameLove, line);
  3139.                             }
  3140.                             final int n = Integer.parseInt(fields[0]);
  3141.                             final int m = Integer.parseInt(fields[1]);
  3142.                             if ((n < 2) || (n > 3) || (m < 0) || (m > n)) {
  3143.                                 // this should never happen with files embedded within Orekit
  3144.                                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  3145.                                                           lineNumber, nameLove, line);

  3147.                             }
  3148.                             real[n][m]      = Double.parseDouble(fields[2]);
  3149.                             imaginary[n][m] = Double.parseDouble(fields[3]);
  3150.                             plus[n][m]      = Double.parseDouble(fields[4]);
  3151.                         }

  3153.                         // next line
  3154.                         lineNumber++;
  3155.                         line = reader.readLine();

  3157.                     }
  3158.                 } catch (NumberFormatException nfe) {
  3159.                     // this should never happen with files embedded within Orekit
  3160.                     throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  3161.                                               lineNumber, nameLove, line);
  3162.                 }
  3163.             }

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

  3167.         } catch (IOException ioe) {
  3168.             // this should never happen with files embedded within Orekit
  3169.             throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, nameLove);
  3170.         }
  3171.     }

  3173.     /** Get a data stream.
  3174.      * @param name file name of the resource stream
  3175.      * @return stream
  3176.      */
  3177.     private static InputStream getStream(final String name) {
  3178.         return IERSConventions.class.getResourceAsStream(name);
  3179.     }

  3181.     /** Correction to equation of equinoxes.
  3182.      * <p>IAU 1994 resolution C7 added two terms to the equation of equinoxes
  3183.      * taking effect since 1997-02-27 for continuity
  3184.      * </p>
  3185.      */
  3186.     private static class IAU1994ResolutionC7 {

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

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


  3195.         /** Evaluate the correction.
  3196.          * @param arguments Delaunay for nutation
  3197.          * @param tai TAI time scale.
  3198.          * @return correction value (0 before 1997-02-27)
  3199.          */
  3200.         public static double value(final DelaunayArguments arguments,
  3201.                                    final TimeScale tai) {
  3202.             /* Start date for applying Moon corrections to the equation of the equinoxes.
  3203.              * This date corresponds to 1997-02-27T00:00:00 UTC, hence the 30s offset from TAI.
  3204.              */
  3205.             final AbsoluteDate modelStart = new AbsoluteDate(1997, 2, 27, 0, 0, 30, tai);
  3206.             if (arguments.getDate().compareTo(modelStart) >= 0) {

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

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

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

  3217.             } else {
  3218.                 return 0.0;
  3219.             }
  3220.         }

  3222.         /** Evaluate the correction.
  3223.          * @param arguments Delaunay for nutation
  3224.          * @param tai TAI time scale.
  3225.          * @param <T> type of the field elements
  3226.          * @return correction value (0 before 1997-02-27)
  3227.          */
  3228.         public static <T extends RealFieldElement<T>> T value(
  3229.                 final FieldDelaunayArguments<T> arguments,
  3230.                 final TimeScale tai) {
  3231.             /* Start date for applying Moon corrections to the equation of the equinoxes.
  3232.              * This date corresponds to 1997-02-27T00:00:00 UTC, hence the 30s offset from TAI.
  3233.              */
  3234.             final AbsoluteDate modelStart = new AbsoluteDate(1997, 2, 27, 0, 0, 30, tai);
  3235.             if (arguments.getDate().toAbsoluteDate().compareTo(modelStart) >= 0) {

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

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

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

  3246.             } else {
  3247.                 return arguments.getDate().getField().getZero();
  3248.             }
  3249.         }

  3251.     };

  3253.     /** Stellar angle model.
  3254.      * <p>
  3255.      * The stellar angle computed here has been defined in the paper "A non-rotating origin on the
  3256.      * instantaneous equator: Definition, properties and use", N. Capitaine, Guinot B., and Souchay J.,
  3257.      * Celestial Mechanics, Volume 39, Issue 3, pp 283-307. It has been proposed as a conventional
  3258.      * conventional relationship between UT1 and Earth rotation in the paper "Definition of the Celestial
  3259.      * Ephemeris origin and of UT1 in the International Celestial Reference Frame", Capitaine, N.,
  3260.      * Guinot, B., and McCarthy, D. D., 2000, Astronomy and Astrophysics, 355(1), pp. 398–405.
  3261.      * </p>
  3262.      * <p>
  3263.      * It is presented simply as stellar angle in IERS conventions 1996 but as since been adopted as
  3264.      * the conventional relationship defining UT1 from ICRF and is called Earth Rotation Angle in
  3265.      * IERS conventions 2003 and 2010.
  3266.      * </p>
  3267.      */
  3268.     private static class StellarAngleCapitaine implements TimeScalarFunction {

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

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

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

  3281.         /** UT1 time scale. */
  3282.         private final TimeScale ut1;

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

  3287.         /** Simple constructor.
  3288.          * @param ut1 UT1 time scale
  3289.          * @param tai TAI time scale
  3290.          */
  3291.         StellarAngleCapitaine(final TimeScale ut1, final TimeScale tai) {
  3292.             this.ut1 = ut1;
  3293.             referenceDate = new AbsoluteDate(
  3294.                     DateComponents.J2000_EPOCH,
  3295.                     TimeComponents.H12,
  3296.                     tai);
  3297.         }

  3299.         /** Get the rotation rate.
  3300.          * @return rotation rate
  3301.          */
  3302.         public double getRate() {
  3303.             return ERA_1A + ERA_1B;
  3304.         }

  3306.         /** {@inheritDoc} */
  3307.         @Override
  3308.         public double value(final AbsoluteDate date) {

  3310.             // split the date offset as a full number of days plus a smaller part
  3311.             final int secondsInDay = 86400;
  3312.             final double dt  = date.durationFrom(referenceDate);
  3313.             final long days  = ((long) dt) / secondsInDay;
  3314.             final double dtA = secondsInDay * days;
  3315.             final double dtB = (dt - dtA) + ut1.offsetFromTAI(date);

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

  3319.         }

  3321.         /** {@inheritDoc} */
  3322.         @Override
  3323.         public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  3325.             // split the date offset as a full number of days plus a smaller part
  3326.             final int secondsInDay = 86400;
  3327.             final T dt  = date.durationFrom(referenceDate);
  3328.             final long days  = ((long) dt.getReal()) / secondsInDay;
  3329.             final double dtA = secondsInDay * days;
  3330.             final T dtB = dt.subtract(dtA).add(ut1.offsetFromTAI(date.toAbsoluteDate()));

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

  3334.         }

  3336.     }

  3338.     /** Mean pole. */
  3339.     private static class MeanPole implements TimeStamped, Serializable {

  3341.         /** Serializable UID. */
  3342.         private static final long serialVersionUID = 20131028l;

  3344.         /** Date. */
  3345.         private final AbsoluteDate date;

  3347.         /** X coordinate. */
  3348.         private double x;

  3350.         /** Y coordinate. */
  3351.         private double y;

  3353.         /** Simple constructor.
  3354.          * @param date date
  3355.          * @param x x coordinate
  3356.          * @param y y coordinate
  3357.          */
  3358.         MeanPole(final AbsoluteDate date, final double x, final double y) {
  3359.             this.date = date;
  3360.             this.x    = x;
  3361.             this.y    = y;
  3362.         }

  3364.         /** {@inheritDoc} */
  3365.         @Override
  3366.         public AbsoluteDate getDate() {
  3367.             return date;
  3368.         }

  3370.         /** Get x coordinate.
  3371.          * @return x coordinate
  3372.          */
  3373.         public double getX() {
  3374.             return x;
  3375.         }

  3377.         /** Get y coordinate.
  3378.          * @return y coordinate
  3379.          */
  3380.         public double getY() {
  3381.             return y;
  3382.         }

  3384.     }

  3386.     /** Local class for precession function. */
  3387.     private class PrecessionFunction implements TimeVectorFunction {

  3389.         /** Polynomial nutation for psiA. */
  3390.         private final PolynomialNutation psiA;

  3392.         /** Polynomial nutation for omegaA. */
  3393.         private final PolynomialNutation omegaA;

  3395.         /** Polynomial nutation for chiA. */
  3396.         private final PolynomialNutation chiA;

  3398.         /** Time scales to use. */
  3399.         private final TimeScales timeScales;

  3401.         /** Simple constructor.
  3402.          * @param psiA polynomial nutation for psiA
  3403.          * @param omegaA polynomial nutation for omegaA
  3404.          * @param chiA polynomial nutation for chiA
  3405.          * @param timeScales used in the computation.
  3406.          */
  3407.         PrecessionFunction(final PolynomialNutation psiA,
  3408.                            final PolynomialNutation omegaA,
  3409.                            final PolynomialNutation chiA,
  3410.                            final TimeScales timeScales) {
  3411.             this.psiA   = psiA;
  3412.             this.omegaA = omegaA;
  3413.             this.chiA   = chiA;
  3414.             this.timeScales = timeScales;
  3415.         }


  3418.         /** {@inheritDoc} */
  3419.         @Override
  3420.         public double[] value(final AbsoluteDate date) {
  3421.             final double tc = evaluateTC(date, timeScales);
  3422.             return new double[] {
  3423.                 psiA.value(tc), omegaA.value(tc), chiA.value(tc)
  3424.             };
  3425.         }

  3427.         /** {@inheritDoc} */
  3428.         @Override
  3429.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  3430.             final T[] a = MathArrays.buildArray(date.getField(), 3);
  3431.             final T tc = evaluateTC(date, timeScales);
  3432.             a[0] = psiA.value(tc);
  3433.             a[1] = omegaA.value(tc);
  3434.             a[2] = chiA.value(tc);
  3435.             return a;
  3436.         }

  3438.     }

  3440.     /** Local class for tides frequency function. */
  3441.     private static class TideFrequencyDependenceFunction implements TimeVectorFunction {

  3443.         /** Nutation arguments. */
  3444.         private final FundamentalNutationArguments arguments;

  3446.         /** Correction series. */
  3447.         private final PoissonSeries.CompiledSeries kSeries;

  3449.         /** Simple constructor.
  3450.          * @param arguments nutation arguments
  3451.          * @param c20Series correction series for the C20 term
  3452.          * @param c21Series correction series for the C21 term
  3453.          * @param s21Series correction series for the S21 term
  3454.          * @param c22Series correction series for the C22 term
  3455.          * @param s22Series correction series for the S22 term
  3456.          */
  3457.         TideFrequencyDependenceFunction(final FundamentalNutationArguments arguments,
  3458.                                         final PoissonSeries c20Series,
  3459.                                         final PoissonSeries c21Series,
  3460.                                         final PoissonSeries s21Series,
  3461.                                         final PoissonSeries c22Series,
  3462.                                         final PoissonSeries s22Series) {
  3463.             this.arguments = arguments;
  3464.             this.kSeries   = PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);
  3465.         }


  3468.         /** {@inheritDoc} */
  3469.         @Override
  3470.         public double[] value(final AbsoluteDate date) {
  3471.             return kSeries.value(arguments.evaluateAll(date));
  3472.         }

  3474.         /** {@inheritDoc} */
  3475.         @Override
  3476.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  3477.             return kSeries.value(arguments.evaluateAll(date));
  3478.         }

  3480.     }

  3482.     /** Local class for EOP tidal corrections. */
  3483.     private static class EOPTidalCorrection implements TimeVectorFunction {

  3485.         /** Nutation arguments. */
  3486.         private final FundamentalNutationArguments arguments;

  3488.         /** Correction series. */
  3489.         private final PoissonSeries.CompiledSeries correctionSeries;

  3491.         /** Simple constructor.
  3492.          * @param arguments nutation arguments
  3493.          * @param xSeries correction series for the x coordinate
  3494.          * @param ySeries correction series for the y coordinate
  3495.          * @param ut1Series correction series for the UT1
  3496.          */
  3497.         EOPTidalCorrection(final FundamentalNutationArguments arguments,
  3498.                            final PoissonSeries xSeries,
  3499.                            final PoissonSeries ySeries,
  3500.                            final PoissonSeries ut1Series) {
  3501.             this.arguments        = arguments;
  3502.             this.correctionSeries = PoissonSeries.compile(xSeries, ySeries, ut1Series);
  3503.         }

  3505.         /** {@inheritDoc} */
  3506.         @Override
  3507.         public double[] value(final AbsoluteDate date) {
  3508.             final BodiesElements elements = arguments.evaluateAll(date);
  3509.             final double[] correction    = correctionSeries.value(elements);
  3510.             final double[] correctionDot = correctionSeries.derivative(elements);
  3511.             return new double[] {
  3512.                 correction[0],
  3513.                 correction[1],
  3514.                 correction[2],
  3515.                 correctionDot[2] * (-Constants.JULIAN_DAY)
  3516.             };
  3517.         }

  3519.         /** {@inheritDoc} */
  3520.         @Override
  3521.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

  3523.             final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  3524.             final T[] correction    = correctionSeries.value(elements);
  3525.             final T[] correctionDot = correctionSeries.derivative(elements);
  3526.             final T[] a = MathArrays.buildArray(date.getField(), 4);
  3527.             a[0] = correction[0];
  3528.             a[1] = correction[1];
  3529.             a[2] = correction[2];
  3530.             a[3] = correctionDot[2].multiply(-Constants.JULIAN_DAY);
  3531.             return a;
  3532.         }

  3534.     }

  3536. }
Created with JaCoCo 0.8.5.201910111838