IERSConventions.java

  1. /* Copyright 2002-2025 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.function.Function;
  27. import java.util.regex.Pattern;

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


  69. /** Supported IERS conventions.
  70.  * @since 6.0
  71.  * @author Luc Maisonobe
  72.  */
  73. public enum IERSConventions {

  75.     /** Constant for IERS 1996 conventions. */
  76.     IERS_1996 {

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

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

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

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

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

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

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

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

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

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

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

  111.         /** {@inheritDoc} */
  112.         @Override
  113.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
  114.                                                                  final TimeScales timeScales) {
  115.             return load(NUTATION_ARGUMENTS, in -> new FundamentalNutationArguments(this, timeScale,
  116.                                                                                             in, NUTATION_ARGUMENTS,
  117.                                                                                             timeScales));
  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 CalculusFieldElement<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 = load(X_Y_SERIES, in -> xParser.parse(in, X_Y_SERIES));

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

  193.             final double fYCosOm    = -0.00231 * Constants.ARC_SECONDS_TO_RADIANS;
  194.             final double fYCos2FDOm = -0.00014 * Constants.ARC_SECONDS_TO_RADIANS;

  196.             final PoissonSeriesParser yParser =
  197.                     baseParser.
  198.                     withSinCos(0, -1, deciMilliAS, 10, deciMilliAS).
  199.                     withSinCos(1, 12, deciMilliAS, 11, deciMilliAS);
  200.             final PoissonSeries ySum = load(X_Y_SERIES, in -> yParser.parse(in, X_Y_SERIES));

  202.             final PoissonSeries.CompiledSeries xySum =
  203.                     PoissonSeries.compile(xSum, ySum);

  205.             // s = -XY/2 + 0.00385″t - 0.07259″t³ - 0.00264″ sin Ω - 0.00006″ sin 2Ω
  206.             //     + 0.00074″t² sin Ω + 0.00006″t² sin 2(F - D + Ω)
  207.             final double fST          =  0.00385 * Constants.ARC_SECONDS_TO_RADIANS;
  208.             final double fST3         = -0.07259 * Constants.ARC_SECONDS_TO_RADIANS;
  209.             final double fSSinOm      = -0.00264 * Constants.ARC_SECONDS_TO_RADIANS;
  210.             final double fSSin2Om     = -0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
  211.             final double fST2SinOm    =  0.00074 * Constants.ARC_SECONDS_TO_RADIANS;
  212.             final double fST2Sin2FDOm =  0.00006 * Constants.ARC_SECONDS_TO_RADIANS;

  214.             return new TimeVectorFunction() {

  216.                 /** {@inheritDoc} */
  217.                 @Override
  218.                 public double[] value(final AbsoluteDate date) {

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

  223.                     final double omega     = elements.getOmega();
  224.                     final double f         = elements.getF();
  225.                     final double d         = elements.getD();
  226.                     final double t         = elements.getTC();

  228.                     final SinCos scOmega   = FastMath.sinCos(omega);
  229.                     final SinCos sc2omega  = SinCos.sum(scOmega, scOmega);
  230.                     final SinCos sc2FD0m   = FastMath.sinCos(2 * (f - d + omega));
  231.                     final double cosOmega  = scOmega.cos();
  232.                     final double sinOmega  = scOmega.sin();
  233.                     final double sin2Omega = sc2omega.sin();
  234.                     final double cos2FDOm  = sc2FD0m.cos();
  235.                     final double sin2FDOm  = sc2FD0m.sin();

  237.                     final double x = xPolynomial.value(t) + sinEps0 * xy[0] +
  238.                             t * t * (fXCosOm * cosOmega + fXSinOm * sinOmega + fXSin2FDOm * cos2FDOm);
  239.                     final double y = yPolynomial.value(t) + xy[1] +
  240.                             t * t * (fYCosOm * cosOmega + fYCos2FDOm * cos2FDOm);
  241.                     final double sPxy2 = fSSinOm * sinOmega + fSSin2Om * sin2Omega +
  242.                             t * (fST + t * (fST2SinOm * sinOmega + fST2Sin2FDOm * sin2FDOm + t * fST3));

  244.                     return new double[] {
  245.                         x, y, sPxy2
  246.                     };

  248.                 }

  250.                 /** {@inheritDoc} */
  251.                 @Override
  252.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

  257.                     final T omega     = elements.getOmega();
  258.                     final T f         = elements.getF();
  259.                     final T d         = elements.getD();
  260.                     final T t         = elements.getTC();
  261.                     final T t2        = t.square();

  263.                     final FieldSinCos<T> scOmega  = FastMath.sinCos(omega);
  264.                     final FieldSinCos<T> sc2omega = FieldSinCos.sum(scOmega, scOmega);
  265.                     final FieldSinCos<T> sc2FD0m  = FastMath.sinCos(f.subtract(d).add(omega).multiply(2));
  266.                     final T cosOmega  = scOmega.cos();
  267.                     final T sinOmega  = scOmega.sin();
  268.                     final T sin2Omega = sc2omega.sin();
  269.                     final T cos2FDOm  = sc2FD0m.cos();
  270.                     final T sin2FDOm  = sc2FD0m.sin();

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

  282.                     final T[] a = MathArrays.buildArray(date.getField(), 3);
  283.                     a[0] = x;
  284.                     a[1] = y;
  285.                     a[2] = sPxy2;
  286.                     return a;

  288.                 }

  290.             };

  292.         }

  294.         /** {@inheritDoc} */
  295.         @Override
  296.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

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

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

  322.         }

  324.         /** {@inheritDoc} */
  325.         @Override
  326.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  328.             // set up nutation arguments
  329.             final FundamentalNutationArguments arguments =
  330.                     getNutationArguments(timeScales);

  332.             // set up Poisson series
  333.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  334.             final PoissonSeriesParser baseParser =
  335.                     new PoissonSeriesParser(10).withFirstDelaunay(1);

  337.             final PoissonSeriesParser psiParser =
  338.                     baseParser.
  339.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  340.                     withSinCos(1, 8, deciMilliAS, -1, deciMilliAS);
  341.             final PoissonSeries psiSeries = load(PSI_EPSILON_SERIES, in -> psiParser.parse(in, PSI_EPSILON_SERIES));

  343.             final PoissonSeriesParser epsilonParser =
  344.                     baseParser.
  345.                     withSinCos(0, -1, deciMilliAS, 9, deciMilliAS).
  346.                     withSinCos(1, -1, deciMilliAS, 10, deciMilliAS);
  347.             final PoissonSeries epsilonSeries = load(PSI_EPSILON_SERIES, in -> epsilonParser.parse(in, PSI_EPSILON_SERIES));

  349.             final PoissonSeries.CompiledSeries psiEpsilonSeries =
  350.                     PoissonSeries.compile(psiSeries, epsilonSeries);

  352.             return new TimeVectorFunction() {

  354.                 /** {@inheritDoc} */
  355.                 @Override
  356.                 public double[] value(final AbsoluteDate date) {
  357.                     final BodiesElements elements = arguments.evaluateAll(date);
  358.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  359.                     return new double[] {
  360.                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  361.                     };
  362.                 }

  364.                 /** {@inheritDoc} */
  365.                 @Override
  366.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  367.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  368.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  369.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  370.                     result[0] = psiEpsilon[0];
  371.                     result[1] = psiEpsilon[1];
  372.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  373.                     return result;
  374.                 }

  376.             };

  378.         }

  380.         /** {@inheritDoc} */
  381.         @Override
  382.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  383.                                                   final TimeScales timeScales) {

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

  388.             // constants from IERS 1996 page 21
  389.             // the underlying model is IAU 1982 GMST-UT1
  390.             final AbsoluteDate gmstReference = new AbsoluteDate(
  391.                     DateComponents.J2000_EPOCH, TimeComponents.H12, timeScales.getTAI());
  392.             final double gmst0 = 24110.54841;
  393.             final double gmst1 = 8640184.812866;
  394.             final double gmst2 = 0.093104;
  395.             final double gmst3 = -6.2e-6;

  397.             return new TimeScalarFunction() {

  399.                 /** {@inheritDoc} */
  400.                 @Override
  401.                 public double value(final AbsoluteDate date) {

  403.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  404.                     final double dtai = date.durationFrom(gmstReference);
  405.                     final double tut1 = dtai + ut1.offsetFromTAI(date).toDouble();
  406.                     final double tt   = tut1 / Constants.JULIAN_CENTURY;

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

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

  415.                 }

  417.                 /** {@inheritDoc} */
  418.                 @Override
  419.                 public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  421.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  422.                     final T dtai = date.durationFrom(gmstReference);
  423.                     final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()).toDouble());
  424.                     final T tt   = tut1.divide(Constants.JULIAN_CENTURY);

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

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

  433.                 }

  435.             };

  437.         }

  439.         /** {@inheritDoc} */
  440.         @Override
  441.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  442.                                                       final TimeScales timeScales) {

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

  447.             // constants from IERS 1996 page 21
  448.             // the underlying model is IAU 1982 GMST-UT1
  449.             final AbsoluteDate gmstReference = new AbsoluteDate(
  450.                     DateComponents.J2000_EPOCH, TimeComponents.H12, timeScales.getTAI());
  451.             final double gmst1 = 8640184.812866;
  452.             final double gmst2 = 0.093104;
  453.             final double gmst3 = -6.2e-6;

  455.             return new TimeScalarFunction() {

  457.                 /** {@inheritDoc} */
  458.                 @Override
  459.                 public double value(final AbsoluteDate date) {

  461.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  462.                     final double dtai = date.durationFrom(gmstReference);
  463.                     final double tut1 = dtai + ut1.offsetFromTAI(date).toDouble();
  464.                     final double tt   = tut1 / Constants.JULIAN_CENTURY;

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

  469.                 }

  471.                 /** {@inheritDoc} */
  472.                 @Override
  473.                 public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  475.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  476.                     final T dtai = date.durationFrom(gmstReference);
  477.                     final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()).toDouble());
  478.                     final T tt   = tut1.divide(Constants.JULIAN_CENTURY);

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

  483.                 }

  485.             };

  487.         }

  489.         /** {@inheritDoc} */
  490.         @Override
  491.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  492.                                                   final EOPHistory eopHistory,
  493.                                                   final TimeScales timeScales) {

  495.             // obliquity
  496.             final TimeScalarFunction epsilonA = getMeanObliquityFunction(timeScales);

  498.             // GMST function
  499.             final TimeScalarFunction gmst = getGMSTFunction(ut1, timeScales);

  501.             // nutation function
  502.             final TimeVectorFunction nutation = getNutationFunction(timeScales);

  504.             return new TimeScalarFunction() {

  506.                 /** {@inheritDoc} */
  507.                 @Override
  508.                 public double value(final AbsoluteDate date) {

  510.                     // compute equation of equinoxes
  511.                     final double[] angles = nutation.value(date);
  512.                     double deltaPsi = angles[0];
  513.                     if (eopHistory != null) {
  514.                         deltaPsi += eopHistory.getEquinoxNutationCorrection(date)[0];
  515.                     }
  516.                     final double eqe = deltaPsi  * FastMath.cos(epsilonA.value(date)) + angles[2];

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

  521.                 }

  523.                 /** {@inheritDoc} */
  524.                 @Override
  525.                 public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  527.                     // compute equation of equinoxes
  528.                     final T[] angles = nutation.value(date);
  529.                     T deltaPsi = angles[0];
  530.                     if (eopHistory != null) {
  531.                         deltaPsi = deltaPsi.add(eopHistory.getEquinoxNutationCorrection(date)[0]);
  532.                     }
  533.                     final T eqe = deltaPsi.multiply(epsilonA.value(date).cos()).add(angles[2]);

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

  538.                 }

  540.             };

  542.         }

  544.         /** {@inheritDoc} */
  545.         @Override
  546.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

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

  556.             // set up Poisson series
  557.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  558.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(17).
  559.                     withOptionalColumn(1).
  560.                     withGamma(7).
  561.                     withFirstDelaunay(2);
  562.             final PoissonSeries xSeries =
  563.                             load(TIDAL_CORRECTION_XP_YP_SERIES, in -> xyParser.
  564.                                                                       withSinCos(0, 14, milliAS, 15, milliAS).
  565.                                                                       parse(in, TIDAL_CORRECTION_XP_YP_SERIES));
  566.             final PoissonSeries ySeries =
  567.                             load(TIDAL_CORRECTION_XP_YP_SERIES, in -> xyParser.
  568.                                                                       withSinCos(0, 16, milliAS, 17, milliAS).
  569.                                                                       parse(in, TIDAL_CORRECTION_XP_YP_SERIES));

  571.             final double deciMilliS = 1.0e-4;
  572.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(17).
  573.                     withOptionalColumn(1).
  574.                     withGamma(7).
  575.                     withFirstDelaunay(2);
  576.             final PoissonSeries ut1Series =
  577.                             load(TIDAL_CORRECTION_UT1_SERIES, in -> ut1Parser.
  578.                                                                     withSinCos(0, 16, deciMilliS, 17, deciMilliS).
  579.                                                                     parse(in, TIDAL_CORRECTION_UT1_SERIES));

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

  583.         }

  585.         /** {@inheritDoc} */
  586.         public LoveNumbers getLoveNumbers() {
  587.             return loadLoveNumbers(LOVE_NUMBERS);
  588.         }

  590.         /** {@inheritDoc} */
  591.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  592.                                                                      final TimeScales timeScales) {

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

  597.             // set up Poisson series
  598.             final PoissonSeriesParser k20Parser =
  599.                     new PoissonSeriesParser(18).
  600.                         withOptionalColumn(1).
  601.                         withDoodson(4, 2).
  602.                         withFirstDelaunay(10);
  603.             final PoissonSeriesParser k21Parser =
  604.                     new PoissonSeriesParser(18).
  605.                         withOptionalColumn(1).
  606.                         withDoodson(4, 3).
  607.                         withFirstDelaunay(10);
  608.             final PoissonSeriesParser k22Parser =
  609.                     new PoissonSeriesParser(16).
  610.                         withOptionalColumn(1).
  611.                         withDoodson(4, 2).
  612.                         withFirstDelaunay(10);

  614.             final double pico = 1.0e-12;
  615.             final PoissonSeries c20Series =
  616.                             load(K20_FREQUENCY_DEPENDENCE, in -> k20Parser.
  617.                                                                  withSinCos(0, 18, -pico, 16, pico).
  618.                                                                  parse(in, K20_FREQUENCY_DEPENDENCE));
  619.             final PoissonSeries c21Series =
  620.                             load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
  621.                                                                  withSinCos(0, 17, pico, 18, pico).
  622.                                                                  parse(in, K21_FREQUENCY_DEPENDENCE));
  623.             final PoissonSeries s21Series =
  624.                             load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
  625.                                                                  withSinCos(0, 18, -pico, 17, pico).
  626.                                                                  parse(in, K21_FREQUENCY_DEPENDENCE));
  627.             final PoissonSeries c22Series =
  628.                             load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
  629.                                                                  withSinCos(0, -1, pico, 16, pico).
  630.                                                                  parse(in, K22_FREQUENCY_DEPENDENCE));
  631.             final PoissonSeries s22Series =
  632.                             load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
  633.                                                                  withSinCos(0, 16, -pico, -1, pico).
  634.                                                                  parse(in, K22_FREQUENCY_DEPENDENCE));

  636.             return new TideFrequencyDependenceFunction(arguments,
  637.                                                        c20Series,
  638.                                                        c21Series, s21Series,
  639.                                                        c22Series, s22Series);

  641.         }

  643.         /** {@inheritDoc} */
  644.         @Override
  645.         public double getPermanentTide() {
  646.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  647.         }

  649.         /** {@inheritDoc} */
  650.         @Override
  651.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  653.             // constants from IERS 1996 page 47
  654.             final double globalFactor = -1.348e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  655.             final double coupling     =  0.0112;

  657.             return new TimeVectorFunction() {

  659.                 /** {@inheritDoc} */
  660.                 @Override
  661.                 public double[] value(final AbsoluteDate date) {
  662.                     final PoleCorrection pole = eopHistory.getPoleCorrection(date);
  663.                     return new double[] {
  664.                         globalFactor * (pole.getXp() + coupling * pole.getYp()),
  665.                         globalFactor * (coupling * pole.getXp() - pole.getYp()),
  666.                     };
  667.                 }

  669.                 /** {@inheritDoc} */
  670.                 @Override
  671.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  672.                     final FieldPoleCorrection<T> pole = eopHistory.getPoleCorrection(date);
  673.                     final T[] a = MathArrays.buildArray(date.getField(), 2);
  674.                     a[0] = pole.getXp().add(pole.getYp().multiply(coupling)).multiply(globalFactor);
  675.                     a[1] = pole.getXp().multiply(coupling).subtract(pole.getYp()).multiply(globalFactor);
  676.                     return a;
  677.                 }

  679.             };
  680.         }

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

  686.             return new TimeVectorFunction() {

  688.                 /** {@inheritDoc} */
  689.                 @Override
  690.                 public double[] value(final AbsoluteDate date) {
  691.                     // there are no model for ocean pole tide prior to conventions 2010
  692.                     return new double[] {
  693.                         0.0, 0.0
  694.                     };
  695.                 }

  697.                 /** {@inheritDoc} */
  698.                 @Override
  699.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  700.                     // there are no model for ocean pole tide prior to conventions 2010
  701.                     return MathArrays.buildArray(date.getField(), 2);
  702.                 }

  704.             };
  705.         }

  707.         /** {@inheritDoc} */
  708.         @Override
  709.         public double[] getNominalTidalDisplacement() {

  711.             //  // elastic Earth values
  712.             //  return new double[] {
  713.             //      // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  714.             //      0.6026,                                  -0.0006, 0.292, -0.0025, -0.0022,
  715.             //      // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  716.             //      0.0831,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  717.             //      // H₀
  718.             //      -0.31460
  719.             //  };

  721.             // anelastic Earth values
  722.             return new double[] {
  723.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  724.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  725.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  726.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  727.                 // H₀
  728.                 -0.31460
  729.             };

  731.         }

  733.         /** {@inheritDoc} */
  734.         @Override
  735.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
  736.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  737.                                                                   18, 17, -1, 18, -1);
  738.         }

  740.         /** {@inheritDoc} */
  741.         @Override
  742.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  743.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  744.                                                                 20, 17, 19, 18, 20);
  745.         }

  747.     },

  749.     /** Constant for IERS 2003 conventions. */
  750.     IERS_2003 {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  800.         /** {@inheritDoc} */
  801.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
  802.                                                                  final TimeScales timeScales) {
  803.             return load(NUTATION_ARGUMENTS, in -> new FundamentalNutationArguments(this, timeScale,
  804.                                                                                    in, NUTATION_ARGUMENTS,
  805.                                                                                    timeScales));
  806.         }

  808.         /** {@inheritDoc} */
  809.         @Override
  810.         public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {

  812.             // epsilon 0 value from chapter 5, page 41, other terms from equation 32 page 45
  813.             final PolynomialNutation epsilonA =
  814.                     new PolynomialNutation(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  815.                                              -46.84024  * Constants.ARC_SECONDS_TO_RADIANS,
  816.                                               -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  817.                                                0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  819.             return new TimeScalarFunction() {

  821.                 /** {@inheritDoc} */
  822.                 @Override
  823.                 public double value(final AbsoluteDate date) {
  824.                     return epsilonA.value(evaluateTC(date, timeScales));
  825.                 }

  827.                 /** {@inheritDoc} */
  828.                 @Override
  829.                 public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  830.                     return epsilonA.value(evaluateTC(date, timeScales));
  831.                 }

  833.             };

  835.         }

  837.         /** {@inheritDoc} */
  838.         @Override
  839.         public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {

  841.             // set up nutation arguments
  842.             final FundamentalNutationArguments arguments =
  843.                     getNutationArguments(timeScales);

  845.             // set up Poisson series
  846.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  847.             final PoissonSeriesParser parser =
  848.                     new PoissonSeriesParser(17).
  849.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  850.                         withFirstDelaunay(4).
  851.                         withFirstPlanetary(9).
  852.                         withSinCos(0, 2, microAS, 3, microAS);

  854.             final PoissonSeries xSeries = load(X_SERIES, in -> parser.parse(in, X_SERIES));
  855.             final PoissonSeries ySeries = load(Y_SERIES, in -> parser.parse(in, Y_SERIES));
  856.             final PoissonSeries sSeries = load(S_SERIES, in -> parser.parse(in, S_SERIES));
  857.             final PoissonSeries.CompiledSeries xys = PoissonSeries.compile(xSeries, ySeries, sSeries);

  859.             // create a function evaluating the series
  860.             return new TimeVectorFunction() {

  862.                 /** {@inheritDoc} */
  863.                 @Override
  864.                 public double[] value(final AbsoluteDate date) {
  865.                     return xys.value(arguments.evaluateAll(date));
  866.                 }

  868.                 /** {@inheritDoc} */
  869.                 @Override
  870.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  871.                     return xys.value(arguments.evaluateAll(date));
  872.                 }

  874.             };

  876.         }


  879.         /** {@inheritDoc} */
  880.         @Override
  881.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

  883.             // set up the conventional polynomials
  884.             // the following values are from equation 32 in IERS 2003 conventions
  885.             final PolynomialNutation psiA =
  886.                     new PolynomialNutation(    0.0,
  887.                                             5038.47875   * Constants.ARC_SECONDS_TO_RADIANS,
  888.                                               -1.07259   * Constants.ARC_SECONDS_TO_RADIANS,
  889.                                               -0.001147  * Constants.ARC_SECONDS_TO_RADIANS);
  890.             final PolynomialNutation omegaA =
  891.                     new PolynomialNutation(getMeanObliquityFunction(timeScales)
  892.                             .value(getNutationReferenceEpoch(timeScales)),
  893.                                            -0.02524   * Constants.ARC_SECONDS_TO_RADIANS,
  894.                                             0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
  895.                                            -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
  896.             final PolynomialNutation chiA =
  897.                     new PolynomialNutation( 0.0,
  898.                                            10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
  899.                                            -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
  900.                                            -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);

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

  904.         }

  906.         /** {@inheritDoc} */
  907.         @Override
  908.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  910.             // set up nutation arguments
  911.             final FundamentalNutationArguments arguments =
  912.                     getNutationArguments(timeScales);

  914.             // set up Poisson series
  915.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  916.             final PoissonSeriesParser luniSolarParser =
  917.                     new PoissonSeriesParser(14).withFirstDelaunay(1);
  918.             final PoissonSeriesParser luniSolarPsiParser =
  919.                     luniSolarParser.
  920.                     withSinCos(0, 7, milliAS, 11, milliAS).
  921.                     withSinCos(1, 8, milliAS, 12, milliAS);
  922.             final PoissonSeries psiLuniSolarSeries =
  923.                             load(LUNI_SOLAR_SERIES, in -> luniSolarPsiParser.parse(in, LUNI_SOLAR_SERIES));
  924.             final PoissonSeriesParser luniSolarEpsilonParser = luniSolarParser.
  925.                                                                withSinCos(0, 13, milliAS, 9, milliAS).
  926.                                                                withSinCos(1, 14, milliAS, 10, milliAS);
  927.             final PoissonSeries epsilonLuniSolarSeries =
  928.                             load(LUNI_SOLAR_SERIES, in -> luniSolarEpsilonParser.parse(in, LUNI_SOLAR_SERIES));

  930.             final PoissonSeriesParser planetaryParser =
  931.                     new PoissonSeriesParser(21).
  932.                         withFirstDelaunay(2).
  933.                         withFirstPlanetary(7);
  934.             final PoissonSeriesParser planetaryPsiParser =
  935.                     planetaryParser.withSinCos(0, 17, milliAS, 18, milliAS);
  936.             final PoissonSeries psiPlanetarySeries =
  937.                             load(PLANETARY_SERIES, in -> planetaryPsiParser.parse(in, PLANETARY_SERIES));
  938.             final PoissonSeriesParser planetaryEpsilonParser =
  939.                     planetaryParser.withSinCos(0, 19, milliAS, 20, milliAS);
  940.             final PoissonSeries epsilonPlanetarySeries =
  941.                             load(PLANETARY_SERIES, in -> planetaryEpsilonParser.parse(in, PLANETARY_SERIES));

  943.             final PoissonSeries.CompiledSeries luniSolarSeries =
  944.                     PoissonSeries.compile(psiLuniSolarSeries, epsilonLuniSolarSeries);
  945.             final PoissonSeries.CompiledSeries planetarySeries =
  946.                     PoissonSeries.compile(psiPlanetarySeries, epsilonPlanetarySeries);

  948.             return new TimeVectorFunction() {

  950.                 /** {@inheritDoc} */
  951.                 @Override
  952.                 public double[] value(final AbsoluteDate date) {
  953.                     final BodiesElements elements = arguments.evaluateAll(date);
  954.                     final double[] luniSolar = luniSolarSeries.value(elements);
  955.                     final double[] planetary = planetarySeries.value(elements);
  956.                     return new double[] {
  957.                         luniSolar[0] + planetary[0], luniSolar[1] + planetary[1],
  958.                         IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  959.                     };
  960.                 }

  962.                 /** {@inheritDoc} */
  963.                 @Override
  964.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  965.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  966.                     final T[] luniSolar = luniSolarSeries.value(elements);
  967.                     final T[] planetary = planetarySeries.value(elements);
  968.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  969.                     result[0] = luniSolar[0].add(planetary[0]);
  970.                     result[1] = luniSolar[1].add(planetary[1]);
  971.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  972.                     return result;
  973.                 }

  975.             };

  977.         }

  979.         /** {@inheritDoc} */
  980.         @Override
  981.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  982.                                                   final TimeScales timeScales) {

  984.             // Earth Rotation Angle
  985.             final StellarAngleCapitaine era =
  986.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  988.             // Polynomial part of the apparent sidereal time series
  989.             // which is the opposite of Equation of Origins (EO)
  990.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  991.             final PoissonSeriesParser parser =
  992.                     new PoissonSeriesParser(17).
  993.                         withFirstDelaunay(4).
  994.                         withFirstPlanetary(9).
  995.                         withSinCos(0, 2, microAS, 3, microAS).
  996.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  997.             final PolynomialNutation minusEO = load(GST_SERIES, in -> parser.parse(in, GST_SERIES).getPolynomial());

  999.             // create a function evaluating the series
  1000.             return new TimeScalarFunction() {

  1002.                 /** {@inheritDoc} */
  1003.                 @Override
  1004.                 public double value(final AbsoluteDate date) {
  1005.                     return era.value(date) + minusEO.value(evaluateTC(date, timeScales));
  1006.                 }

  1008.                 /** {@inheritDoc} */
  1009.                 @Override
  1010.                 public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1011.                     return era.value(date).add(minusEO.value(evaluateTC(date, timeScales)));
  1012.                 }

  1014.             };

  1016.         }

  1018.         /** {@inheritDoc} */
  1019.         @Override
  1020.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  1021.                                                       final TimeScales timeScales) {

  1023.             // Earth Rotation Angle
  1024.             final StellarAngleCapitaine era =
  1025.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  1027.             // Polynomial part of the apparent sidereal time series
  1028.             // which is the opposite of Equation of Origins (EO)
  1029.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1030.             final PoissonSeriesParser parser =
  1031.                     new PoissonSeriesParser(17).
  1032.                         withFirstDelaunay(4).
  1033.                         withFirstPlanetary(9).
  1034.                         withSinCos(0, 2, microAS, 3, microAS).
  1035.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1036.             final PolynomialNutation minusEO = load(GST_SERIES, in -> parser.parse(in, GST_SERIES).getPolynomial());

  1038.             // create a function evaluating the series
  1039.             return new TimeScalarFunction() {

  1041.                 /** {@inheritDoc} */
  1042.                 @Override
  1043.                 public double value(final AbsoluteDate date) {
  1044.                     return era.getRate() + minusEO.derivative(evaluateTC(date, timeScales));
  1045.                 }

  1047.                 /** {@inheritDoc} */
  1048.                 @Override
  1049.                 public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1050.                     return minusEO.derivative(evaluateTC(date, timeScales)).add(era.getRate());
  1051.                 }

  1053.             };

  1055.         }

  1057.         /** {@inheritDoc} */
  1058.         @Override
  1059.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  1060.                                                   final EOPHistory eopHistory,
  1061.                                                   final TimeScales timeScales) {

  1063.             // set up nutation arguments
  1064.             final FundamentalNutationArguments arguments =
  1065.                     getNutationArguments(timeScales);

  1067.             // mean obliquity function
  1068.             final TimeScalarFunction epsilon = getMeanObliquityFunction(timeScales);

  1070.             // set up Poisson series
  1071.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  1072.             final PoissonSeriesParser luniSolarPsiParser =
  1073.                     new PoissonSeriesParser(14).
  1074.                         withFirstDelaunay(1).
  1075.                         withSinCos(0, 7, milliAS, 11, milliAS).
  1076.                         withSinCos(1, 8, milliAS, 12, milliAS);
  1077.             final PoissonSeries psiLuniSolarSeries =
  1078.                             load(LUNI_SOLAR_SERIES, in -> luniSolarPsiParser.parse(in, LUNI_SOLAR_SERIES));

  1080.             final PoissonSeriesParser planetaryPsiParser =
  1081.                     new PoissonSeriesParser(21).
  1082.                         withFirstDelaunay(2).
  1083.                         withFirstPlanetary(7).
  1084.                         withSinCos(0, 17, milliAS, 18, milliAS);
  1085.             final PoissonSeries psiPlanetarySeries =
  1086.                             load(PLANETARY_SERIES, in -> planetaryPsiParser.parse(in, PLANETARY_SERIES));


  1089.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1090.             final PoissonSeriesParser gstParser =
  1091.                     new PoissonSeriesParser(17).
  1092.                         withFirstDelaunay(4).
  1093.                         withFirstPlanetary(9).
  1094.                         withSinCos(0, 2, microAS, 3, microAS).
  1095.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1096.             final PoissonSeries gstSeries = load(GST_SERIES, in -> gstParser.parse(in, GST_SERIES));
  1097.             final PoissonSeries.CompiledSeries psiGstSeries =
  1098.                     PoissonSeries.compile(psiLuniSolarSeries, psiPlanetarySeries, gstSeries);

  1100.             // ERA function
  1101.             final TimeScalarFunction era =
  1102.                     getEarthOrientationAngleFunction(ut1, timeScales.getTAI());

  1104.             return new TimeScalarFunction() {

  1106.                 /** {@inheritDoc} */
  1107.                 @Override
  1108.                 public double value(final AbsoluteDate date) {

  1110.                     // evaluate equation of origins
  1111.                     final BodiesElements elements = arguments.evaluateAll(date);
  1112.                     final double[] angles = psiGstSeries.value(elements);
  1113.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  1114.                     final double deltaPsi = angles[0] + angles[1] + ddPsi;
  1115.                     final double epsilonA = epsilon.value(date);

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

  1121.                 }

  1123.                 /** {@inheritDoc} */
  1124.                 @Override
  1125.                 public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  1127.                     // evaluate equation of origins
  1128.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1129.                     final T[] angles = psiGstSeries.value(elements);
  1130.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  1131.                     final T deltaPsi = angles[0].add(angles[1]).add(ddPsi);
  1132.                     final T epsilonA = epsilon.value(date);

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

  1138.                 }

  1140.             };

  1142.         }

  1144.         /** {@inheritDoc} */
  1145.         @Override
  1146.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

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

  1156.             // set up Poisson series
  1157.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1158.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
  1159.                     withOptionalColumn(1).
  1160.                     withGamma(2).
  1161.                     withFirstDelaunay(3);
  1162.             final double microS = 1.0e-6;
  1163.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  1164.                             withOptionalColumn(1).
  1165.                             withGamma(2).
  1166.                             withFirstDelaunay(3);
  1167.             final PoissonSeries xSeries =
  1168.                             load(TIDAL_CORRECTION_XP_YP_SERIES, in -> xyParser.
  1169.                                                                       withSinCos(0, 10, microAS, 11, microAS).
  1170.                                                                       parse(in, TIDAL_CORRECTION_XP_YP_SERIES));
  1171.             final PoissonSeries ySeries =
  1172.                             load(TIDAL_CORRECTION_XP_YP_SERIES, in -> xyParser.
  1173.                                                                       withSinCos(0, 12, microAS, 13, microAS).
  1174.                                                                       parse(in, TIDAL_CORRECTION_XP_YP_SERIES));

  1176.             final PoissonSeries ut1Series =
  1177.                             load(TIDAL_CORRECTION_UT1_SERIES, in -> ut1Parser.
  1178.                                                                     withSinCos(0, 10, microS, 11, microS).
  1179.                                                                     parse(in, TIDAL_CORRECTION_UT1_SERIES));

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

  1183.         }

  1185.         /** {@inheritDoc} */
  1186.         public LoveNumbers getLoveNumbers() {
  1187.             return loadLoveNumbers(LOVE_NUMBERS);
  1188.         }

  1190.         /** {@inheritDoc} */
  1191.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  1192.                                                                      final TimeScales timeScales) {

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

  1197.             // set up Poisson series
  1198.             final PoissonSeriesParser k20Parser =
  1199.                     new PoissonSeriesParser(18).
  1200.                         withOptionalColumn(1).
  1201.                         withDoodson(4, 2).
  1202.                         withFirstDelaunay(10);
  1203.             final PoissonSeriesParser k21Parser =
  1204.                     new PoissonSeriesParser(18).
  1205.                         withOptionalColumn(1).
  1206.                         withDoodson(4, 3).
  1207.                         withFirstDelaunay(10);
  1208.             final PoissonSeriesParser k22Parser =
  1209.                     new PoissonSeriesParser(16).
  1210.                         withOptionalColumn(1).
  1211.                         withDoodson(4, 2).
  1212.                         withFirstDelaunay(10);

  1214.             final double pico = 1.0e-12;
  1215.             final PoissonSeries c20Series =
  1216.                             load(K20_FREQUENCY_DEPENDENCE, in -> k20Parser.
  1217.                                                                  withSinCos(0, 18, -pico, 16, pico).
  1218.                                                                  parse(in, K20_FREQUENCY_DEPENDENCE));
  1219.             final PoissonSeries c21Series =
  1220.                             load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
  1221.                                                                  withSinCos(0, 17, pico, 18, pico).
  1222.                                                                  parse(in, K21_FREQUENCY_DEPENDENCE));
  1223.             final PoissonSeries s21Series =
  1224.                             load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
  1225.                                                                  withSinCos(0, 18, -pico, 17, pico).
  1226.                                                                  parse(in, K21_FREQUENCY_DEPENDENCE));
  1227.             final PoissonSeries c22Series =
  1228.                             load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
  1229.                                                                  withSinCos(0, -1, pico, 16, pico).
  1230.                                                                  parse(in, K22_FREQUENCY_DEPENDENCE));
  1231.             final PoissonSeries s22Series =
  1232.                             load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
  1233.                                                                  withSinCos(0, 16, -pico, -1, pico).
  1234.                                                                  parse(in, K22_FREQUENCY_DEPENDENCE));

  1236.             return new TideFrequencyDependenceFunction(arguments,
  1237.                                                        c20Series,
  1238.                                                        c21Series, s21Series,
  1239.                                                        c22Series, s22Series);

  1241.         }

  1243.         /** {@inheritDoc} */
  1244.         @Override
  1245.         public double getPermanentTide() {
  1246.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1247.         }

  1249.         /** {@inheritDoc} */
  1250.         @Override
  1251.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  1253.             // annual pole from ftp://tai.bipm.org/iers/conv2003/chapter7/annual.pole
  1254.             final TimeScale utc = eopHistory.getTimeScales().getUTC();
  1255.             final SimpleTimeStampedTableParser.RowConverter<MeanPole> converter =
  1256.                 new SimpleTimeStampedTableParser.RowConverter<MeanPole>() {
  1257.                     /** {@inheritDoc} */
  1258.                     @Override
  1259.                     public MeanPole convert(final double[] rawFields) {
  1260.                         return new MeanPole(new AbsoluteDate((int) rawFields[0], 1, 1, utc),
  1261.                                             rawFields[1] * Constants.ARC_SECONDS_TO_RADIANS,
  1262.                                             rawFields[2] * Constants.ARC_SECONDS_TO_RADIANS);
  1263.                     }
  1264.                 };
  1265.             final SimpleTimeStampedTableParser<MeanPole> parser =
  1266.                     new SimpleTimeStampedTableParser<>(3, converter);
  1267.             final List<MeanPole> annualPoleList = load(ANNUAL_POLE, in -> parser.parse(in, ANNUAL_POLE));
  1268.             final AbsoluteDate firstAnnualPoleDate = annualPoleList.get(0).getDate();
  1269.             final AbsoluteDate lastAnnualPoleDate  = annualPoleList.get(annualPoleList.size() - 1).getDate();
  1270.             final ImmutableTimeStampedCache<MeanPole> annualCache =
  1271.                     new ImmutableTimeStampedCache<>(2, annualPoleList);

  1273.             // polynomial extension from IERS 2003, section 7.1.4, equations 23a and 23b
  1274.             final double xp0    = 0.054   * Constants.ARC_SECONDS_TO_RADIANS;
  1275.             final double xp0Dot = 0.00083 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
  1276.             final double yp0    = 0.357   * Constants.ARC_SECONDS_TO_RADIANS;
  1277.             final double yp0Dot = 0.00395 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;

  1279.             // constants from IERS 2003, section 6.2
  1280.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1281.             final double ratio        =  0.0115;

  1283.             return new TimeVectorFunction() {

  1285.                 /** {@inheritDoc} */
  1286.                 @Override
  1287.                 public double[] value(final AbsoluteDate date) {

  1289.                     // we can't compute anything before the range covered by the annual pole file
  1290.                     if (date.compareTo(firstAnnualPoleDate) <= 0) {
  1291.                         return new double[] {
  1292.                             0.0, 0.0
  1293.                         };
  1294.                     }

  1296.                     // evaluate mean pole
  1297.                     double meanPoleX = 0;
  1298.                     double meanPoleY = 0;
  1299.                     if (date.compareTo(lastAnnualPoleDate) <= 0) {
  1300.                         // we are within the range covered by the annual pole file,
  1301.                         // we interpolate within it
  1302.                         try {
  1303.                             final HermiteInterpolator interpolator = new HermiteInterpolator();
  1304.                             annualCache.getNeighbors(date).forEach(neighbor ->
  1305.                                 interpolator.addSamplePoint(neighbor.getDate().durationFrom(date),
  1306.                                                             new double[] {
  1307.                                                                 neighbor.getX(), neighbor.getY()
  1308.                                                             }));
  1309.                             final double[] interpolated = interpolator.value(0);
  1310.                             meanPoleX = interpolated[0];
  1311.                             meanPoleY = interpolated[1];
  1312.                         } catch (TimeStampedCacheException tsce) {
  1313.                             // this should never happen
  1314.                             throw new OrekitInternalError(tsce);
  1315.                         }
  1316.                     } else {

  1318.                         // we are after the range covered by the annual pole file,
  1319.                         // we use the polynomial extension
  1320.                         final double t = date.durationFrom(
  1321.                                 eopHistory.getTimeScales().getJ2000Epoch());
  1322.                         meanPoleX = xp0 + t * xp0Dot;
  1323.                         meanPoleY = yp0 + t * yp0Dot;

  1325.                     }

  1327.                     // evaluate wobble variables
  1328.                     final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1329.                     final double m1 = correction.getXp() - meanPoleX;
  1330.                     final double m2 = meanPoleY - correction.getYp();

  1332.                     return new double[] {
  1333.                         // the following correspond to the equations published in IERS 2003 conventions,
  1334.                         // section 6.2 page 65. In the publication, the equations read:
  1335.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1336.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1337.                         // However, it seems there are sign errors in these equations, which have
  1338.                         // been fixed in IERS 2010 conventions, section 6.4 page 94. In these newer
  1339.                         // publication, the equations read:
  1340.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1341.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1342.                         // the newer equations seem more consistent with the premises as the
  1343.                         // deformation due to the centrifugal potential has the form:
  1344.                         // −Ω²r²/2 sin 2θ Re [k₂(m₁ − im₂) exp(iλ)] where k₂ is the complex
  1345.                         // number 0.3077 + 0.0036i, so the real part in the previous equation is:
  1346.                         // A[Re(k₂) m₁ + Im(k₂) m₂)] cos λ + A[Re(k₂) m₂ - Im(k₂) m₁] sin λ
  1347.                         // identifying this with ∆C₂₁ cos λ + ∆S₂₁ sin λ we get:
  1348.                         // ∆C₂₁ = A Re(k₂) [m₁ + Im(k₂)/Re(k₂) m₂)]
  1349.                         // ∆S₂₁ = A Re(k₂) [m₂ - Im(k₂)/Re(k₂) m₁)]
  1350.                         // and Im(k₂)/Re(k₂) is very close to +0.0115
  1351.                         // As the equation as written in the IERS 2003 conventions are used in
  1352.                         // legacy systems, we have reproduced this alleged error here (and fixed it in
  1353.                         // the IERS 2010 conventions below) for validation purposes. We don't recommend
  1354.                         // using the IERS 2003 conventions for solid pole tide computation other than
  1355.                         // for validation or reproducibility of legacy applications behavior.
  1356.                         // As solid pole tide is small and as the sign change is on the smallest coefficient,
  1357.                         // the effect is quite small. A test case on a propagated orbit showed a position change
  1358.                         // slightly below 0.4m after a 30 days propagation on a Low Earth Orbit
  1359.                         globalFactor * (m1 - ratio * m2),
  1360.                         globalFactor * (m2 + ratio * m1),
  1361.                     };

  1363.                 }

  1365.                 /** {@inheritDoc} */
  1366.                 @Override
  1367.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

  1376.                     // evaluate mean pole
  1377.                     T meanPoleX = date.getField().getZero();
  1378.                     T meanPoleY = date.getField().getZero();
  1379.                     if (aDate.compareTo(lastAnnualPoleDate) <= 0) {
  1380.                         // we are within the range covered by the annual pole file,
  1381.                         // we interpolate within it
  1382.                         try {
  1383.                             final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
  1384.                             final T[] y = MathArrays.buildArray(date.getField(), 2);
  1385.                             final T zero = date.getField().getZero();
  1386.                             final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
  1387.                                                                                                        // for example removing derivatives
  1388.                                                                                                        // if T was DerivativeStructure
  1389.                             annualCache.getNeighbors(aDate).forEach(neighbor -> {
  1390.                                 y[0] = zero.newInstance(neighbor.getX());
  1391.                                 y[1] = zero.newInstance(neighbor.getY());
  1392.                                 interpolator.addSamplePoint(central.durationFrom(neighbor.getDate()).negate(), y);
  1393.                             });
  1394.                             final T[] interpolated = interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
  1395.                             meanPoleX = interpolated[0];
  1396.                             meanPoleY = interpolated[1];
  1397.                         } catch (TimeStampedCacheException tsce) {
  1398.                             // this should never happen
  1399.                             throw new OrekitInternalError(tsce);
  1400.                         }
  1401.                     } else {

  1403.                         // we are after the range covered by the annual pole file,
  1404.                         // we use the polynomial extension
  1405.                         final T t = date.durationFrom(
  1406.                                 eopHistory.getTimeScales().getJ2000Epoch());
  1407.                         meanPoleX = t.multiply(xp0Dot).add(xp0);
  1408.                         meanPoleY = t.multiply(yp0Dot).add(yp0);

  1410.                     }

  1412.                     // evaluate wobble variables
  1413.                     final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1414.                     final T m1 = correction.getXp().subtract(meanPoleX);
  1415.                     final T m2 = meanPoleY.subtract(correction.getYp());

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

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

  1448.                     return a;

  1450.                 }

  1452.             };

  1454.         }

  1456.         /** {@inheritDoc} */
  1457.         @Override
  1458.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {

  1460.             return new TimeVectorFunction() {

  1462.                 /** {@inheritDoc} */
  1463.                 @Override
  1464.                 public double[] value(final AbsoluteDate date) {
  1465.                     // there are no model for ocean pole tide prior to conventions 2010
  1466.                     return new double[] {
  1467.                         0.0, 0.0
  1468.                     };
  1469.                 }

  1471.                 /** {@inheritDoc} */
  1472.                 @Override
  1473.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1474.                     // there are no model for ocean pole tide prior to conventions 2010
  1475.                     return MathArrays.buildArray(date.getField(), 2);
  1476.                 }

  1478.             };
  1479.         }

  1481.         /** {@inheritDoc} */
  1482.         @Override
  1483.         public double[] getNominalTidalDisplacement() {
  1484.             return new double[] {
  1485.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  1486.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  1487.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  1488.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  1489.                 // H₀
  1490.                 -0.31460
  1491.             };
  1492.         }

  1494.         /** {@inheritDoc} */
  1495.         @Override
  1496.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
  1497.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  1498.                                                                   18, 15, 16, 17, 18);
  1499.         }

  1501.         /** {@inheritDoc} */
  1502.         @Override
  1503.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  1504.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  1505.                                                                 18, 15, 16, 17, 18);
  1506.         }

  1508.     },

  1510.     /** Constant for IERS 2010 conventions. */
  1511.     IERS_2010 {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  1558.         /** {@inheritDoc} */
  1559.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
  1560.                                                                  final TimeScales timeScales) {
  1561.             return load(NUTATION_ARGUMENTS, in -> new FundamentalNutationArguments(this, timeScale,
  1562.                                                                                             in, NUTATION_ARGUMENTS,
  1563.                                                                                             timeScales));
  1564.         }

  1566.         /** {@inheritDoc} */
  1567.         @Override
  1568.         public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {

  1570.             // epsilon 0 value from chapter 5, page 56, other terms from equation 5.40 page 65
  1571.             final PolynomialNutation epsilonA =
  1572.                     new PolynomialNutation(84381.406        * Constants.ARC_SECONDS_TO_RADIANS,
  1573.                                              -46.836769     * Constants.ARC_SECONDS_TO_RADIANS,
  1574.                                               -0.0001831    * Constants.ARC_SECONDS_TO_RADIANS,
  1575.                                                0.00200340   * Constants.ARC_SECONDS_TO_RADIANS,
  1576.                                               -0.000000576  * Constants.ARC_SECONDS_TO_RADIANS,
  1577.                                               -0.0000000434 * Constants.ARC_SECONDS_TO_RADIANS);

  1579.             return new TimeScalarFunction() {

  1581.                 /** {@inheritDoc} */
  1582.                 @Override
  1583.                 public double value(final AbsoluteDate date) {
  1584.                     return epsilonA.value(evaluateTC(date, timeScales));
  1585.                 }

  1587.                 /** {@inheritDoc} */
  1588.                 @Override
  1589.                 public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1590.                     return epsilonA.value(evaluateTC(date, timeScales));
  1591.                 }

  1593.             };

  1595.         }

  1597.         /** {@inheritDoc} */
  1598.         @Override
  1599.         public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {

  1601.             // set up nutation arguments
  1602.             final FundamentalNutationArguments arguments =
  1603.                     getNutationArguments(timeScales);

  1605.             // set up Poisson series
  1606.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1607.             final PoissonSeriesParser parser =
  1608.                     new PoissonSeriesParser(17).
  1609.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  1610.                         withFirstDelaunay(4).
  1611.                         withFirstPlanetary(9).
  1612.                         withSinCos(0, 2, microAS, 3, microAS);
  1613.             final PoissonSeries                xSeries = load(X_SERIES, in -> parser.parse(in, X_SERIES));
  1614.             final PoissonSeries                ySeries = load(Y_SERIES, in -> parser.parse(in, Y_SERIES));
  1615.             final PoissonSeries                sSeries = load(S_SERIES, in -> parser.parse(in, S_SERIES));
  1616.             final PoissonSeries.CompiledSeries xys     = PoissonSeries.compile(xSeries, ySeries, sSeries);

  1618.             // create a function evaluating the series
  1619.             return new TimeVectorFunction() {

  1621.                 /** {@inheritDoc} */
  1622.                 @Override
  1623.                 public double[] value(final AbsoluteDate date) {
  1624.                     return xys.value(arguments.evaluateAll(date));
  1625.                 }

  1627.                 /** {@inheritDoc} */
  1628.                 @Override
  1629.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1630.                     return xys.value(arguments.evaluateAll(date));
  1631.                 }

  1633.             };

  1635.         }

  1637.         /** {@inheritDoc} */
  1638.         public LoveNumbers getLoveNumbers() {
  1639.             return loadLoveNumbers(LOVE_NUMBERS);
  1640.         }

  1642.         /** {@inheritDoc} */
  1643.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  1644.                                                                      final TimeScales timeScales) {

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

  1649.             // set up Poisson series
  1650.             final PoissonSeriesParser k20Parser =
  1651.                     new PoissonSeriesParser(18).
  1652.                         withOptionalColumn(1).
  1653.                         withDoodson(4, 2).
  1654.                         withFirstDelaunay(10);
  1655.             final PoissonSeriesParser k21Parser =
  1656.                     new PoissonSeriesParser(18).
  1657.                         withOptionalColumn(1).
  1658.                         withDoodson(4, 3).
  1659.                         withFirstDelaunay(10);
  1660.             final PoissonSeriesParser k22Parser =
  1661.                     new PoissonSeriesParser(16).
  1662.                         withOptionalColumn(1).
  1663.                         withDoodson(4, 2).
  1664.                         withFirstDelaunay(10);

  1666.             final double pico = 1.0e-12;
  1667.             final PoissonSeries c20Series = load(K20_FREQUENCY_DEPENDENCE, in -> k20Parser.
  1668.                                                                                  withSinCos(0, 18, -pico, 16, pico).
  1669.                                                                                  parse(in, K20_FREQUENCY_DEPENDENCE));
  1670.             final PoissonSeries c21Series = load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
  1671.                                                                                  withSinCos(0, 17, pico, 18, pico).
  1672.                                                                                  parse(in, K21_FREQUENCY_DEPENDENCE));
  1673.             final PoissonSeries s21Series = load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
  1674.                                                                                  withSinCos(0, 18, -pico, 17, pico).
  1675.                                                                                  parse(in, K21_FREQUENCY_DEPENDENCE));
  1676.             final PoissonSeries c22Series = load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
  1677.                                                                                  withSinCos(0, -1, pico, 16, pico).
  1678.                                                                                  parse(in, K22_FREQUENCY_DEPENDENCE));
  1679.             final PoissonSeries s22Series = load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
  1680.                                                                                  withSinCos(0, 16, -pico, -1, pico).
  1681.                                                                                  parse(in, K22_FREQUENCY_DEPENDENCE));

  1683.             return new TideFrequencyDependenceFunction(arguments,
  1684.                                                        c20Series,
  1685.                                                        c21Series, s21Series,
  1686.                                                        c22Series, s22Series);

  1688.         }

  1690.         /** {@inheritDoc} */
  1691.         @Override
  1692.         public double getPermanentTide() {
  1693.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1694.         }

  1696.         /** Compute pole wobble variables m₁ and m₂.
  1697.          * @param date current date
  1698.          * @param eopHistory EOP history
  1699.          * @return array containing m₁ and m₂
  1700.          */
  1701.         private double[] computePoleWobble(final AbsoluteDate date, final EOPHistory eopHistory) {

  1703.             // polynomial model from IERS 2010, table 7.7
  1704.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1705.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1706.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1707.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1708.             final AbsoluteDate changeDate =
  1709.                     new AbsoluteDate(2010, 1, 1, eopHistory.getTimeScales().getTT());

  1711.             // evaluate mean pole
  1712.             final double[] xPolynomial;
  1713.             final double[] yPolynomial;
  1714.             if (date.compareTo(changeDate) <= 0) {
  1715.                 xPolynomial = new double[] {
  1716.                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
  1717.                 };
  1718.                 yPolynomial = new double[] {
  1719.                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
  1720.                 };
  1721.             } else {
  1722.                 xPolynomial = new double[] {
  1723.                     23.513 * f0, 7.6141 * f1
  1724.                 };
  1725.                 yPolynomial = new double[] {
  1726.                     358.891 * f0,  -0.6287 * f1
  1727.                 };
  1728.             }
  1729.             double meanPoleX = 0;
  1730.             double meanPoleY = 0;
  1731.             final double t = date.durationFrom(
  1732.                     eopHistory.getTimeScales().getJ2000Epoch());
  1733.             for (int i = xPolynomial.length - 1; i >= 0; --i) {
  1734.                 meanPoleX = meanPoleX * t + xPolynomial[i];
  1735.             }
  1736.             for (int i = yPolynomial.length - 1; i >= 0; --i) {
  1737.                 meanPoleY = meanPoleY * t + yPolynomial[i];
  1738.             }

  1740.             // evaluate wobble variables
  1741.             final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1742.             final double m1 = correction.getXp() - meanPoleX;
  1743.             final double m2 = meanPoleY - correction.getYp();

  1745.             return new double[] {
  1746.                 m1, m2
  1747.             };

  1749.         }

  1751.         /** Compute pole wobble variables m₁ and m₂.
  1752.          * @param date current date
  1753.          * @param <T> type of the field elements
  1754.          * @param eopHistory EOP history
  1755.          * @return array containing m₁ and m₂
  1756.          */
  1757.         private <T extends CalculusFieldElement<T>> T[] computePoleWobble(final FieldAbsoluteDate<T> date, final EOPHistory eopHistory) {

  1759.             // polynomial model from IERS 2010, table 7.7
  1760.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1761.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1762.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1763.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1764.             final AbsoluteDate changeDate =
  1765.                     new AbsoluteDate(2010, 1, 1, eopHistory.getTimeScales().getTT());

  1767.             // evaluate mean pole
  1768.             final double[] xPolynomial;
  1769.             final double[] yPolynomial;
  1770.             if (date.toAbsoluteDate().compareTo(changeDate) <= 0) {
  1771.                 xPolynomial = new double[] {
  1772.                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
  1773.                 };
  1774.                 yPolynomial = new double[] {
  1775.                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
  1776.                 };
  1777.             } else {
  1778.                 xPolynomial = new double[] {
  1779.                     23.513 * f0, 7.6141 * f1
  1780.                 };
  1781.                 yPolynomial = new double[] {
  1782.                     358.891 * f0,  -0.6287 * f1
  1783.                 };
  1784.             }
  1785.             T meanPoleX = date.getField().getZero();
  1786.             T meanPoleY = date.getField().getZero();
  1787.             final T t = date.durationFrom(
  1788.                     eopHistory.getTimeScales().getJ2000Epoch());
  1789.             for (int i = xPolynomial.length - 1; i >= 0; --i) {
  1790.                 meanPoleX = meanPoleX.multiply(t).add(xPolynomial[i]);
  1791.             }
  1792.             for (int i = yPolynomial.length - 1; i >= 0; --i) {
  1793.                 meanPoleY = meanPoleY.multiply(t).add(yPolynomial[i]);
  1794.             }

  1796.             // evaluate wobble variables
  1797.             final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1798.             final T[] m = MathArrays.buildArray(date.getField(), 2);
  1799.             m[0] = correction.getXp().subtract(meanPoleX);
  1800.             m[1] = meanPoleY.subtract(correction.getYp());

  1802.             return m;

  1804.         }

  1806.         /** {@inheritDoc} */
  1807.         @Override
  1808.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  1810.             // constants from IERS 2010, section 6.4
  1811.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1812.             final double ratio        =  0.0115;

  1814.             return new TimeVectorFunction() {

  1816.                 /** {@inheritDoc} */
  1817.                 @Override
  1818.                 public double[] value(final AbsoluteDate date) {

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

  1823.                     return new double[] {
  1824.                         // the following correspond to the equations published in IERS 2010 conventions,
  1825.                         // section 6.4 page 94. The equations read:
  1826.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1827.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1828.                         // These equations seem to fix what was probably a sign error in IERS 2003
  1829.                         // conventions section 6.2 page 65. In this older publication, the equations read:
  1830.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1831.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1832.                         globalFactor * (wobbleM[0] + ratio * wobbleM[1]),
  1833.                         globalFactor * (wobbleM[1] - ratio * wobbleM[0])
  1834.                     };

  1836.                 }

  1838.                 /** {@inheritDoc} */
  1839.                 @Override
  1840.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

  1847.                     // the following correspond to the equations published in IERS 2010 conventions,
  1848.                     // section 6.4 page 94. The equations read:
  1849.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1850.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1851.                     // These equations seem to fix what was probably a sign error in IERS 2003
  1852.                     // conventions section 6.2 page 65. In this older publication, the equations read:
  1853.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1854.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1855.                     a[0] = wobbleM[0].add(wobbleM[1].multiply( ratio)).multiply(globalFactor);
  1856.                     a[1] = wobbleM[1].add(wobbleM[0].multiply(-ratio)).multiply(globalFactor);

  1858.                     return a;

  1860.                 }

  1862.             };

  1864.         }

  1866.         /** {@inheritDoc} */
  1867.         @Override
  1868.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {

  1870.             return new TimeVectorFunction() {

  1872.                 /** {@inheritDoc} */
  1873.                 @Override
  1874.                 public double[] value(final AbsoluteDate date) {

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

  1879.                     return new double[] {
  1880.                         // the following correspond to the equations published in IERS 2010 conventions,
  1881.                         // section 6.4 page 94 equation 6.24:
  1882.                         // ∆C₂₁ = −2.1778 × 10⁻¹⁰ (m₁ − 0.01724m₂)
  1883.                         // ∆S₂₁ = −1.7232 × 10⁻¹⁰ (m₂ − 0.03365m₁)
  1884.                         -2.1778e-10 * (wobbleM[0] - 0.01724 * wobbleM[1]) / Constants.ARC_SECONDS_TO_RADIANS,
  1885.                         -1.7232e-10 * (wobbleM[1] - 0.03365 * wobbleM[0]) / Constants.ARC_SECONDS_TO_RADIANS
  1886.                     };

  1888.                 }

  1890.                 /** {@inheritDoc} */
  1891.                 @Override
  1892.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

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

  1906.                     return v;

  1908.                 }

  1910.             };

  1912.         }

  1914.         /** {@inheritDoc} */
  1915.         @Override
  1916.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

  1918.             // set up the conventional polynomials
  1919.             // the following values are from equation 5.40 in IERS 2010 conventions
  1920.             final PolynomialNutation psiA =
  1921.                     new PolynomialNutation(   0.0,
  1922.                                            5038.481507     * Constants.ARC_SECONDS_TO_RADIANS,
  1923.                                              -1.0790069    * Constants.ARC_SECONDS_TO_RADIANS,
  1924.                                              -0.00114045   * Constants.ARC_SECONDS_TO_RADIANS,
  1925.                                               0.000132851  * Constants.ARC_SECONDS_TO_RADIANS,
  1926.                                              -0.0000000951 * Constants.ARC_SECONDS_TO_RADIANS);
  1927.             final PolynomialNutation omegaA =
  1928.                     new PolynomialNutation(getMeanObliquityFunction(timeScales)
  1929.                             .value(getNutationReferenceEpoch(timeScales)),
  1930.                                            -0.025754     * Constants.ARC_SECONDS_TO_RADIANS,
  1931.                                             0.0512623    * Constants.ARC_SECONDS_TO_RADIANS,
  1932.                                            -0.00772503   * Constants.ARC_SECONDS_TO_RADIANS,
  1933.                                            -0.000000467  * Constants.ARC_SECONDS_TO_RADIANS,
  1934.                                             0.0000003337 * Constants.ARC_SECONDS_TO_RADIANS);
  1935.             final PolynomialNutation chiA =
  1936.                     new PolynomialNutation( 0.0,
  1937.                                            10.556403     * Constants.ARC_SECONDS_TO_RADIANS,
  1938.                                            -2.3814292    * Constants.ARC_SECONDS_TO_RADIANS,
  1939.                                            -0.00121197   * Constants.ARC_SECONDS_TO_RADIANS,
  1940.                                             0.000170663  * Constants.ARC_SECONDS_TO_RADIANS,
  1941.                                            -0.0000000560 * Constants.ARC_SECONDS_TO_RADIANS);

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

  1945.         }

  1947.         /** {@inheritDoc} */
  1948.         @Override
  1949.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  1951.             // set up nutation arguments
  1952.             final FundamentalNutationArguments arguments =
  1953.                     getNutationArguments(timeScales);

  1955.             // set up Poisson series
  1956.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1957.             final PoissonSeriesParser parser =
  1958.                     new PoissonSeriesParser(17).
  1959.                         withFirstDelaunay(4).
  1960.                         withFirstPlanetary(9).
  1961.                         withSinCos(0, 2, microAS, 3, microAS);
  1962.             final PoissonSeries psiSeries     = load(PSI_SERIES,     in -> parser.parse(in, PSI_SERIES));
  1963.             final PoissonSeries epsilonSeries = load(EPSILON_SERIES, in -> parser.parse(in, EPSILON_SERIES));
  1964.             final PoissonSeries.CompiledSeries psiEpsilonSeries = PoissonSeries.compile(psiSeries, epsilonSeries);

  1966.             return new TimeVectorFunction() {

  1968.                 /** {@inheritDoc} */
  1969.                 @Override
  1970.                 public double[] value(final AbsoluteDate date) {
  1971.                     final BodiesElements elements = arguments.evaluateAll(date);
  1972.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  1973.                     return new double[] {
  1974.                         psiEpsilon[0], psiEpsilon[1],
  1975.                         IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  1976.                     };
  1977.                 }

  1979.                 /** {@inheritDoc} */
  1980.                 @Override
  1981.                 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1982.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1983.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  1984.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  1985.                     result[0] = psiEpsilon[0];
  1986.                     result[1] = psiEpsilon[1];
  1987.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  1988.                     return result;
  1989.                 }

  1991.             };

  1993.         }

  1995.         /** {@inheritDoc} */
  1996.         @Override
  1997.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  1998.                                                   final TimeScales timeScales) {

  2000.             // Earth Rotation Angle
  2001.             final StellarAngleCapitaine era =
  2002.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  2004.             // Polynomial part of the apparent sidereal time series
  2005.             // which is the opposite of Equation of Origins (EO)
  2006.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2007.             final PoissonSeriesParser parser =
  2008.                     new PoissonSeriesParser(17).
  2009.                         withFirstDelaunay(4).
  2010.                         withFirstPlanetary(9).
  2011.                         withSinCos(0, 2, microAS, 3, microAS).
  2012.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  2013.             final PolynomialNutation minusEO = load(GST_SERIES, in ->  parser.parse(in, GST_SERIES).getPolynomial());

  2015.             // create a function evaluating the series
  2016.             return new TimeScalarFunction() {

  2018.                 /** {@inheritDoc} */
  2019.                 @Override
  2020.                 public double value(final AbsoluteDate date) {
  2021.                     return era.value(date) + minusEO.value(evaluateTC(date, timeScales));
  2022.                 }

  2024.                 /** {@inheritDoc} */
  2025.                 @Override
  2026.                 public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  2027.                     return era.value(date).add(minusEO.value(evaluateTC(date, timeScales)));
  2028.                 }

  2030.             };

  2032.         }

  2034.         /** {@inheritDoc} */
  2035.         @Override
  2036.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  2037.                                                       final TimeScales timeScales) {

  2039.             // Earth Rotation Angle
  2040.             final StellarAngleCapitaine era =
  2041.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

  2043.             // Polynomial part of the apparent sidereal time series
  2044.             // which is the opposite of Equation of Origins (EO)
  2045.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2046.             final PoissonSeriesParser parser =
  2047.                     new PoissonSeriesParser(17).
  2048.                         withFirstDelaunay(4).
  2049.                         withFirstPlanetary(9).
  2050.                         withSinCos(0, 2, microAS, 3, microAS).
  2051.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  2052.             final PolynomialNutation minusEO = load(GST_SERIES, in -> parser.parse(in, GST_SERIES).getPolynomial());

  2054.             // create a function evaluating the series
  2055.             return new TimeScalarFunction() {

  2057.                 /** {@inheritDoc} */
  2058.                 @Override
  2059.                 public double value(final AbsoluteDate date) {
  2060.                     return era.getRate() + minusEO.derivative(evaluateTC(date, timeScales));
  2061.                 }

  2063.                 /** {@inheritDoc} */
  2064.                 @Override
  2065.                 public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  2066.                     return minusEO.derivative(evaluateTC(date, timeScales)).add(era.getRate());
  2067.                 }

  2069.             };

  2071.         }

  2073.         /** {@inheritDoc} */
  2074.         @Override
  2075.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  2076.                                                   final EOPHistory eopHistory,
  2077.                                                   final TimeScales timeScales) {

  2079.             // set up nutation arguments
  2080.             final FundamentalNutationArguments arguments =
  2081.                     getNutationArguments(timeScales);

  2083.             // mean obliquity function
  2084.             final TimeScalarFunction epsilon = getMeanObliquityFunction(timeScales);

  2086.             // set up Poisson series
  2087.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2088.             final PoissonSeriesParser baseParser =
  2089.                     new PoissonSeriesParser(17).
  2090.                         withFirstDelaunay(4).
  2091.                         withFirstPlanetary(9).
  2092.                         withSinCos(0, 2, microAS, 3, microAS);
  2093.             final PoissonSeriesParser          gstParser    = baseParser.withPolynomialPart('t', Unit.ARC_SECONDS);
  2094.             final PoissonSeries                psiSeries    = load(PSI_SERIES, in -> baseParser.parse(in, PSI_SERIES));
  2095.             final PoissonSeries                gstSeries    = load(GST_SERIES, in -> gstParser.parse(in, GST_SERIES));
  2096.             final PoissonSeries.CompiledSeries psiGstSeries = PoissonSeries.compile(psiSeries, gstSeries);

  2098.             // ERA function
  2099.             final TimeScalarFunction era =
  2100.                     getEarthOrientationAngleFunction(ut1, timeScales.getTAI());

  2102.             return new TimeScalarFunction() {

  2104.                 /** {@inheritDoc} */
  2105.                 @Override
  2106.                 public double value(final AbsoluteDate date) {

  2108.                     // evaluate equation of origins
  2109.                     final BodiesElements elements = arguments.evaluateAll(date);
  2110.                     final double[] angles = psiGstSeries.value(elements);
  2111.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  2112.                     final double deltaPsi = angles[0] + ddPsi;
  2113.                     final double epsilonA = epsilon.value(date);

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

  2119.                 }

  2121.                 /** {@inheritDoc} */
  2122.                 @Override
  2123.                 public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  2125.                     // evaluate equation of origins
  2126.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  2127.                     final T[] angles = psiGstSeries.value(elements);
  2128.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  2129.                     final T deltaPsi = angles[0].add(ddPsi);
  2130.                     final T epsilonA = epsilon.value(date);

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

  2136.                 }

  2138.             };

  2140.         }

  2142.         /** {@inheritDoc} */
  2143.         @Override
  2144.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

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

  2154.             // set up Poisson series
  2155.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2156.             final PoissonSeriesParser xyParser  = new PoissonSeriesParser(13).
  2157.                                                   withOptionalColumn(1).
  2158.                                                   withGamma(2).
  2159.                                                   withFirstDelaunay(3);
  2160.             final double microS = 1.0e-6;
  2161.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  2162.                                                   withOptionalColumn(1).
  2163.                                                   withGamma(2).
  2164.                                                   withFirstDelaunay(3);
  2165.             final PoissonSeries xSeries =
  2166.                             load(TIDAL_CORRECTION_XP_YP_SERIES, xpIn -> xyParser.
  2167.                                                                         withSinCos(0, 10, microAS, 11, microAS).
  2168.                                                                         parse(xpIn, TIDAL_CORRECTION_XP_YP_SERIES));
  2169.             final PoissonSeries ySeries =
  2170.                             load(TIDAL_CORRECTION_XP_YP_SERIES, ypIn -> xyParser.
  2171.                                                                         withSinCos(0, 12, microAS, 13, microAS).
  2172.                                                                         parse(ypIn, TIDAL_CORRECTION_XP_YP_SERIES));
  2173.             final PoissonSeries ut1Series =
  2174.                             load(TIDAL_CORRECTION_UT1_SERIES, ut1In -> ut1Parser.
  2175.                                                                        withSinCos(0, 10, microS, 11, microS).
  2176.                                                                        parse(ut1In, TIDAL_CORRECTION_UT1_SERIES));

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

  2180.         }

  2182.         /** {@inheritDoc} */
  2183.         @Override
  2184.         public double[] getNominalTidalDisplacement() {
  2185.             return new double[] {
  2186.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  2187.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  2188.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  2189.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  2190.                 // H₀
  2191.                 -0.31460
  2192.             };
  2193.         }

  2195.         /** {@inheritDoc} */
  2196.         @Override
  2197.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
  2198.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  2199.                                                                   18, 15, 16, 17, 18);
  2200.         }

  2202.         /** {@inheritDoc} */
  2203.         @Override
  2204.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  2205.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  2206.                                                                 18, 15, 16, 17, 18);
  2207.         }

  2209.     };

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

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

  2217.     /** Get the reference epoch for fundamental nutation arguments.
  2218.      *
  2219.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2220.      *
  2221.      * @return reference epoch for fundamental nutation arguments
  2222.      * @since 6.1
  2223.      * @see #getNutationReferenceEpoch(TimeScales)
  2224.      */
  2225.     @DefaultDataContext
  2226.     public AbsoluteDate getNutationReferenceEpoch() {
  2227.         return getNutationReferenceEpoch(DataContext.getDefault().getTimeScales());
  2228.     }

  2230.     /**
  2231.      * Get the reference epoch for fundamental nutation arguments.
  2232.      *
  2233.      * @param timeScales to use for the reference epoch.
  2234.      * @return reference epoch for fundamental nutation arguments
  2235.      * @since 10.1
  2236.      */
  2237.     public AbsoluteDate getNutationReferenceEpoch(final TimeScales timeScales) {
  2238.         // IERS 1996, IERS 2003 and IERS 2010 use the same J2000.0 reference date
  2239.         return timeScales.getJ2000Epoch();
  2240.     }

  2242.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2243.      *
  2244.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2245.      *
  2246.      * @param date current date
  2247.      * @return date offset in Julian centuries
  2248.      * @since 6.1
  2249.      * @see #evaluateTC(AbsoluteDate, TimeScales)
  2250.      */
  2251.     @DefaultDataContext
  2252.     public double evaluateTC(final AbsoluteDate date) {
  2253.         return evaluateTC(date, DataContext.getDefault().getTimeScales());
  2254.     }

  2256.     /**
  2257.      * Evaluate the date offset between the current date and the {@link
  2258.      * #getNutationReferenceEpoch() reference date}.
  2259.      *
  2260.      * @param date       current date
  2261.      * @param timeScales used in the evaluation.
  2262.      * @return date offset in Julian centuries
  2263.      * @since 10.1
  2264.      */
  2265.     public double evaluateTC(final AbsoluteDate date, final TimeScales timeScales) {
  2266.         return date.durationFrom(getNutationReferenceEpoch(timeScales)) /
  2267.                 Constants.JULIAN_CENTURY;
  2268.     }

  2270.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2271.      *
  2272.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2273.      *
  2274.      * @param date current date
  2275.      * @param <T> type of the field elements
  2276.      * @return date offset in Julian centuries
  2277.      * @since 9.0
  2278.      * @see #evaluateTC(FieldAbsoluteDate, TimeScales)
  2279.      */
  2280.     @DefaultDataContext
  2281.     public <T extends CalculusFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date) {
  2282.         return evaluateTC(date, DataContext.getDefault().getTimeScales());
  2283.     }

  2285.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2286.      * @param <T> type of the field elements
  2287.      * @param date current date
  2288.      * @param timeScales used in the evaluation.
  2289.      * @return date offset in Julian centuries
  2290.      * @since 10.1
  2291.      */
  2292.     public <T extends CalculusFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date,
  2293.                                                         final TimeScales timeScales) {
  2294.         return date.durationFrom(getNutationReferenceEpoch(timeScales))
  2295.                 .divide(Constants.JULIAN_CENTURY);
  2296.     }

  2298.     /**
  2299.      * Get the fundamental nutation arguments. Does not compute GMST based values: gamma,
  2300.      * gammaDot.
  2301.      *
  2302.      * @param timeScales other time scales used in the computation including TAI and TT.
  2303.      * @return fundamental nutation arguments
  2304.      * @see #getNutationArguments(TimeScale, TimeScales)
  2305.      * @since 10.1
  2306.      */
  2307.     protected FundamentalNutationArguments getNutationArguments(
  2308.             final TimeScales timeScales) {

  2310.         return getNutationArguments(null, timeScales);
  2311.     }

  2313.     /** Get the fundamental nutation arguments.
  2314.      *
  2315.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2316.      *
  2317.      * @param timeScale time scale for computing Greenwich Mean Sidereal Time
  2318.      * (typically {@link TimeScales#getUT1(IERSConventions, boolean) UT1})
  2319.      * @return fundamental nutation arguments
  2320.      * @since 6.1
  2321.      * @see #getNutationArguments(TimeScale, TimeScales)
  2322.      * @see #getNutationArguments(TimeScales)
  2323.      */
  2324.     @DefaultDataContext
  2325.     public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale) {
  2326.         return getNutationArguments(timeScale, DataContext.getDefault().getTimeScales());
  2327.     }

  2329.     /**
  2330.      * Get the fundamental nutation arguments.
  2331.      *
  2332.      * @param timeScale time scale for computing Greenwich Mean Sidereal Time (typically
  2333.      *                  {@link TimeScales#getUT1(IERSConventions, boolean) UT1})
  2334.      * @param timeScales other time scales used in the computation including TAI and TT.
  2335.      * @return fundamental nutation arguments
  2336.      * @since 10.1
  2337.      */
  2338.     public abstract FundamentalNutationArguments getNutationArguments(
  2339.             TimeScale timeScale,
  2340.             TimeScales timeScales);

  2342.     /** Get the function computing mean obliquity of the ecliptic.
  2343.      *
  2344.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2345.      *
  2346.      * @return function computing mean obliquity of the ecliptic
  2347.      * @since 6.1
  2348.      * @see #getMeanObliquityFunction(TimeScales)
  2349.      */
  2350.     @DefaultDataContext
  2351.     public TimeScalarFunction getMeanObliquityFunction() {
  2352.         return getMeanObliquityFunction(DataContext.getDefault().getTimeScales());
  2353.     }

  2355.     /**
  2356.      * Get the function computing mean obliquity of the ecliptic.
  2357.      *
  2358.      * @param timeScales used in computing the function.
  2359.      * @return function computing mean obliquity of the ecliptic
  2360.      * @since 10.1
  2361.      */
  2362.     public abstract TimeScalarFunction getMeanObliquityFunction(TimeScales timeScales);

  2364.     /** Get the function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components.
  2365.      * <p>
  2366.      * The returned function computes the two X, Y components of CIP and the S+XY/2 component of the non-rotating CIO.
  2367.      * </p>
  2368.      *
  2369.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2370.      *
  2371.      * @return function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components
  2372.      * @since 6.1
  2373.      * @see #getXYSpXY2Function(TimeScales)
  2374.      */
  2375.     @DefaultDataContext
  2376.     public TimeVectorFunction getXYSpXY2Function() {
  2377.         return getXYSpXY2Function(DataContext.getDefault().getTimeScales());
  2378.     }

  2380.     /**
  2381.      * Get the function computing the Celestial Intermediate Pole and Celestial
  2382.      * Intermediate Origin components.
  2383.      * <p>
  2384.      * The returned function computes the two X, Y components of CIP and the S+XY/2
  2385.      * component of the non-rotating CIO.
  2386.      * </p>
  2387.      *
  2388.      * @param timeScales used to define the function.
  2389.      * @return function computing the Celestial Intermediate Pole and Celestial
  2390.      * Intermediate Origin components
  2391.      * @since 10.1
  2392.      */
  2393.     public abstract TimeVectorFunction getXYSpXY2Function(TimeScales timeScales);

  2395.     /** Get the function computing the raw Earth Orientation Angle.
  2396.      *
  2397.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2398.      *
  2399.      * <p>
  2400.      * The raw angle does not contain any correction. If for example dTU1 correction
  2401.      * due to tidal effect is desired, it must be added afterward by the caller.
  2402.      * The returned value contain the angle as the value and the angular rate as
  2403.      * the first derivative.
  2404.      * </p>
  2405.      * @param ut1 UT1 time scale
  2406.      * @return function computing the rawEarth Orientation Angle, in the non-rotating origin paradigm
  2407.      * @since 6.1
  2408.      * @see #getEarthOrientationAngleFunction(TimeScale, TimeScale)
  2409.      */
  2410.     @DefaultDataContext
  2411.     public TimeScalarFunction getEarthOrientationAngleFunction(final TimeScale ut1) {
  2412.         return getEarthOrientationAngleFunction(ut1,
  2413.                 DataContext.getDefault().getTimeScales().getTAI());
  2414.     }

  2416.     /** Get the function computing the raw Earth Orientation Angle.
  2417.      * <p>
  2418.      * The raw angle does not contain any correction. If for example dTU1 correction
  2419.      * due to tidal effect is desired, it must be added afterward by the caller.
  2420.      * The returned value contain the angle as the value and the angular rate as
  2421.      * the first derivative.
  2422.      * </p>
  2423.      * @param ut1 UT1 time scale
  2424.      * @param tai TAI time scale
  2425.      * @return function computing the rawEarth Orientation Angle, in the non-rotating origin paradigm
  2426.      * @since 10.1
  2427.      */
  2428.     public TimeScalarFunction getEarthOrientationAngleFunction(final TimeScale ut1,
  2429.                                                                final TimeScale tai) {
  2430.         return new StellarAngleCapitaine(ut1, tai);
  2431.     }

  2433.     /** Get the function computing the precession angles.
  2434.      * <p>
  2435.      * The function returned computes the three precession angles
  2436.      * ψ<sub>A</sub> (around Z axis), ω<sub>A</sub> (around X axis)
  2437.      * and χ<sub>A</sub> (around Z axis). The constant angle ε₀
  2438.      * for the fourth rotation (around X axis) can be retrieved by evaluating the
  2439.      * function returned by {@link #getMeanObliquityFunction()} at {@link
  2440.      * #getNutationReferenceEpoch() nutation reference epoch}.
  2441.      * </p>
  2442.      *
  2443.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2444.      *
  2445.      * @return function computing the precession angle
  2446.      * @since 6.1
  2447.      * @see #getPrecessionFunction(TimeScales)
  2448.      */
  2449.     @DefaultDataContext
  2450.     public TimeVectorFunction getPrecessionFunction()
  2451.     {
  2452.         return getPrecessionFunction(DataContext.getDefault().getTimeScales());
  2453.     }

  2455.     /** Get the function computing the precession angles.
  2456.      * <p>
  2457.      * The function returned computes the three precession angles
  2458.      * ψ<sub>A</sub> (around Z axis), ω<sub>A</sub> (around X axis)
  2459.      * and χ<sub>A</sub> (around Z axis). The constant angle ε₀
  2460.      * for the fourth rotation (around X axis) can be retrieved by evaluating the
  2461.      * function returned by {@link #getMeanObliquityFunction()} at {@link
  2462.      * #getNutationReferenceEpoch() nutation reference epoch}.
  2463.      * </p>
  2464.      * @return function computing the precession angle
  2465.      * @since 10.1
  2466.      * @param timeScales used to define the function.
  2467.      */
  2468.     public abstract TimeVectorFunction getPrecessionFunction(TimeScales timeScales);

  2470.     /** Get the function computing the nutation angles.
  2471.      *
  2472.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2473.      *
  2474.      * <p>
  2475.      * The function returned computes the two classical angles ΔΨ and Δε,
  2476.      * and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
  2477.      * resolution C7 (the correction is forced to 0 before this date)
  2478.      * </p>
  2479.      * @return function computing the nutation in longitude ΔΨ and Δε
  2480.      * and the correction of equation of equinoxes
  2481.      * @since 6.1
  2482.      */
  2483.     @DefaultDataContext
  2484.     public TimeVectorFunction getNutationFunction() {
  2485.         return getNutationFunction(DataContext.getDefault().getTimeScales());
  2486.     }

  2488.     /** Get the function computing the nutation angles.
  2489.      * <p>
  2490.      * The function returned computes the two classical angles ΔΨ and Δε,
  2491.      * and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
  2492.      * resolution C7 (the correction is forced to 0 before this date)
  2493.      * </p>
  2494.      * @return function computing the nutation in longitude ΔΨ and Δε
  2495.      * and the correction of equation of equinoxes
  2496.      * @param timeScales used in the computation including TAI and TT.
  2497.      * @since 10.1
  2498.      */
  2499.     public abstract TimeVectorFunction getNutationFunction(TimeScales timeScales);

  2501.     /** Get the function computing Greenwich mean sidereal time, in radians.
  2502.      *
  2503.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2504.      *
  2505.      * @param ut1 UT1 time scale
  2506.      * @return function computing Greenwich mean sidereal time
  2507.      * @since 6.1
  2508.      * @see #getGMSTFunction(TimeScale, TimeScales)
  2509.      */
  2510.     @DefaultDataContext
  2511.     public TimeScalarFunction getGMSTFunction(final TimeScale ut1) {
  2512.         return getGMSTFunction(ut1, DataContext.getDefault().getTimeScales());
  2513.     }

  2515.     /**
  2516.      * Get the function computing Greenwich mean sidereal time, in radians.
  2517.      *
  2518.      * @param ut1 UT1 time scale
  2519.      * @param timeScales other time scales used in the computation including TAI and TT.
  2520.      * @return function computing Greenwich mean sidereal time
  2521.      * @since 10.1
  2522.      */
  2523.     public abstract TimeScalarFunction getGMSTFunction(TimeScale ut1,
  2524.                                                        TimeScales timeScales);

  2526.     /** Get the function computing Greenwich mean sidereal time rate, in radians per second.
  2527.      *
  2528.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2529.      *
  2530.      * @param ut1 UT1 time scale
  2531.      * @return function computing Greenwich mean sidereal time rate
  2532.      * @since 9.0
  2533.      * @see #getGMSTRateFunction(TimeScale, TimeScales)
  2534.      */
  2535.     @DefaultDataContext
  2536.     public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1) {
  2537.         return getGMSTRateFunction(ut1,
  2538.                 DataContext.getDefault().getTimeScales());
  2539.     }

  2541.     /**
  2542.      * Get the function computing Greenwich mean sidereal time rate, in radians per
  2543.      * second.
  2544.      *
  2545.      * @param ut1 UT1 time scale
  2546.      * @param timeScales other time scales used in the computation including TAI and TT.
  2547.      * @return function computing Greenwich mean sidereal time rate
  2548.      * @since 10.1
  2549.      */
  2550.     public abstract TimeScalarFunction getGMSTRateFunction(TimeScale ut1,
  2551.                                                            TimeScales timeScales);

  2553.     /**
  2554.      * Get the function computing Greenwich apparent sidereal time, in radians.
  2555.      *
  2556.      * <p>This method uses the {@link DataContext#getDefault() default data context} if
  2557.      * {@code eopHistory == null}.
  2558.      *
  2559.      * @param ut1        UT1 time scale
  2560.      * @param eopHistory EOP history. If {@code null} then no nutation correction is
  2561.      *                   applied for EOP.
  2562.      * @return function computing Greenwich apparent sidereal time
  2563.      * @since 6.1
  2564.      * @see #getGASTFunction(TimeScale, EOPHistory, TimeScales)
  2565.      */
  2566.     @DefaultDataContext
  2567.     public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  2568.                                               final EOPHistory eopHistory) {
  2569.         final TimeScales timeScales = eopHistory != null ?
  2570.                 eopHistory.getTimeScales() :
  2571.                 DataContext.getDefault().getTimeScales();
  2572.         return getGASTFunction(ut1, eopHistory, timeScales);
  2573.     }

  2575.     /**
  2576.      * Get the function computing Greenwich apparent sidereal time, in radians.
  2577.      *
  2578.      * @param ut1        UT1 time scale
  2579.      * @param eopHistory EOP history. If {@code null} then no nutation correction is
  2580.      *                   applied for EOP.
  2581.      * @param timeScales        TAI time scale.
  2582.      * @return function computing Greenwich apparent sidereal time
  2583.      * @since 10.1
  2584.      */
  2585.     public abstract TimeScalarFunction getGASTFunction(TimeScale ut1,
  2586.                                                        EOPHistory eopHistory,
  2587.                                                        TimeScales timeScales);

  2589.     /** Get the function computing tidal corrections for Earth Orientation Parameters.
  2590.      *
  2591.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2592.      *
  2593.      * @return function computing tidal corrections for Earth Orientation Parameters,
  2594.      * for xp, yp, ut1 and lod respectively
  2595.      * @since 6.1
  2596.      * @see #getEOPTidalCorrection(TimeScales)
  2597.      */
  2598.     @DefaultDataContext
  2599.     public TimeVectorFunction getEOPTidalCorrection() {
  2600.         return getEOPTidalCorrection(DataContext.getDefault().getTimeScales());
  2601.     }

  2603.     /**
  2604.      * Get the function computing tidal corrections for Earth Orientation Parameters.
  2605.      *
  2606.      * @param timeScales used in the computation. The TT and TAI scales are used.
  2607.      * @return function computing tidal corrections for Earth Orientation Parameters, for
  2608.      * xp, yp, ut1 and lod respectively
  2609.      * @since 10.1
  2610.      */
  2611.     public abstract TimeVectorFunction getEOPTidalCorrection(TimeScales timeScales);

  2613.     /** Get the Love numbers.
  2614.      * @return Love numbers
  2615.           * @since 6.1
  2616.      */
  2617.     public abstract LoveNumbers getLoveNumbers();

  2619.     /** Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2620.      *
  2621.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2622.      *
  2623.      * @param ut1 UT1 time scale
  2624.      * @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2625.      * @since 6.1
  2626.      * @see #getTideFrequencyDependenceFunction(TimeScale, TimeScales)
  2627.      */
  2628.     @DefaultDataContext
  2629.     public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1) {
  2630.         return getTideFrequencyDependenceFunction(ut1,
  2631.                 DataContext.getDefault().getTimeScales());
  2632.     }

  2634.     /**
  2635.      * Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂,
  2636.      * ΔS₂₂).
  2637.      *
  2638.      * @param ut1 UT1 time scale
  2639.      * @param timeScales other time scales used in the computation including TAI and TT.
  2640.      * @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂,
  2641.      * ΔS₂₂).
  2642.      * @since 10.1
  2643.      */
  2644.     public abstract TimeVectorFunction getTideFrequencyDependenceFunction(
  2645.             TimeScale ut1,
  2646.             TimeScales timeScales);

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

  2653.     /** Get the function computing solid pole tide (ΔC₂₁, ΔS₂₁).
  2654.      * @param eopHistory EOP history
  2655.      * @return model for solid pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2656.           * @since 6.1
  2657.      */
  2658.     public abstract TimeVectorFunction getSolidPoleTide(EOPHistory eopHistory);

  2660.     /** Get the function computing ocean pole tide (ΔC₂₁, ΔS₂₁).
  2661.      * @param eopHistory EOP history
  2662.      * @return model for ocean pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2663.           * @since 6.1
  2664.      */
  2665.     public abstract TimeVectorFunction getOceanPoleTide(EOPHistory eopHistory);

  2667.     /** Get the nominal values of the displacement numbers.
  2668.      * @return an array containing h⁽⁰⁾, h⁽²⁾, h₃, hI diurnal, hI semi-diurnal,
  2669.      * l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾, l₃, lI diurnal, lI semi-diurnal,
  2670.      * H₀ permanent deformation amplitude
  2671.      * @since 9.1
  2672.      */
  2673.     public abstract double[] getNominalTidalDisplacement();

  2675.     /** Get the correction function for tidal displacement for diurnal tides.
  2676.      * <ul>
  2677.      *  <li>f[0]: radial correction, longitude cosine part</li>
  2678.      *  <li>f[1]: radial correction, longitude sine part</li>
  2679.      *  <li>f[2]: North correction, longitude cosine part</li>
  2680.      *  <li>f[3]: North correction, longitude sine part</li>
  2681.      *  <li>f[4]: East correction, longitude cosine part</li>
  2682.      *  <li>f[5]: East correction, longitude sine part</li>
  2683.      * </ul>
  2684.      * @return correction function for tidal displacement
  2685.      * @since 9.1
  2686.      */
  2687.     public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal();

  2689.     /** Get the correction function for tidal displacement for diurnal tides.
  2690.      * <ul>
  2691.      *  <li>f[0]: radial correction, longitude cosine part</li>
  2692.      *  <li>f[1]: radial correction, longitude sine part</li>
  2693.      *  <li>f[2]: North correction, longitude cosine part</li>
  2694.      *  <li>f[3]: North correction, longitude sine part</li>
  2695.      *  <li>f[4]: East correction, longitude cosine part</li>
  2696.      *  <li>f[5]: East correction, longitude sine part</li>
  2697.      * </ul>
  2698.      * @param tableName name for the diurnal tides table
  2699.      * @param cols total number of columns of the diurnal tides table
  2700.      * @param rIp column holding ∆Rf(ip) in the diurnal tides table, counting from 1
  2701.      * @param rOp column holding ∆Rf(op) in the diurnal tides table, counting from 1
  2702.      * @param tIp column holding ∆Tf(ip) in the diurnal tides table, counting from 1
  2703.      * @param tOp column holding ∆Tf(op) in the diurnal tides table, counting from 1
  2704.      * @return correction function for tidal displacement for diurnal tides
  2705.           * @since 9.1
  2706.      */
  2707.     protected static CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal(final String tableName, final int cols,
  2708.                                                                                    final int rIp, final int rOp,
  2709.                                                                                    final int tIp, final int tOp) {

  2711.         // radial component, missing the sin 2φ factor; this corresponds to:
  2712.         //  - equation 15a in IERS conventions 1996, chapter 7
  2713.         //  - equation 16a in IERS conventions 2003, chapter 7
  2714.         //  - equation 7.12a in IERS conventions 2010, chapter 7
  2715.         final PoissonSeries drCos = load(tableName, in -> new PoissonSeriesParser(cols).
  2716.                                                           withOptionalColumn(1).
  2717.                                                           withDoodson(4, 3).
  2718.                                                           withFirstDelaunay(10).
  2719.                                                           withSinCos(0, rIp, +1.0e-3, rOp, +1.0e-3).
  2720.                                                           parse(in, tableName));
  2721.         final PoissonSeries drSin = load(tableName, in -> new PoissonSeriesParser(cols).
  2722.                                                           withOptionalColumn(1).
  2723.                                                           withDoodson(4, 3).
  2724.                                                           withFirstDelaunay(10).
  2725.                                                           withSinCos(0, rOp, -1.0e-3, rIp, +1.0e-3).
  2726.                                                           parse(in, tableName));

  2728.         // North component, missing the cos 2φ factor; this corresponds to:
  2729.         //  - equation 15b in IERS conventions 1996, chapter 7
  2730.         //  - equation 16b in IERS conventions 2003, chapter 7
  2731.         //  - equation 7.12b in IERS conventions 2010, chapter 7
  2732.         final PoissonSeries dnCos = load(tableName, in -> new PoissonSeriesParser(cols).
  2733.                                                           withOptionalColumn(1).
  2734.                                                           withDoodson(4, 3).
  2735.                                                           withFirstDelaunay(10).
  2736.                                                           withSinCos(0, tIp, +1.0e-3, tOp, +1.0e-3).
  2737.                                                           parse(in, tableName));
  2738.         final PoissonSeries dnSin = load(tableName, in -> new PoissonSeriesParser(cols).
  2739.                                                           withOptionalColumn(1).
  2740.                                                           withDoodson(4, 3).
  2741.                                                           withFirstDelaunay(10).
  2742.                                                           withSinCos(0, tOp, -1.0e-3, tIp, +1.0e-3).
  2743.                                                           parse(in, tableName));

  2745.         // East component, missing the sin φ factor; this corresponds to:
  2746.         //  - equation 15b in IERS conventions 1996, chapter 7
  2747.         //  - equation 16b in IERS conventions 2003, chapter 7
  2748.         //  - equation 7.12b in IERS conventions 2010, chapter 7
  2749.         final PoissonSeries deCos = load(tableName, in -> new PoissonSeriesParser(cols).
  2750.                                                           withOptionalColumn(1).
  2751.                                                           withDoodson(4, 3).
  2752.                                                           withFirstDelaunay(10).
  2753.                                                           withSinCos(0, tOp, -1.0e-3, tIp, +1.0e-3).
  2754.                                                           parse(in, tableName));
  2755.         final PoissonSeries deSin = load(tableName, in -> new PoissonSeriesParser(cols).
  2756.                                                           withOptionalColumn(1).
  2757.                                                           withDoodson(4, 3).
  2758.                                                           withFirstDelaunay(10).
  2759.                                                           withSinCos(0, tIp, -1.0e-3, tOp, -1.0e-3).
  2760.                                                           parse(in, tableName));

  2762.         return PoissonSeries.compile(drCos, drSin, dnCos, dnSin, deCos, deSin);

  2764.     }

  2766.     /** Get the correction function for tidal displacement for zonal tides.
  2767.      * <ul>
  2768.      *  <li>f[0]: radial correction</li>
  2769.      *  <li>f[1]: North correction</li>
  2770.      * </ul>
  2771.      * @return correction function for tidal displacement
  2772.      * @since 9.1
  2773.      */
  2774.     public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionZonal();

  2776.     /** Get the correction function for tidal displacement for zonal tides.
  2777.      * <ul>
  2778.      *  <li>f[0]: radial correction</li>
  2779.      *  <li>f[1]: North correction</li>
  2780.      * </ul>
  2781.      * @param tableName name for the zonal tides table
  2782.      * @param cols total number of columns of the table
  2783.      * @param rIp column holding ∆Rf(ip) in the table, counting from 1
  2784.      * @param rOp column holding ∆Rf(op) in the table, counting from 1
  2785.      * @param tIp column holding ∆Tf(ip) in the table, counting from 1
  2786.      * @param tOp column holding ∆Tf(op) in the table, counting from 1
  2787.      * @return correction function for tidal displacement for zonal tides
  2788.           * @since 9.1
  2789.      */
  2790.     protected static CompiledSeries getTidalDisplacementFrequencyCorrectionZonal(final String tableName, final int cols,
  2791.                                                                                  final int rIp, final int rOp,
  2792.                                                                                  final int tIp, final int tOp) {

  2794.         // radial component, missing the 3⁄2 sin² φ - 1⁄2 factor; this corresponds to:
  2795.         //  - equation 16a in IERS conventions 1996, chapter 7
  2796.         //  - equation 17a in IERS conventions 2003, chapter 7
  2797.         //  - equation 7.13a in IERS conventions 2010, chapter 7
  2798.         final PoissonSeries dr = load(tableName, in -> new PoissonSeriesParser(cols).
  2799.                                                        withOptionalColumn(1).
  2800.                                                        withDoodson(4, 3).
  2801.                                                        withFirstDelaunay(10).
  2802.                                                        withSinCos(0, rOp, +1.0e-3, rIp, +1.0e-3).
  2803.                                                        parse(in, tableName));

  2805.         // North component, missing the sin 2φ factor; this corresponds to:
  2806.         //  - equation 16b in IERS conventions 1996, chapter 7
  2807.         //  - equation 17b in IERS conventions 2003, chapter 7
  2808.         //  - equation 7.13b in IERS conventions 2010, chapter 7
  2809.         final PoissonSeries dn = load(tableName, in -> new PoissonSeriesParser(cols).
  2810.                                                        withOptionalColumn(1).
  2811.                                                        withDoodson(4, 3).
  2812.                                                        withFirstDelaunay(10).
  2813.                                                        withSinCos(0, tOp, +1.0e-3, tIp, +1.0e-3).
  2814.                                                        parse(in, tableName));

  2816.         return PoissonSeries.compile(dr, dn);

  2818.     }

  2820.     /** Interface for functions converting nutation corrections between
  2821.      * δΔψ/δΔε to δX/δY.
  2822.      * <ul>
  2823.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  2824.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  2825.      * </ul>
  2826.      * @since 6.1
  2827.      */
  2828.     public interface NutationCorrectionConverter {

  2830.         /** Convert nutation corrections.
  2831.          * @param date current date
  2832.          * @param ddPsi δΔψ part of the nutation correction
  2833.          * @param ddEpsilon δΔε part of the nutation correction
  2834.          * @return array containing δX and δY
  2835.          */
  2836.         double[] toNonRotating(AbsoluteDate date, double ddPsi, double ddEpsilon);

  2838.         /** Convert nutation corrections.
  2839.          * @param date current date
  2840.          * @param dX δX part of the nutation correction
  2841.          * @param dY δY part of the nutation correction
  2842.          * @return array containing δΔψ and δΔε
  2843.          */
  2844.         double[] toEquinox(AbsoluteDate date, double dX, double dY);

  2846.     }

  2848.     /** Create a function converting nutation corrections between
  2849.      * δX/δY and δΔψ/δΔε.
  2850.      * <ul>
  2851.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  2852.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  2853.      * </ul>
  2854.      *
  2855.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2856.      *
  2857.      * @return a new converter
  2858.      * @since 6.1
  2859.      * @see #getNutationCorrectionConverter(TimeScales)
  2860.      */
  2861.     @DefaultDataContext
  2862.     public NutationCorrectionConverter getNutationCorrectionConverter() {
  2863.         return getNutationCorrectionConverter(DataContext.getDefault().getTimeScales());
  2864.     }

  2866.     /** Create a function converting nutation corrections between
  2867.      * δX/δY and δΔψ/δΔε.
  2868.      * <ul>
  2869.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  2870.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  2871.      * </ul>
  2872.      * @return a new converter
  2873.      * @since 10.1
  2874.      * @param timeScales used to define the conversion.
  2875.      */
  2876.     public NutationCorrectionConverter getNutationCorrectionConverter(
  2877.             final TimeScales timeScales) {

  2879.         // get models parameters
  2880.         final TimeVectorFunction precessionFunction = getPrecessionFunction(timeScales);
  2881.         final TimeScalarFunction epsilonAFunction = getMeanObliquityFunction(timeScales);
  2882.         final AbsoluteDate date0 = getNutationReferenceEpoch(timeScales);
  2883.         final double cosE0 = FastMath.cos(epsilonAFunction.value(date0));

  2885.         return new NutationCorrectionConverter() {

  2887.             /** {@inheritDoc} */
  2888.             @Override
  2889.             public double[] toNonRotating(final AbsoluteDate date,
  2890.                                           final double ddPsi, final double ddEpsilon) {
  2891.                 // compute precession angles psiA, omegaA and chiA
  2892.                 final double[] angles = precessionFunction.value(date);

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

  2898.                 // convert nutation corrections (equation 23/IERS-2003 or 5.25/IERS-2010)
  2899.                 return new double[] {
  2900.                     sinEA * ddPsi + c * ddEpsilon,
  2901.                     ddEpsilon - c * sinEA * ddPsi
  2902.                 };

  2904.             }

  2906.             /** {@inheritDoc} */
  2907.             @Override
  2908.             public double[] toEquinox(final AbsoluteDate date,
  2909.                                       final double dX, final double dY) {
  2910.                 // compute precession angles psiA, omegaA and chiA
  2911.                 final double[] angles   = precessionFunction.value(date);

  2913.                 // conversion coefficients
  2914.                 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
  2915.                 final double c     = angles[0] * cosE0 - angles[2];
  2916.                 final double opC2  = 1 + c * c;

  2918.                 // convert nutation corrections (inverse of equation 23/IERS-2003 or 5.25/IERS-2010)
  2919.                 return new double[] {
  2920.                     (dX - c * dY) / (sinEA * opC2),
  2921.                     (dY + c * dX) / opC2
  2922.                 };
  2923.             }

  2925.         };

  2927.     }

  2929.     /** Load the Love numbers.
  2930.      * @param nameLove name of the Love number resource
  2931.      * @return Love numbers
  2932.      */
  2933.     protected LoveNumbers loadLoveNumbers(final String nameLove) {
  2934.         try {

  2936.             // allocate the three triangular arrays for real, imaginary and time-dependent numbers
  2937.             final double[][] real      = new double[4][];
  2938.             final double[][] imaginary = new double[4][];
  2939.             final double[][] plus      = new double[4][];
  2940.             for (int i = 0; i < real.length; ++i) {
  2941.                 real[i]      = new double[i + 1];
  2942.                 imaginary[i] = new double[i + 1];
  2943.                 plus[i]      = new double[i + 1];
  2944.             }

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

  2948.                 if (stream == null) {
  2949.                     // this should never happen with files embedded within Orekit
  2950.                     throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, nameLove);
  2951.                 }

  2953.                 int lineNumber = 1;
  2954.                 String line = null;
  2955.                 // setup the reader
  2956.                 try (BufferedReader reader = new BufferedReader(new InputStreamReader(stream, StandardCharsets.UTF_8))) {

  2958.                     line = reader.readLine();

  2960.                     // look for the Love numbers
  2961.                     while (line != null) {

  2963.                         line = line.trim();
  2964.                         if (!(line.isEmpty() || line.startsWith("#"))) {
  2965.                             final String[] fields = SEPARATOR.split(line);
  2966.                             if (fields.length != 5) {
  2967.                                 // this should never happen with files embedded within Orekit
  2968.                                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  2969.                                                           lineNumber, nameLove, line);
  2970.                             }
  2971.                             final int n = Integer.parseInt(fields[0]);
  2972.                             final int m = Integer.parseInt(fields[1]);
  2973.                             if (n < 2 || n > 3 || m < 0 || m > n) {
  2974.                                 // this should never happen with files embedded within Orekit
  2975.                                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  2976.                                                           lineNumber, nameLove, line);

  2978.                             }
  2979.                             real[n][m]      = Double.parseDouble(fields[2]);
  2980.                             imaginary[n][m] = Double.parseDouble(fields[3]);
  2981.                             plus[n][m]      = Double.parseDouble(fields[4]);
  2982.                         }

  2984.                         // next line
  2985.                         lineNumber++;
  2986.                         line = reader.readLine();

  2988.                     }
  2989.                 } catch (NumberFormatException nfe) {
  2990.                     // this should never happen with files embedded within Orekit
  2991.                     throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  2992.                                               lineNumber, nameLove, line);
  2993.                 }
  2994.             }

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

  2998.         } catch (IOException ioe) {
  2999.             // this should never happen with files embedded within Orekit
  3000.             throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, nameLove);
  3001.         }
  3002.     }

  3004.     /** Load resources.
  3005.      * @param name name of the resource
  3006.      * @param loader loaader for the resource
  3007.      * @param <T> type of the processed data
  3008.      * @return processed data
  3009.      */
  3010.     private static <T> T load(final String name, final Function<InputStream, T> loader) {
  3011.         try (InputStream is = IERSConventions.class.getResourceAsStream(name)) {
  3012.             return loader.apply(is);
  3013.         } catch (IOException ioe) {
  3014.             // this should never happen with internal streams
  3015.             throw new OrekitException(OrekitMessages.INTERNAL_ERROR, ioe);
  3016.         }
  3017.     }

  3019.     /** Correction to equation of equinoxes.
  3020.      * <p>IAU 1994 resolution C7 added two terms to the equation of equinoxes
  3021.      * taking effect since 1997-02-27 for continuity
  3022.      * </p>
  3023.      */
  3024.     private static class IAU1994ResolutionC7 {

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

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


  3033.         /** Evaluate the correction.
  3034.          * @param arguments Delaunay for nutation
  3035.          * @param tai TAI time scale.
  3036.          * @return correction value (0 before 1997-02-27)
  3037.          */
  3038.         public static double value(final DelaunayArguments arguments,
  3039.                                    final TimeScale tai) {
  3040.             /* Start date for applying Moon corrections to the equation of the equinoxes.
  3041.              * This date corresponds to 1997-02-27T00:00:00 UTC, hence the 30s offset from TAI.
  3042.              */
  3043.             final AbsoluteDate modelStart = new AbsoluteDate(1997, 2, 27, 0, 0, 30, tai);
  3044.             if (arguments.getDate().compareTo(modelStart) >= 0) {

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

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

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

  3055.             } else {
  3056.                 return 0.0;
  3057.             }
  3058.         }

  3060.         /** Evaluate the correction.
  3061.          * @param arguments Delaunay for nutation
  3062.          * @param tai TAI time scale.
  3063.          * @param <T> type of the field elements
  3064.          * @return correction value (0 before 1997-02-27)
  3065.          */
  3066.         public static <T extends CalculusFieldElement<T>> T value(
  3067.                 final FieldDelaunayArguments<T> arguments,
  3068.                 final TimeScale tai) {
  3069.             /* Start date for applying Moon corrections to the equation of the equinoxes.
  3070.              * This date corresponds to 1997-02-27T00:00:00 UTC, hence the 30s offset from TAI.
  3071.              */
  3072.             final AbsoluteDate modelStart = new AbsoluteDate(1997, 2, 27, 0, 0, 30, tai);
  3073.             if (arguments.getDate().toAbsoluteDate().compareTo(modelStart) >= 0) {

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

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

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

  3084.             } else {
  3085.                 return arguments.getDate().getField().getZero();
  3086.             }
  3087.         }

  3089.     }

  3091.     /** Stellar angle model.
  3092.      * <p>
  3093.      * The stellar angle computed here has been defined in the paper "A non-rotating origin on the
  3094.      * instantaneous equator: Definition, properties and use", N. Capitaine, Guinot B., and Souchay J.,
  3095.      * Celestial Mechanics, Volume 39, Issue 3, pp 283-307. It has been proposed as a conventional
  3096.      * conventional relationship between UT1 and Earth rotation in the paper "Definition of the Celestial
  3097.      * Ephemeris origin and of UT1 in the International Celestial Reference Frame", Capitaine, N.,
  3098.      * Guinot, B., and McCarthy, D. D., 2000, Astronomy and Astrophysics, 355(1), pp. 398–405.
  3099.      * </p>
  3100.      * <p>
  3101.      * It is presented simply as stellar angle in IERS conventions 1996 but as since been adopted as
  3102.      * the conventional relationship defining UT1 from ICRF and is called Earth Rotation Angle in
  3103.      * IERS conventions 2003 and 2010.
  3104.      * </p>
  3105.      */
  3106.     private static class StellarAngleCapitaine implements TimeScalarFunction {

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

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

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

  3119.         /** UT1 time scale. */
  3120.         private final TimeScale ut1;

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

  3125.         /** Simple constructor.
  3126.          * @param ut1 UT1 time scale
  3127.          * @param tai TAI time scale
  3128.          */
  3129.         StellarAngleCapitaine(final TimeScale ut1, final TimeScale tai) {
  3130.             this.ut1 = ut1;
  3131.             referenceDate = new AbsoluteDate(
  3132.                     DateComponents.J2000_EPOCH,
  3133.                     TimeComponents.H12,
  3134.                     tai);
  3135.         }

  3137.         /** Get the rotation rate.
  3138.          * @return rotation rate
  3139.          */
  3140.         public double getRate() {
  3141.             return ERA_1A + ERA_1B;
  3142.         }

  3144.         /** {@inheritDoc} */
  3145.         @Override
  3146.         public double value(final AbsoluteDate date) {

  3148.             // split the date offset as a full number of days plus a smaller part
  3149.             final int secondsInDay = 86400;
  3150.             final double dt  = date.durationFrom(referenceDate);
  3151.             final long days  = ((long) dt) / secondsInDay;
  3152.             final double dtA = secondsInDay * days;
  3153.             final double dtB = (dt - dtA) + ut1.offsetFromTAI(date).toDouble();

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

  3157.         }

  3159.         /** {@inheritDoc} */
  3160.         @Override
  3161.         public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  3163.             // split the date offset as a full number of days plus a smaller part
  3164.             final int secondsInDay = 86400;
  3165.             final T dt  = date.durationFrom(referenceDate);
  3166.             final long days  = ((long) dt.getReal()) / secondsInDay;
  3167.             final double dtA = secondsInDay * days;
  3168.             final T dtB = dt.subtract(dtA).add(ut1.offsetFromTAI(date.toAbsoluteDate()).toDouble());

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

  3172.         }

  3174.     }

  3176.     /** Mean pole. */
  3177.     private static class MeanPole implements TimeStamped, Serializable {

  3179.         /** Serializable UID. */
  3180.         private static final long serialVersionUID = 20131028L;

  3182.         /** Date. */
  3183.         private final AbsoluteDate date;

  3185.         /** X coordinate. */
  3186.         private final double x;

  3188.         /** Y coordinate. */
  3189.         private final double y;

  3191.         /** Simple constructor.
  3192.          * @param date date
  3193.          * @param x x coordinate
  3194.          * @param y y coordinate
  3195.          */
  3196.         MeanPole(final AbsoluteDate date, final double x, final double y) {
  3197.             this.date = date;
  3198.             this.x    = x;
  3199.             this.y    = y;
  3200.         }

  3202.         /** {@inheritDoc} */
  3203.         @Override
  3204.         public AbsoluteDate getDate() {
  3205.             return date;
  3206.         }

  3208.         /** Get x coordinate.
  3209.          * @return x coordinate
  3210.          */
  3211.         public double getX() {
  3212.             return x;
  3213.         }

  3215.         /** Get y coordinate.
  3216.          * @return y coordinate
  3217.          */
  3218.         public double getY() {
  3219.             return y;
  3220.         }

  3222.     }

  3224.     /** Local class for precession function. */
  3225.     private class PrecessionFunction implements TimeVectorFunction {

  3227.         /** Polynomial nutation for psiA. */
  3228.         private final PolynomialNutation psiA;

  3230.         /** Polynomial nutation for omegaA. */
  3231.         private final PolynomialNutation omegaA;

  3233.         /** Polynomial nutation for chiA. */
  3234.         private final PolynomialNutation chiA;

  3236.         /** Time scales to use. */
  3237.         private final TimeScales timeScales;

  3239.         /** Simple constructor.
  3240.          * @param psiA polynomial nutation for psiA
  3241.          * @param omegaA polynomial nutation for omegaA
  3242.          * @param chiA polynomial nutation for chiA
  3243.          * @param timeScales used in the computation.
  3244.          */
  3245.         PrecessionFunction(final PolynomialNutation psiA,
  3246.                            final PolynomialNutation omegaA,
  3247.                            final PolynomialNutation chiA,
  3248.                            final TimeScales timeScales) {
  3249.             this.psiA   = psiA;
  3250.             this.omegaA = omegaA;
  3251.             this.chiA   = chiA;
  3252.             this.timeScales = timeScales;
  3253.         }


  3256.         /** {@inheritDoc} */
  3257.         @Override
  3258.         public double[] value(final AbsoluteDate date) {
  3259.             final double tc = evaluateTC(date, timeScales);
  3260.             return new double[] {
  3261.                 psiA.value(tc), omegaA.value(tc), chiA.value(tc)
  3262.             };
  3263.         }

  3265.         /** {@inheritDoc} */
  3266.         @Override
  3267.         public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  3268.             final T[] a = MathArrays.buildArray(date.getField(), 3);
  3269.             final T tc = evaluateTC(date, timeScales);
  3270.             a[0] = psiA.value(tc);
  3271.             a[1] = omegaA.value(tc);
  3272.             a[2] = chiA.value(tc);
  3273.             return a;
  3274.         }

  3276.     }

  3278.     /** Local class for tides frequency function. */
  3279.     private static class TideFrequencyDependenceFunction implements TimeVectorFunction {

  3281.         /** Nutation arguments. */
  3282.         private final FundamentalNutationArguments arguments;

  3284.         /** Correction series. */
  3285.         private final PoissonSeries.CompiledSeries kSeries;

  3287.         /** Simple constructor.
  3288.          * @param arguments nutation arguments
  3289.          * @param c20Series correction series for the C20 term
  3290.          * @param c21Series correction series for the C21 term
  3291.          * @param s21Series correction series for the S21 term
  3292.          * @param c22Series correction series for the C22 term
  3293.          * @param s22Series correction series for the S22 term
  3294.          */
  3295.         TideFrequencyDependenceFunction(final FundamentalNutationArguments arguments,
  3296.                                         final PoissonSeries c20Series,
  3297.                                         final PoissonSeries c21Series,
  3298.                                         final PoissonSeries s21Series,
  3299.                                         final PoissonSeries c22Series,
  3300.                                         final PoissonSeries s22Series) {
  3301.             this.arguments = arguments;
  3302.             this.kSeries   = PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);
  3303.         }


  3306.         /** {@inheritDoc} */
  3307.         @Override
  3308.         public double[] value(final AbsoluteDate date) {
  3309.             return kSeries.value(arguments.evaluateAll(date));
  3310.         }

  3312.         /** {@inheritDoc} */
  3313.         @Override
  3314.         public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  3315.             return kSeries.value(arguments.evaluateAll(date));
  3316.         }

  3318.     }

  3320.     /** Local class for EOP tidal corrections. */
  3321.     private static class EOPTidalCorrection implements TimeVectorFunction {

  3323.         /** Nutation arguments. */
  3324.         private final FundamentalNutationArguments arguments;

  3326.         /** Correction series. */
  3327.         private final PoissonSeries.CompiledSeries correctionSeries;

  3329.         /** Simple constructor.
  3330.          * @param arguments nutation arguments
  3331.          * @param xSeries correction series for the x coordinate
  3332.          * @param ySeries correction series for the y coordinate
  3333.          * @param ut1Series correction series for the UT1
  3334.          */
  3335.         EOPTidalCorrection(final FundamentalNutationArguments arguments,
  3336.                            final PoissonSeries xSeries,
  3337.                            final PoissonSeries ySeries,
  3338.                            final PoissonSeries ut1Series) {
  3339.             this.arguments        = arguments;
  3340.             this.correctionSeries = PoissonSeries.compile(xSeries, ySeries, ut1Series);
  3341.         }

  3343.         /** {@inheritDoc} */
  3344.         @Override
  3345.         public double[] value(final AbsoluteDate date) {
  3346.             final BodiesElements elements = arguments.evaluateAll(date);
  3347.             final double[] correction    = correctionSeries.value(elements);
  3348.             final double[] correctionDot = correctionSeries.derivative(elements);
  3349.             return new double[] {
  3350.                 correction[0],
  3351.                 correction[1],
  3352.                 correction[2],
  3353.                 correctionDot[2] * (-Constants.JULIAN_DAY)
  3354.             };
  3355.         }

  3357.         /** {@inheritDoc} */
  3358.         @Override
  3359.         public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

  3361.             final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  3362.             final T[] correction    = correctionSeries.value(elements);
  3363.             final T[] correctionDot = correctionSeries.derivative(elements);
  3364.             final T[] a = MathArrays.buildArray(date.getField(), 4);
  3365.             a[0] = correction[0];
  3366.             a[1] = correction[1];
  3367.             a[2] = correction[2];
  3368.             a[3] = correctionDot[2].multiply(-Constants.JULIAN_DAY);
  3369.             return a;
  3370.         }

  3372.     }

  3374. }
Created with JaCoCo 0.8.12.202403310830