IERSConventions.java

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

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

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


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

  68.     /** Constant for IERS 1996 conventions. */
  69.     IERS_1996 {

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

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

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

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

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

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

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

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

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

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

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

  104.         /** {@inheritDoc} */
  105.         @Override
  106.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
  107.             throws OrekitException {
  108.             return new FundamentalNutationArguments(this, timeScale,
  109.                                                     getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
  110.         }

  112.         /** {@inheritDoc} */
  113.         @Override
  114.         public TimeScalarFunction getMeanObliquityFunction() throws OrekitException {

  116.             // value from chapter 5, page 22
  117.             final PolynomialNutation epsilonA =
  118.                     new PolynomialNutation(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  119.                                              -46.8150   * Constants.ARC_SECONDS_TO_RADIANS,
  120.                                               -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  121.                                                0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  123.             return new TimeScalarFunction() {

  125.                 /** {@inheritDoc} */
  126.                 @Override
  127.                 public double value(final AbsoluteDate date) {
  128.                     return epsilonA.value(evaluateTC(date));
  129.                 }

  131.                 /** {@inheritDoc} */
  132.                 @Override
  133.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  134.                     return epsilonA.value(evaluateTC(date));
  135.                 }

  137.             };

  139.         }

  141.         /** {@inheritDoc} */
  142.         @Override
  143.         public TimeVectorFunction getXYSpXY2Function()
  144.             throws OrekitException {

  146.             // set up nutation arguments
  147.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  149.             // X = 2004.3109″t - 0.42665″t² - 0.198656″t³ + 0.0000140″t⁴
  150.             //     + 0.00006″t² cos Ω + sin ε0 { Σ [(Ai + Ai' t) sin(ARGUMENT) + Ai'' t cos(ARGUMENT)]}
  151.             //     + 0.00204″t² sin Ω + 0.00016″t² sin 2(F - D + Ω),
  152.             final PolynomialNutation xPolynomial =
  153.                     new PolynomialNutation(0,
  154.                                            2004.3109 * Constants.ARC_SECONDS_TO_RADIANS,
  155.                                            -0.42665  * Constants.ARC_SECONDS_TO_RADIANS,
  156.                                            -0.198656 * Constants.ARC_SECONDS_TO_RADIANS,
  157.                                            0.0000140 * Constants.ARC_SECONDS_TO_RADIANS);

  159.             final double fXCosOm    = 0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
  160.             final double fXSinOm    = 0.00204 * Constants.ARC_SECONDS_TO_RADIANS;
  161.             final double fXSin2FDOm = 0.00016 * Constants.ARC_SECONDS_TO_RADIANS;
  162.             final double sinEps0   = FastMath.sin(getMeanObliquityFunction().value(getNutationReferenceEpoch()));

  164.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  165.             final PoissonSeriesParser baseParser =
  166.                     new PoissonSeriesParser(12).withFirstDelaunay(1);

  168.             final PoissonSeriesParser xParser =
  169.                     baseParser.
  170.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  171.                     withSinCos(1, 8, deciMilliAS,  9, deciMilliAS);
  172.             final PoissonSeries xSum = xParser.parse(getStream(X_Y_SERIES), X_Y_SERIES);

  174.             // Y = -0.00013″ - 22.40992″t² + 0.001836″t³ + 0.0011130″t⁴
  175.             //     + Σ [(Bi + Bi' t) cos(ARGUMENT) + Bi'' t sin(ARGUMENT)]
  176.             //    - 0.00231″t² cos Ω − 0.00014″t² cos 2(F - D + Ω)
  177.             final PolynomialNutation yPolynomial =
  178.                     new PolynomialNutation(-0.00013  * Constants.ARC_SECONDS_TO_RADIANS,
  179.                                            0.0,
  180.                                            -22.40992 * Constants.ARC_SECONDS_TO_RADIANS,
  181.                                            0.001836  * Constants.ARC_SECONDS_TO_RADIANS,
  182.                                            0.0011130 * Constants.ARC_SECONDS_TO_RADIANS);

  184.             final double fYCosOm    = -0.00231 * Constants.ARC_SECONDS_TO_RADIANS;
  185.             final double fYCos2FDOm = -0.00014 * Constants.ARC_SECONDS_TO_RADIANS;

  187.             final PoissonSeriesParser yParser =
  188.                     baseParser.
  189.                     withSinCos(0, -1, deciMilliAS, 10, deciMilliAS).
  190.                     withSinCos(1, 12, deciMilliAS, 11, deciMilliAS);
  191.             final PoissonSeries ySum = yParser.parse(getStream(X_Y_SERIES), X_Y_SERIES);

  193.             final PoissonSeries.CompiledSeries xySum =
  194.                     PoissonSeries.compile(xSum, ySum);

  196.             // s = -XY/2 + 0.00385″t - 0.07259″t³ - 0.00264″ sin Ω - 0.00006″ sin 2Ω
  197.             //     + 0.00074″t² sin Ω + 0.00006″t² sin 2(F - D + Ω)
  198.             final double fST          =  0.00385 * Constants.ARC_SECONDS_TO_RADIANS;
  199.             final double fST3         = -0.07259 * Constants.ARC_SECONDS_TO_RADIANS;
  200.             final double fSSinOm      = -0.00264 * Constants.ARC_SECONDS_TO_RADIANS;
  201.             final double fSSin2Om     = -0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
  202.             final double fST2SinOm    =  0.00074 * Constants.ARC_SECONDS_TO_RADIANS;
  203.             final double fST2Sin2FDOm =  0.00006 * Constants.ARC_SECONDS_TO_RADIANS;

  205.             return new TimeVectorFunction() {

  207.                 /** {@inheritDoc} */
  208.                 @Override
  209.                 public double[] value(final AbsoluteDate date) {

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

  214.                     final double omega     = elements.getOmega();
  215.                     final double f         = elements.getF();
  216.                     final double d         = elements.getD();
  217.                     final double t         = elements.getTC();

  219.                     final double cosOmega  = FastMath.cos(omega);
  220.                     final double sinOmega  = FastMath.sin(omega);
  221.                     final double sin2Omega = FastMath.sin(2 * omega);
  222.                     final double cos2FDOm  = FastMath.cos(2 * (f - d + omega));
  223.                     final double sin2FDOm  = FastMath.sin(2 * (f - d + omega));

  225.                     final double x = xPolynomial.value(t) + sinEps0 * xy[0] +
  226.                             t * t * (fXCosOm * cosOmega + fXSinOm * sinOmega + fXSin2FDOm * cos2FDOm);
  227.                     final double y = yPolynomial.value(t) + xy[1] +
  228.                             t * t * (fYCosOm * cosOmega + fYCos2FDOm * cos2FDOm);
  229.                     final double sPxy2 = fSSinOm * sinOmega + fSSin2Om * sin2Omega +
  230.                             t * (fST + t * (fST2SinOm * sinOmega + fST2Sin2FDOm * sin2FDOm + t * fST3));

  232.                     return new double[] {
  233.                         x, y, sPxy2
  234.                     };

  236.                 }

  238.                 /** {@inheritDoc} */
  239.                 @Override
  240.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

  245.                     final T omega     = elements.getOmega();
  246.                     final T f         = elements.getF();
  247.                     final T d         = elements.getD();
  248.                     final T t         = elements.getTC();
  249.                     final T t2        = t.multiply(t);

  251.                     final T cosOmega  = omega.cos();
  252.                     final T sinOmega  = omega.sin();
  253.                     final T sin2Omega = omega.multiply(2).sin();
  254.                     final T fMDPO2 = f.subtract(d).add(omega).multiply(2);
  255.                     final T cos2FDOm  = fMDPO2.cos();
  256.                     final T sin2FDOm  = fMDPO2.sin();

  258.                     final T x = xPolynomial.value(t).
  259.                                 add(xy[0].multiply(sinEps0)).
  260.                                 add(t2.multiply(cosOmega.multiply(fXCosOm).add(sinOmega.multiply(fXSinOm)).add(cos2FDOm.multiply(fXSin2FDOm))));
  261.                     final T y = yPolynomial.value(t).
  262.                                 add(xy[1]).
  263.                                 add(t2.multiply(cosOmega.multiply(fYCosOm).add(cos2FDOm.multiply(fYCos2FDOm))));
  264.                     final T sPxy2 = sinOmega.multiply(fSSinOm).
  265.                                     add(sin2Omega.multiply(fSSin2Om)).
  266.                                     add(t.multiply(fST3).add(sinOmega.multiply(fST2SinOm)).add(sin2FDOm.multiply(fST2Sin2FDOm)).multiply(t).add(fST).multiply(t));

  268.                     final T[] a = MathArrays.buildArray(date.getField(), 3);
  269.                     a[0] = x;
  270.                     a[1] = y;
  271.                     a[2] = sPxy2;
  272.                     return a;

  274.                 }

  276.             };

  278.         }

  280.         /** {@inheritDoc} */
  281.         @Override
  282.         public TimeVectorFunction getPrecessionFunction() throws OrekitException {

  284.             // set up the conventional polynomials
  285.             // the following values are from Lieske et al. paper:
  286.             // Expressions for the precession quantities based upon the IAU(1976) system of astronomical constants
  287.             // http://articles.adsabs.harvard.edu/full/1977A%26A....58....1L
  288.             // also available as equation 30 in IERS 2003 conventions
  289.             final PolynomialNutation psiA =
  290.                     new PolynomialNutation(   0.0,
  291.                                            5038.7784   * Constants.ARC_SECONDS_TO_RADIANS,
  292.                                              -1.07259  * Constants.ARC_SECONDS_TO_RADIANS,
  293.                                              -0.001147 * Constants.ARC_SECONDS_TO_RADIANS);
  294.             final PolynomialNutation omegaA =
  295.                     new PolynomialNutation(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
  296.                                             0.0,
  297.                                             0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
  298.                                            -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
  299.             final PolynomialNutation chiA =
  300.                     new PolynomialNutation( 0.0,
  301.                                            10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
  302.                                            -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
  303.                                            -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);

  305.             return new PrecessionFunction(psiA, omegaA, chiA);

  307.         }

  309.         /** {@inheritDoc} */
  310.         @Override
  311.         public TimeVectorFunction getNutationFunction()
  312.             throws OrekitException {

  314.             // set up nutation arguments
  315.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  317.             // set up Poisson series
  318.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  319.             final PoissonSeriesParser baseParser =
  320.                     new PoissonSeriesParser(10).withFirstDelaunay(1);

  322.             final PoissonSeriesParser psiParser =
  323.                     baseParser.
  324.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  325.                     withSinCos(1, 8, deciMilliAS, -1, deciMilliAS);
  326.             final PoissonSeries psiSeries = psiParser.parse(getStream(PSI_EPSILON_SERIES), PSI_EPSILON_SERIES);

  328.             final PoissonSeriesParser epsilonParser =
  329.                     baseParser.
  330.                     withSinCos(0, -1, deciMilliAS, 9, deciMilliAS).
  331.                     withSinCos(1, -1, deciMilliAS, 10, deciMilliAS);
  332.             final PoissonSeries epsilonSeries = epsilonParser.parse(getStream(PSI_EPSILON_SERIES), PSI_EPSILON_SERIES);

  334.             final PoissonSeries.CompiledSeries psiEpsilonSeries =
  335.                     PoissonSeries.compile(psiSeries, epsilonSeries);

  337.             return new TimeVectorFunction() {

  339.                 /** {@inheritDoc} */
  340.                 @Override
  341.                 public double[] value(final AbsoluteDate date) {
  342.                     final BodiesElements elements = arguments.evaluateAll(date);
  343.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  344.                     return new double[] {
  345.                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements)
  346.                     };
  347.                 }

  349.                 /** {@inheritDoc} */
  350.                 @Override
  351.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  352.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  353.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  354.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  355.                     result[0] = psiEpsilon[0];
  356.                     result[1] = psiEpsilon[1];
  357.                     result[2] = IAU1994ResolutionC7.value(elements);
  358.                     return result;
  359.                 }

  361.             };

  363.         }

  365.         /** {@inheritDoc} */
  366.         @Override
  367.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1)
  368.             throws OrekitException {

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

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

  382.             return new TimeScalarFunction() {

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

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

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

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

  400.                 }

  402.                 /** {@inheritDoc} */
  403.                 @Override
  404.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  406.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  407.                     final T dtai = date.durationFrom(gmstReference);
  408.                     final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
  409.                     final T tt   = tut1.divide(Constants.JULIAN_CENTURY);

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

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

  418.                 }

  420.             };

  422.         }

  424.         /** {@inheritDoc} */
  425.         @Override
  426.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1)
  427.             throws OrekitException {

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

  432.             // constants from IERS 1996 page 21
  433.             // the underlying model is IAU 1982 GMST-UT1
  434.             final AbsoluteDate gmstReference =
  435.                 new AbsoluteDate(DateComponents.J2000_EPOCH, TimeComponents.H12, TimeScalesFactory.getTAI());
  436.             final double gmst1 = 8640184.812866;
  437.             final double gmst2 = 0.093104;
  438.             final double gmst3 = -6.2e-6;

  440.             return new TimeScalarFunction() {

  442.                 /** {@inheritDoc} */
  443.                 @Override
  444.                 public double value(final AbsoluteDate date) {

  446.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  447.                     final double dtai = date.durationFrom(gmstReference);
  448.                     final double tut1 = dtai + ut1.offsetFromTAI(date);
  449.                     final double tt   = tut1 / Constants.JULIAN_CENTURY;

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

  454.                 }

  456.                 /** {@inheritDoc} */
  457.                 @Override
  458.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  460.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  461.                     final T dtai = date.durationFrom(gmstReference);
  462.                     final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
  463.                     final T tt   = tut1.divide(Constants.JULIAN_CENTURY);

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

  468.                 }

  470.             };

  472.         }

  474.         /** {@inheritDoc} */
  475.         @Override
  476.         public TimeScalarFunction getGASTFunction(final TimeScale ut1, final EOPHistory eopHistory)
  477.             throws OrekitException {

  479.             // obliquity
  480.             final TimeScalarFunction epsilonA = getMeanObliquityFunction();

  482.             // GMST function
  483.             final TimeScalarFunction gmst = getGMSTFunction(ut1);

  485.             // nutation function
  486.             final TimeVectorFunction nutation = getNutationFunction();

  488.             return new TimeScalarFunction() {

  490.                 /** {@inheritDoc} */
  491.                 @Override
  492.                 public double value(final AbsoluteDate date) {

  494.                     // compute equation of equinoxes
  495.                     final double[] angles = nutation.value(date);
  496.                     double deltaPsi = angles[0];
  497.                     if (eopHistory != null) {
  498.                         deltaPsi += eopHistory.getEquinoxNutationCorrection(date)[0];
  499.                     }
  500.                     final double eqe = deltaPsi  * FastMath.cos(epsilonA.value(date)) + angles[2];

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

  505.                 }

  507.                 /** {@inheritDoc} */
  508.                 @Override
  509.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

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

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

  522.                 }

  524.             };

  526.         }

  528.         /** {@inheritDoc} */
  529.         @Override
  530.         public TimeVectorFunction getEOPTidalCorrection()
  531.             throws OrekitException {

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

  540.             // set up Poisson series
  541.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  542.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(17).
  543.                     withOptionalColumn(1).
  544.                     withGamma(7).
  545.                     withFirstDelaunay(2);
  546.             final PoissonSeries xSeries =
  547.                     xyParser.
  548.                     withSinCos(0, 14, milliAS, 15, milliAS).
  549.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
  550.             final PoissonSeries ySeries =
  551.                     xyParser.
  552.                     withSinCos(0, 16, milliAS, 17, milliAS).
  553.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES),
  554.                                                          TIDAL_CORRECTION_XP_YP_SERIES);

  556.             final double deciMilliS = 1.0e-4;
  557.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(17).
  558.                     withOptionalColumn(1).
  559.                     withGamma(7).
  560.                     withFirstDelaunay(2).
  561.                     withSinCos(0, 16, deciMilliS, 17, deciMilliS);
  562.             final PoissonSeries ut1Series =
  563.                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);

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

  567.         }

  569.         /** {@inheritDoc} */
  570.         public LoveNumbers getLoveNumbers() throws OrekitException {
  571.             return loadLoveNumbers(LOVE_NUMBERS);
  572.         }

  574.         /** {@inheritDoc} */
  575.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1)
  576.             throws OrekitException {

  578.             // set up nutation arguments
  579.             final FundamentalNutationArguments arguments = getNutationArguments(ut1);

  581.             // set up Poisson series
  582.             final PoissonSeriesParser k20Parser =
  583.                     new PoissonSeriesParser(18).
  584.                         withOptionalColumn(1).
  585.                         withDoodson(4, 2).
  586.                         withFirstDelaunay(10);
  587.             final PoissonSeriesParser k21Parser =
  588.                     new PoissonSeriesParser(18).
  589.                         withOptionalColumn(1).
  590.                         withDoodson(4, 3).
  591.                         withFirstDelaunay(10);
  592.             final PoissonSeriesParser k22Parser =
  593.                     new PoissonSeriesParser(16).
  594.                         withOptionalColumn(1).
  595.                         withDoodson(4, 2).
  596.                         withFirstDelaunay(10);

  598.             final double pico = 1.0e-12;
  599.             final PoissonSeries c20Series =
  600.                     k20Parser.
  601.                   withSinCos(0, 18, -pico, 16, pico).
  602.                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
  603.             final PoissonSeries c21Series =
  604.                     k21Parser.
  605.                     withSinCos(0, 17, pico, 18, pico).
  606.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  607.             final PoissonSeries s21Series =
  608.                     k21Parser.
  609.                     withSinCos(0, 18, -pico, 17, pico).
  610.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  611.             final PoissonSeries c22Series =
  612.                     k22Parser.
  613.                     withSinCos(0, -1, pico, 16, pico).
  614.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
  615.             final PoissonSeries s22Series =
  616.                     k22Parser.
  617.                     withSinCos(0, 16, -pico, -1, pico).
  618.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);

  620.             return new TideFrequencyDependenceFunction(arguments,
  621.                                                        c20Series,
  622.                                                        c21Series, s21Series,
  623.                                                        c22Series, s22Series);

  625.         }

  627.         /** {@inheritDoc} */
  628.         @Override
  629.         public double getPermanentTide() throws OrekitException {
  630.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  631.         }

  633.         /** {@inheritDoc} */
  634.         @Override
  635.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  637.             // constants from IERS 1996 page 47
  638.             final double globalFactor = -1.348e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  639.             final double coupling     =  0.00112;

  641.             return new TimeVectorFunction() {

  643.                 /** {@inheritDoc} */
  644.                 @Override
  645.                 public double[] value(final AbsoluteDate date) {
  646.                     final PoleCorrection pole = eopHistory.getPoleCorrection(date);
  647.                     return new double[] {
  648.                         globalFactor * (pole.getXp() + coupling * pole.getYp()),
  649.                         globalFactor * (coupling * pole.getXp() - pole.getYp()),
  650.                     };
  651.                 }

  653.                 /** {@inheritDoc} */
  654.                 @Override
  655.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  656.                     final FieldPoleCorrection<T> pole = eopHistory.getPoleCorrection(date);
  657.                     final T[] a = MathArrays.buildArray(date.getField(), 2);
  658.                     a[0] = pole.getXp().add(pole.getYp().multiply(coupling)).multiply(globalFactor);
  659.                     a[1] = pole.getXp().multiply(coupling).subtract(pole.getYp()).multiply(globalFactor);
  660.                     return a;
  661.                 }

  663.             };
  664.         }

  666.         /** {@inheritDoc} */
  667.         @Override
  668.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory)
  669.             throws OrekitException {

  671.             return new TimeVectorFunction() {

  673.                 /** {@inheritDoc} */
  674.                 @Override
  675.                 public double[] value(final AbsoluteDate date) {
  676.                     // there are no model for ocean pole tide prior to conventions 2010
  677.                     return new double[] {
  678.                         0.0, 0.0
  679.                     };
  680.                 }

  682.                 /** {@inheritDoc} */
  683.                 @Override
  684.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  685.                     // there are no model for ocean pole tide prior to conventions 2010
  686.                     return MathArrays.buildArray(date.getField(), 2);
  687.                 }

  689.             };
  690.         }

  692.         /** {@inheritDoc} */
  693.         @Override
  694.         public double[] getNominalTidalDisplacement() {

  696.             //  // elastic Earth values
  697.             //  return new double[] {
  698.             //      // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  699.             //      0.6026,                                  -0.0006, 0.292, -0.0025, -0.0022,
  700.             //      // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  701.             //      0.0831,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  702.             //      // H₀
  703.             //      -0.31460
  704.             //  };

  706.             // anelastic Earth values
  707.             return new double[] {
  708.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  709.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  710.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  711.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  712.                 // H₀
  713.                 -0.31460
  714.             };

  716.         }

  718.         /** {@inheritDoc} */
  719.         @Override
  720.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal()
  721.             throws OrekitException {
  722.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  723.                                                                   18, 17, -1, 18, -1);
  724.         }

  726.         /** {@inheritDoc} */
  727.         @Override
  728.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal()
  729.             throws OrekitException {
  730.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  731.                                                                 20, 17, 19, 18, 20);
  732.         }

  734.     },

  736.     /** Constant for IERS 2003 conventions. */
  737.     IERS_2003 {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  787.         /** {@inheritDoc} */
  788.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
  789.             throws OrekitException {
  790.             return new FundamentalNutationArguments(this, timeScale,
  791.                                                     getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
  792.         }

  794.         /** {@inheritDoc} */
  795.         @Override
  796.         public TimeScalarFunction getMeanObliquityFunction() throws OrekitException {

  798.             // epsilon 0 value from chapter 5, page 41, other terms from equation 32 page 45
  799.             final PolynomialNutation epsilonA =
  800.                     new PolynomialNutation(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  801.                                              -46.84024  * Constants.ARC_SECONDS_TO_RADIANS,
  802.                                               -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  803.                                                0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  805.             return new TimeScalarFunction() {

  807.                 /** {@inheritDoc} */
  808.                 @Override
  809.                 public double value(final AbsoluteDate date) {
  810.                     return epsilonA.value(evaluateTC(date));
  811.                 }

  813.                 /** {@inheritDoc} */
  814.                 @Override
  815.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  816.                     return epsilonA.value(evaluateTC(date));
  817.                 }

  819.             };

  821.         }

  823.         /** {@inheritDoc} */
  824.         @Override
  825.         public TimeVectorFunction getXYSpXY2Function()
  826.             throws OrekitException {

  828.             // set up nutation arguments
  829.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  831.             // set up Poisson series
  832.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  833.             final PoissonSeriesParser parser =
  834.                     new PoissonSeriesParser(17).
  835.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  836.                         withFirstDelaunay(4).
  837.                         withFirstPlanetary(9).
  838.                         withSinCos(0, 2, microAS, 3, microAS);

  840.             final PoissonSeries xSeries = parser.parse(getStream(X_SERIES), X_SERIES);
  841.             final PoissonSeries ySeries = parser.parse(getStream(Y_SERIES), Y_SERIES);
  842.             final PoissonSeries sSeries = parser.parse(getStream(S_SERIES), S_SERIES);
  843.             final PoissonSeries.CompiledSeries xys = PoissonSeries.compile(xSeries, ySeries, sSeries);

  845.             // create a function evaluating the series
  846.             return new TimeVectorFunction() {

  848.                 /** {@inheritDoc} */
  849.                 @Override
  850.                 public double[] value(final AbsoluteDate date) {
  851.                     return xys.value(arguments.evaluateAll(date));
  852.                 }

  854.                 /** {@inheritDoc} */
  855.                 @Override
  856.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  857.                     return xys.value(arguments.evaluateAll(date));
  858.                 }

  860.             };

  862.         }


  865.         /** {@inheritDoc} */
  866.         @Override
  867.         public TimeVectorFunction getPrecessionFunction() throws OrekitException {

  869.             // set up the conventional polynomials
  870.             // the following values are from equation 32 in IERS 2003 conventions
  871.             final PolynomialNutation psiA =
  872.                     new PolynomialNutation(    0.0,
  873.                                             5038.47875   * Constants.ARC_SECONDS_TO_RADIANS,
  874.                                               -1.07259   * Constants.ARC_SECONDS_TO_RADIANS,
  875.                                               -0.001147  * Constants.ARC_SECONDS_TO_RADIANS);
  876.             final PolynomialNutation omegaA =
  877.                     new PolynomialNutation(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
  878.                                            -0.02524   * Constants.ARC_SECONDS_TO_RADIANS,
  879.                                             0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
  880.                                            -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
  881.             final PolynomialNutation chiA =
  882.                     new PolynomialNutation( 0.0,
  883.                                            10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
  884.                                            -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
  885.                                            -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);

  887.             return new PrecessionFunction(psiA, omegaA, chiA);

  889.         }

  891.         /** {@inheritDoc} */
  892.         @Override
  893.         public TimeVectorFunction getNutationFunction()
  894.             throws OrekitException {

  896.             // set up nutation arguments
  897.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  899.             // set up Poisson series
  900.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  901.             final PoissonSeriesParser luniSolarParser =
  902.                     new PoissonSeriesParser(14).withFirstDelaunay(1);
  903.             final PoissonSeriesParser luniSolarPsiParser =
  904.                     luniSolarParser.
  905.                     withSinCos(0, 7, milliAS, 11, milliAS).
  906.                     withSinCos(1, 8, milliAS, 12, milliAS);
  907.             final PoissonSeries psiLuniSolarSeries =
  908.                     luniSolarPsiParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);
  909.             final PoissonSeriesParser luniSolarEpsilonParser =
  910.                     luniSolarParser.
  911.                     withSinCos(0, 13, milliAS, 9, milliAS).
  912.                     withSinCos(1, 14, milliAS, 10, milliAS);
  913.             final PoissonSeries epsilonLuniSolarSeries =
  914.                     luniSolarEpsilonParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);

  916.             final PoissonSeriesParser planetaryParser =
  917.                     new PoissonSeriesParser(21).
  918.                         withFirstDelaunay(2).
  919.                         withFirstPlanetary(7);
  920.             final PoissonSeriesParser planetaryPsiParser =
  921.                     planetaryParser.withSinCos(0, 17, milliAS, 18, milliAS);
  922.             final PoissonSeries psiPlanetarySeries =
  923.                     planetaryPsiParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);
  924.             final PoissonSeriesParser planetaryEpsilonParser =
  925.                     planetaryParser.withSinCos(0, 19, milliAS, 20, milliAS);
  926.             final PoissonSeries epsilonPlanetarySeries =
  927.                     planetaryEpsilonParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);

  929.             final PoissonSeries.CompiledSeries luniSolarSeries =
  930.                     PoissonSeries.compile(psiLuniSolarSeries, epsilonLuniSolarSeries);
  931.             final PoissonSeries.CompiledSeries planetarySeries =
  932.                     PoissonSeries.compile(psiPlanetarySeries, epsilonPlanetarySeries);

  934.             return new TimeVectorFunction() {

  936.                 /** {@inheritDoc} */
  937.                 @Override
  938.                 public double[] value(final AbsoluteDate date) {
  939.                     final BodiesElements elements = arguments.evaluateAll(date);
  940.                     final double[] luniSolar = luniSolarSeries.value(elements);
  941.                     final double[] planetary = planetarySeries.value(elements);
  942.                     return new double[] {
  943.                         luniSolar[0] + planetary[0], luniSolar[1] + planetary[1],
  944.                         IAU1994ResolutionC7.value(elements)
  945.                     };
  946.                 }

  948.                 /** {@inheritDoc} */
  949.                 @Override
  950.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  951.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  952.                     final T[] luniSolar = luniSolarSeries.value(elements);
  953.                     final T[] planetary = planetarySeries.value(elements);
  954.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  955.                     result[0] = luniSolar[0].add(planetary[0]);
  956.                     result[1] = luniSolar[1].add(planetary[1]);
  957.                     result[2] = IAU1994ResolutionC7.value(elements);
  958.                     return result;
  959.                 }

  961.             };

  963.         }

  965.         /** {@inheritDoc} */
  966.         @Override
  967.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1)
  968.             throws OrekitException {

  970.             // Earth Rotation Angle
  971.             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);

  973.             // Polynomial part of the apparent sidereal time series
  974.             // which is the opposite of Equation of Origins (EO)
  975.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  976.             final PoissonSeriesParser parser =
  977.                     new PoissonSeriesParser(17).
  978.                         withFirstDelaunay(4).
  979.                         withFirstPlanetary(9).
  980.                         withSinCos(0, 2, microAS, 3, microAS).
  981.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  982.             final PolynomialNutation minusEO =
  983.                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();

  985.             // create a function evaluating the series
  986.             return new TimeScalarFunction() {

  988.                 /** {@inheritDoc} */
  989.                 @Override
  990.                 public double value(final AbsoluteDate date) {
  991.                     return era.value(date) + minusEO.value(evaluateTC(date));
  992.                 }

  994.                 /** {@inheritDoc} */
  995.                 @Override
  996.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  997.                     return era.value(date).add(minusEO.value(evaluateTC(date)));
  998.                 }

  1000.             };

  1002.         }

  1004.         /** {@inheritDoc} */
  1005.         @Override
  1006.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1)
  1007.             throws OrekitException {

  1009.             // Earth Rotation Angle
  1010.             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);

  1012.             // Polynomial part of the apparent sidereal time series
  1013.             // which is the opposite of Equation of Origins (EO)
  1014.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1015.             final PoissonSeriesParser parser =
  1016.                     new PoissonSeriesParser(17).
  1017.                         withFirstDelaunay(4).
  1018.                         withFirstPlanetary(9).
  1019.                         withSinCos(0, 2, microAS, 3, microAS).
  1020.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1021.             final PolynomialNutation minusEO =
  1022.                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();

  1024.             // create a function evaluating the series
  1025.             return new TimeScalarFunction() {

  1027.                 /** {@inheritDoc} */
  1028.                 @Override
  1029.                 public double value(final AbsoluteDate date) {
  1030.                     return era.getRate() + minusEO.derivative(evaluateTC(date));
  1031.                 }

  1033.                 /** {@inheritDoc} */
  1034.                 @Override
  1035.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1036.                     return minusEO.derivative(evaluateTC(date)).add(era.getRate());
  1037.                 }

  1039.             };

  1041.         }

  1043.         /** {@inheritDoc} */
  1044.         @Override
  1045.         public TimeScalarFunction getGASTFunction(final TimeScale ut1, final EOPHistory eopHistory)
  1046.             throws OrekitException {

  1048.             // set up nutation arguments
  1049.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  1051.             // mean obliquity function
  1052.             final TimeScalarFunction epsilon = getMeanObliquityFunction();

  1054.             // set up Poisson series
  1055.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  1056.             final PoissonSeriesParser luniSolarPsiParser =
  1057.                     new PoissonSeriesParser(14).
  1058.                         withFirstDelaunay(1).
  1059.                         withSinCos(0, 7, milliAS, 11, milliAS).
  1060.                         withSinCos(1, 8, milliAS, 12, milliAS);
  1061.             final PoissonSeries psiLuniSolarSeries =
  1062.                     luniSolarPsiParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);

  1064.             final PoissonSeriesParser planetaryPsiParser =
  1065.                     new PoissonSeriesParser(21).
  1066.                         withFirstDelaunay(2).
  1067.                         withFirstPlanetary(7).
  1068.                         withSinCos(0, 17, milliAS, 18, milliAS);
  1069.             final PoissonSeries psiPlanetarySeries =
  1070.                     planetaryPsiParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);

  1072.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1073.             final PoissonSeriesParser gstParser =
  1074.                     new PoissonSeriesParser(17).
  1075.                         withFirstDelaunay(4).
  1076.                         withFirstPlanetary(9).
  1077.                         withSinCos(0, 2, microAS, 3, microAS).
  1078.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1079.             final PoissonSeries gstSeries = gstParser.parse(getStream(GST_SERIES), GST_SERIES);
  1080.             final PoissonSeries.CompiledSeries psiGstSeries =
  1081.                     PoissonSeries.compile(psiLuniSolarSeries, psiPlanetarySeries, gstSeries);

  1083.             // ERA function
  1084.             final TimeScalarFunction era = getEarthOrientationAngleFunction(ut1);

  1086.             return new TimeScalarFunction() {

  1088.                 /** {@inheritDoc} */
  1089.                 @Override
  1090.                 public double value(final AbsoluteDate date) {

  1092.                     // evaluate equation of origins
  1093.                     final BodiesElements elements = arguments.evaluateAll(date);
  1094.                     final double[] angles = psiGstSeries.value(elements);
  1095.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  1096.                     final double deltaPsi = angles[0] + angles[1] + ddPsi;
  1097.                     final double epsilonA = epsilon.value(date);

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

  1103.                 }

  1105.                 /** {@inheritDoc} */
  1106.                 @Override
  1107.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  1109.                     // evaluate equation of origins
  1110.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1111.                     final T[] angles = psiGstSeries.value(elements);
  1112.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  1113.                     final T deltaPsi = angles[0].add(angles[1]).add(ddPsi);
  1114.                     final T epsilonA = epsilon.value(date);

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

  1120.                 }

  1122.             };

  1124.         }

  1126.         /** {@inheritDoc} */
  1127.         @Override
  1128.         public TimeVectorFunction getEOPTidalCorrection()
  1129.             throws OrekitException {

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

  1138.             // set up Poisson series
  1139.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1140.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
  1141.                     withOptionalColumn(1).
  1142.                     withGamma(2).
  1143.                     withFirstDelaunay(3);
  1144.             final PoissonSeries xSeries =
  1145.                     xyParser.
  1146.                     withSinCos(0, 10, microAS, 11, microAS).
  1147.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
  1148.             final PoissonSeries ySeries =
  1149.                     xyParser.
  1150.                     withSinCos(0, 12, microAS, 13, microAS).
  1151.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);

  1153.             final double microS = 1.0e-6;
  1154.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  1155.                     withOptionalColumn(1).
  1156.                     withGamma(2).
  1157.                     withFirstDelaunay(3).
  1158.                     withSinCos(0, 10, microS, 11, microS);
  1159.             final PoissonSeries ut1Series =
  1160.                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);

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

  1164.         }

  1166.         /** {@inheritDoc} */
  1167.         public LoveNumbers getLoveNumbers() throws OrekitException {
  1168.             return loadLoveNumbers(LOVE_NUMBERS);
  1169.         }

  1171.         /** {@inheritDoc} */
  1172.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1)
  1173.             throws OrekitException {

  1175.             // set up nutation arguments
  1176.             final FundamentalNutationArguments arguments = getNutationArguments(ut1);

  1178.             // set up Poisson series
  1179.             final PoissonSeriesParser k20Parser =
  1180.                     new PoissonSeriesParser(18).
  1181.                         withOptionalColumn(1).
  1182.                         withDoodson(4, 2).
  1183.                         withFirstDelaunay(10);
  1184.             final PoissonSeriesParser k21Parser =
  1185.                     new PoissonSeriesParser(18).
  1186.                         withOptionalColumn(1).
  1187.                         withDoodson(4, 3).
  1188.                         withFirstDelaunay(10);
  1189.             final PoissonSeriesParser k22Parser =
  1190.                     new PoissonSeriesParser(16).
  1191.                         withOptionalColumn(1).
  1192.                         withDoodson(4, 2).
  1193.                         withFirstDelaunay(10);

  1195.             final double pico = 1.0e-12;
  1196.             final PoissonSeries c20Series =
  1197.                     k20Parser.
  1198.                   withSinCos(0, 18, -pico, 16, pico).
  1199.                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
  1200.             final PoissonSeries c21Series =
  1201.                     k21Parser.
  1202.                     withSinCos(0, 17, pico, 18, pico).
  1203.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  1204.             final PoissonSeries s21Series =
  1205.                     k21Parser.
  1206.                     withSinCos(0, 18, -pico, 17, pico).
  1207.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  1208.             final PoissonSeries c22Series =
  1209.                     k22Parser.
  1210.                     withSinCos(0, -1, pico, 16, pico).
  1211.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
  1212.             final PoissonSeries s22Series =
  1213.                     k22Parser.
  1214.                     withSinCos(0, 16, -pico, -1, pico).
  1215.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);

  1217.             return new TideFrequencyDependenceFunction(arguments,
  1218.                                                        c20Series,
  1219.                                                        c21Series, s21Series,
  1220.                                                        c22Series, s22Series);

  1222.         }

  1224.         /** {@inheritDoc} */
  1225.         @Override
  1226.         public double getPermanentTide() throws OrekitException {
  1227.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1228.         }

  1230.         /** {@inheritDoc} */
  1231.         @Override
  1232.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory)
  1233.             throws OrekitException {

  1235.             // annual pole from ftp://tai.bipm.org/iers/conv2003/chapter7/annual.pole
  1236.             final TimeScale utc = TimeScalesFactory.getUTC();
  1237.             final SimpleTimeStampedTableParser.RowConverter<MeanPole> converter =
  1238.                 new SimpleTimeStampedTableParser.RowConverter<MeanPole>() {
  1239.                     /** {@inheritDoc} */
  1240.                     @Override
  1241.                     public MeanPole convert(final double[] rawFields) throws OrekitException {
  1242.                         return new MeanPole(new AbsoluteDate((int) rawFields[0], 1, 1, utc),
  1243.                                             rawFields[1] * Constants.ARC_SECONDS_TO_RADIANS,
  1244.                                             rawFields[2] * Constants.ARC_SECONDS_TO_RADIANS);
  1245.                     }
  1246.                 };
  1247.             final SimpleTimeStampedTableParser<MeanPole> parser =
  1248.                     new SimpleTimeStampedTableParser<MeanPole>(3, converter);
  1249.             final List<MeanPole> annualPoleList = parser.parse(getStream(ANNUAL_POLE), ANNUAL_POLE);
  1250.             final AbsoluteDate firstAnnualPoleDate = annualPoleList.get(0).getDate();
  1251.             final AbsoluteDate lastAnnualPoleDate  = annualPoleList.get(annualPoleList.size() - 1).getDate();
  1252.             final ImmutableTimeStampedCache<MeanPole> annualCache =
  1253.                     new ImmutableTimeStampedCache<MeanPole>(2, annualPoleList);

  1255.             // polynomial extension from IERS 2003, section 7.1.4, equations 23a and 23b
  1256.             final double xp0    = 0.054   * Constants.ARC_SECONDS_TO_RADIANS;
  1257.             final double xp0Dot = 0.00083 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
  1258.             final double yp0    = 0.357   * Constants.ARC_SECONDS_TO_RADIANS;
  1259.             final double yp0Dot = 0.00395 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;

  1261.             // constants from IERS 2003, section 6.2
  1262.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1263.             final double ratio        =  0.00115;

  1265.             return new TimeVectorFunction() {

  1267.                 /** {@inheritDoc} */
  1268.                 @Override
  1269.                 public double[] value(final AbsoluteDate date) {

  1271.                     // we can't compute anything before the range covered by the annual pole file
  1272.                     if (date.compareTo(firstAnnualPoleDate) <= 0) {
  1273.                         return new double[] {
  1274.                             0.0, 0.0
  1275.                         };
  1276.                     }

  1278.                     // evaluate mean pole
  1279.                     double meanPoleX = 0;
  1280.                     double meanPoleY = 0;
  1281.                     if (date.compareTo(lastAnnualPoleDate) <= 0) {
  1282.                         // we are within the range covered by the annual pole file,
  1283.                         // we interpolate within it
  1284.                         try {
  1285.                             final HermiteInterpolator interpolator = new HermiteInterpolator();
  1286.                             annualCache.getNeighbors(date).forEach(neighbor ->
  1287.                                 interpolator.addSamplePoint(neighbor.getDate().durationFrom(date),
  1288.                                                             new double[] {
  1289.                                                                 neighbor.getX(), neighbor.getY()
  1290.                                                             }));
  1291.                             final double[] interpolated = interpolator.value(0);
  1292.                             meanPoleX = interpolated[0];
  1293.                             meanPoleY = interpolated[1];
  1294.                         } catch (TimeStampedCacheException tsce) {
  1295.                             // this should never happen
  1296.                             throw new OrekitInternalError(tsce);
  1297.                         }
  1298.                     } else {

  1300.                         // we are after the range covered by the annual pole file,
  1301.                         // we use the polynomial extension
  1302.                         final double t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
  1303.                         meanPoleX = xp0 + t * xp0Dot;
  1304.                         meanPoleY = yp0 + t * yp0Dot;

  1306.                     }

  1308.                     // evaluate wobble variables
  1309.                     final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1310.                     final double m1 = correction.getXp() - meanPoleX;
  1311.                     final double m2 = meanPoleY - correction.getYp();

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

  1344.                 }

  1346.                 /** {@inheritDoc} */
  1347.                 @Override
  1348.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

  1357.                     // evaluate mean pole
  1358.                     T meanPoleX = date.getField().getZero();
  1359.                     T meanPoleY = date.getField().getZero();
  1360.                     if (aDate.compareTo(lastAnnualPoleDate) <= 0) {
  1361.                         // we are within the range covered by the annual pole file,
  1362.                         // we interpolate within it
  1363.                         try {
  1364.                             final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
  1365.                             final T[] y = MathArrays.buildArray(date.getField(), 2);
  1366.                             final T zero = date.getField().getZero();
  1367.                             final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
  1368.                                                                                                        // for example removing derivatives
  1369.                                                                                                        // if T was DerivativeStructure
  1370.                             annualCache.getNeighbors(aDate).forEach(neighbor -> {
  1371.                                 y[0] = zero.add(neighbor.getX());
  1372.                                 y[1] = zero.add(neighbor.getY());
  1373.                                 interpolator.addSamplePoint(central.durationFrom(neighbor.getDate()).negate(), y);
  1374.                             });
  1375.                             final T[] interpolated = interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
  1376.                             meanPoleX = interpolated[0];
  1377.                             meanPoleY = interpolated[1];
  1378.                         } catch (TimeStampedCacheException tsce) {
  1379.                             // this should never happen
  1380.                             throw new OrekitInternalError(tsce);
  1381.                         }
  1382.                     } else {

  1384.                         // we are after the range covered by the annual pole file,
  1385.                         // we use the polynomial extension
  1386.                         final T t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
  1387.                         meanPoleX = t.multiply(xp0Dot).add(xp0);
  1388.                         meanPoleY = t.multiply(yp0Dot).add(yp0);

  1390.                     }

  1392.                     // evaluate wobble variables
  1393.                     final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1394.                     final T m1 = correction.getXp().subtract(meanPoleX);
  1395.                     final T m2 = meanPoleY.subtract(correction.getYp());

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

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

  1428.                     return a;

  1430.                 }

  1432.             };

  1434.         }

  1436.         /** {@inheritDoc} */
  1437.         @Override
  1438.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory)
  1439.             throws OrekitException {

  1441.             return new TimeVectorFunction() {

  1443.                 /** {@inheritDoc} */
  1444.                 @Override
  1445.                 public double[] value(final AbsoluteDate date) {
  1446.                     // there are no model for ocean pole tide prior to conventions 2010
  1447.                     return new double[] {
  1448.                         0.0, 0.0
  1449.                     };
  1450.                 }

  1452.                 /** {@inheritDoc} */
  1453.                 @Override
  1454.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1455.                     // there are no model for ocean pole tide prior to conventions 2010
  1456.                     return MathArrays.buildArray(date.getField(), 2);
  1457.                 }

  1459.             };
  1460.         }

  1462.         /** {@inheritDoc} */
  1463.         @Override
  1464.         public double[] getNominalTidalDisplacement() {
  1465.             return new double[] {
  1466.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  1467.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  1468.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  1469.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  1470.                 // H₀
  1471.                 -0.31460
  1472.             };
  1473.         }

  1475.         /** {@inheritDoc} */
  1476.         @Override
  1477.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal()
  1478.             throws OrekitException {
  1479.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  1480.                                                                   18, 15, 16, 17, 18);
  1481.         }

  1483.         /** {@inheritDoc} */
  1484.         @Override
  1485.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal()
  1486.             throws OrekitException {
  1487.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  1488.                                                                 18, 15, 16, 17, 18);
  1489.         }

  1491.     },

  1493.     /** Constant for IERS 2010 conventions. */
  1494.     IERS_2010 {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  1541.         /** {@inheritDoc} */
  1542.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
  1543.             throws OrekitException {
  1544.             return new FundamentalNutationArguments(this, timeScale,
  1545.                                                     getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
  1546.         }

  1548.         /** {@inheritDoc} */
  1549.         @Override
  1550.         public TimeScalarFunction getMeanObliquityFunction() throws OrekitException {

  1552.             // epsilon 0 value from chapter 5, page 56, other terms from equation 5.40 page 65
  1553.             final PolynomialNutation epsilonA =
  1554.                     new PolynomialNutation(84381.406        * Constants.ARC_SECONDS_TO_RADIANS,
  1555.                                              -46.836769     * Constants.ARC_SECONDS_TO_RADIANS,
  1556.                                               -0.0001831    * Constants.ARC_SECONDS_TO_RADIANS,
  1557.                                                0.00200340   * Constants.ARC_SECONDS_TO_RADIANS,
  1558.                                               -0.000000576  * Constants.ARC_SECONDS_TO_RADIANS,
  1559.                                               -0.0000000434 * Constants.ARC_SECONDS_TO_RADIANS);

  1561.             return new TimeScalarFunction() {

  1563.                 /** {@inheritDoc} */
  1564.                 @Override
  1565.                 public double value(final AbsoluteDate date) {
  1566.                     return epsilonA.value(evaluateTC(date));
  1567.                 }

  1569.                 /** {@inheritDoc} */
  1570.                 @Override
  1571.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1572.                     return epsilonA.value(evaluateTC(date));
  1573.                 }

  1575.             };

  1577.         }

  1579.         /** {@inheritDoc} */
  1580.         @Override
  1581.         public TimeVectorFunction getXYSpXY2Function() throws OrekitException {

  1583.             // set up nutation arguments
  1584.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  1586.             // set up Poisson series
  1587.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1588.             final PoissonSeriesParser parser =
  1589.                     new PoissonSeriesParser(17).
  1590.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  1591.                         withFirstDelaunay(4).
  1592.                         withFirstPlanetary(9).
  1593.                         withSinCos(0, 2, microAS, 3, microAS);
  1594.             final PoissonSeries xSeries = parser.parse(getStream(X_SERIES), X_SERIES);
  1595.             final PoissonSeries ySeries = parser.parse(getStream(Y_SERIES), Y_SERIES);
  1596.             final PoissonSeries sSeries = parser.parse(getStream(S_SERIES), S_SERIES);
  1597.             final PoissonSeries.CompiledSeries xys = PoissonSeries.compile(xSeries, ySeries, sSeries);

  1599.             // create a function evaluating the series
  1600.             return new TimeVectorFunction() {

  1602.                 /** {@inheritDoc} */
  1603.                 @Override
  1604.                 public double[] value(final AbsoluteDate date) {
  1605.                     return xys.value(arguments.evaluateAll(date));
  1606.                 }

  1608.                 /** {@inheritDoc} */
  1609.                 @Override
  1610.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1611.                     return xys.value(arguments.evaluateAll(date));
  1612.                 }

  1614.             };

  1616.         }

  1618.         /** {@inheritDoc} */
  1619.         public LoveNumbers getLoveNumbers() throws OrekitException {
  1620.             return loadLoveNumbers(LOVE_NUMBERS);
  1621.         }

  1623.         /** {@inheritDoc} */
  1624.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1)
  1625.             throws OrekitException {

  1627.             // set up nutation arguments
  1628.             final FundamentalNutationArguments arguments = getNutationArguments(ut1);

  1630.             // set up Poisson series
  1631.             final PoissonSeriesParser k20Parser =
  1632.                     new PoissonSeriesParser(18).
  1633.                         withOptionalColumn(1).
  1634.                         withDoodson(4, 2).
  1635.                         withFirstDelaunay(10);
  1636.             final PoissonSeriesParser k21Parser =
  1637.                     new PoissonSeriesParser(18).
  1638.                         withOptionalColumn(1).
  1639.                         withDoodson(4, 3).
  1640.                         withFirstDelaunay(10);
  1641.             final PoissonSeriesParser k22Parser =
  1642.                     new PoissonSeriesParser(16).
  1643.                         withOptionalColumn(1).
  1644.                         withDoodson(4, 2).
  1645.                         withFirstDelaunay(10);

  1647.             final double pico = 1.0e-12;
  1648.             final PoissonSeries c20Series =
  1649.                     k20Parser.
  1650.                     withSinCos(0, 18, -pico, 16, pico).
  1651.                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
  1652.             final PoissonSeries c21Series =
  1653.                     k21Parser.
  1654.                     withSinCos(0, 17, pico, 18, pico).
  1655.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  1656.             final PoissonSeries s21Series =
  1657.                     k21Parser.
  1658.                     withSinCos(0, 18, -pico, 17, pico).
  1659.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  1660.             final PoissonSeries c22Series =
  1661.                     k22Parser.
  1662.                     withSinCos(0, -1, pico, 16, pico).
  1663.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
  1664.             final PoissonSeries s22Series =
  1665.                     k22Parser.
  1666.                     withSinCos(0, 16, -pico, -1, pico).
  1667.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);

  1669.             return new TideFrequencyDependenceFunction(arguments,
  1670.                                                        c20Series,
  1671.                                                        c21Series, s21Series,
  1672.                                                        c22Series, s22Series);

  1674.         }

  1676.         /** {@inheritDoc} */
  1677.         @Override
  1678.         public double getPermanentTide() throws OrekitException {
  1679.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1680.         }

  1682.         /** Compute pole wobble variables m₁ and m₂.
  1683.          * @param date current date
  1684.          * @param eopHistory EOP history
  1685.          * @return array containing m₁ and m₂
  1686.          */
  1687.         private double[] computePoleWobble(final AbsoluteDate date, final EOPHistory eopHistory) {

  1689.             // polynomial model from IERS 2010, table 7.7
  1690.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1691.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1692.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1693.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1694.             final AbsoluteDate changeDate = new AbsoluteDate(2010, 1, 1, TimeScalesFactory.getTT());

  1696.             // evaluate mean pole
  1697.             final double[] xPolynomial;
  1698.             final double[] yPolynomial;
  1699.             if (date.compareTo(changeDate) <= 0) {
  1700.                 xPolynomial = new double[] {
  1701.                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
  1702.                 };
  1703.                 yPolynomial = new double[] {
  1704.                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
  1705.                 };
  1706.             } else {
  1707.                 xPolynomial = new double[] {
  1708.                     23.513 * f0, 7.6141 * f1
  1709.                 };
  1710.                 yPolynomial = new double[] {
  1711.                     358.891 * f0,  -0.6287 * f1
  1712.                 };
  1713.             }
  1714.             double meanPoleX = 0;
  1715.             double meanPoleY = 0;
  1716.             final double t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
  1717.             for (int i = xPolynomial.length - 1; i >= 0; --i) {
  1718.                 meanPoleX = meanPoleX * t + xPolynomial[i];
  1719.             }
  1720.             for (int i = yPolynomial.length - 1; i >= 0; --i) {
  1721.                 meanPoleY = meanPoleY * t + yPolynomial[i];
  1722.             }

  1724.             // evaluate wobble variables
  1725.             final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1726.             final double m1 = correction.getXp() - meanPoleX;
  1727.             final double m2 = meanPoleY - correction.getYp();

  1729.             return new double[] {
  1730.                 m1, m2
  1731.             };

  1733.         }

  1735.         /** Compute pole wobble variables m₁ and m₂.
  1736.          * @param date current date
  1737.          * @param <T> type of the field elements
  1738.          * @param eopHistory EOP history
  1739.          * @return array containing m₁ and m₂
  1740.          */
  1741.         private <T extends RealFieldElement<T>> T[] computePoleWobble(final FieldAbsoluteDate<T> date, final EOPHistory eopHistory) {

  1743.             // polynomial model from IERS 2010, table 7.7
  1744.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1745.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1746.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1747.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1748.             final AbsoluteDate changeDate = new AbsoluteDate(2010, 1, 1, TimeScalesFactory.getTT());

  1750.             // evaluate mean pole
  1751.             final double[] xPolynomial;
  1752.             final double[] yPolynomial;
  1753.             if (date.toAbsoluteDate().compareTo(changeDate) <= 0) {
  1754.                 xPolynomial = new double[] {
  1755.                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
  1756.                 };
  1757.                 yPolynomial = new double[] {
  1758.                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
  1759.                 };
  1760.             } else {
  1761.                 xPolynomial = new double[] {
  1762.                     23.513 * f0, 7.6141 * f1
  1763.                 };
  1764.                 yPolynomial = new double[] {
  1765.                     358.891 * f0,  -0.6287 * f1
  1766.                 };
  1767.             }
  1768.             T meanPoleX = date.getField().getZero();
  1769.             T meanPoleY = date.getField().getZero();
  1770.             final T t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
  1771.             for (int i = xPolynomial.length - 1; i >= 0; --i) {
  1772.                 meanPoleX = meanPoleX.multiply(t).add(xPolynomial[i]);
  1773.             }
  1774.             for (int i = yPolynomial.length - 1; i >= 0; --i) {
  1775.                 meanPoleY = meanPoleY.multiply(t).add(yPolynomial[i]);
  1776.             }

  1778.             // evaluate wobble variables
  1779.             final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1780.             final T[] m = MathArrays.buildArray(date.getField(), 2);
  1781.             m[0] = correction.getXp().subtract(meanPoleX);
  1782.             m[1] = meanPoleY.subtract(correction.getYp());

  1784.             return m;

  1786.         }

  1788.         /** {@inheritDoc} */
  1789.         @Override
  1790.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory)
  1791.             throws OrekitException {

  1793.             // constants from IERS 2010, section 6.4
  1794.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1795.             final double ratio        =  0.00115;

  1797.             return new TimeVectorFunction() {

  1799.                 /** {@inheritDoc} */
  1800.                 @Override
  1801.                 public double[] value(final AbsoluteDate date) {

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

  1806.                     return new double[] {
  1807.                         // the following correspond to the equations published in IERS 2010 conventions,
  1808.                         // section 6.4 page 94. The equations read:
  1809.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1810.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1811.                         // These equations seem to fix what was probably a sign error in IERS 2003
  1812.                         // conventions section 6.2 page 65. In this older publication, the equations read:
  1813.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1814.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1815.                         globalFactor * (wobbleM[0] + ratio * wobbleM[1]),
  1816.                         globalFactor * (wobbleM[1] - ratio * wobbleM[0])
  1817.                     };

  1819.                 }

  1821.                 /** {@inheritDoc} */
  1822.                 @Override
  1823.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

  1830.                     // the following correspond to the equations published in IERS 2010 conventions,
  1831.                     // section 6.4 page 94. The equations read:
  1832.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1833.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1834.                     // These equations seem to fix what was probably a sign error in IERS 2003
  1835.                     // conventions section 6.2 page 65. In this older publication, the equations read:
  1836.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1837.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1838.                     a[0] = wobbleM[0].add(wobbleM[1].multiply( ratio)).multiply(globalFactor);
  1839.                     a[1] = wobbleM[1].add(wobbleM[0].multiply(-ratio)).multiply(globalFactor);

  1841.                     return a;

  1843.                 }

  1845.             };

  1847.         }

  1849.         /** {@inheritDoc} */
  1850.         @Override
  1851.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory)
  1852.             throws OrekitException {

  1854.             return new TimeVectorFunction() {

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

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

  1863.                     return new double[] {
  1864.                         // the following correspond to the equations published in IERS 2010 conventions,
  1865.                         // section 6.4 page 94 equation 6.24:
  1866.                         // ∆C₂₁ = −2.1778 × 10⁻¹⁰ (m₁ − 0.01724m₂)
  1867.                         // ∆S₂₁ = −1.7232 × 10⁻¹⁰ (m₂ − 0.03365m₁)
  1868.                         -2.1778e-10 * (wobbleM[0] - 0.01724 * wobbleM[1]) / Constants.ARC_SECONDS_TO_RADIANS,
  1869.                         -1.7232e-10 * (wobbleM[1] - 0.03365 * wobbleM[0]) / Constants.ARC_SECONDS_TO_RADIANS
  1870.                     };

  1872.                 }

  1874.                 /** {@inheritDoc} */
  1875.                 @Override
  1876.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

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

  1890.                     return v;

  1892.                 }

  1894.             };

  1896.         }

  1898.         /** {@inheritDoc} */
  1899.         @Override
  1900.         public TimeVectorFunction getPrecessionFunction() throws OrekitException {

  1902.             // set up the conventional polynomials
  1903.             // the following values are from equation 5.40 in IERS 2010 conventions
  1904.             final PolynomialNutation psiA =
  1905.                     new PolynomialNutation(   0.0,
  1906.                                            5038.481507     * Constants.ARC_SECONDS_TO_RADIANS,
  1907.                                              -1.0790069    * Constants.ARC_SECONDS_TO_RADIANS,
  1908.                                              -0.00114045   * Constants.ARC_SECONDS_TO_RADIANS,
  1909.                                               0.000132851  * Constants.ARC_SECONDS_TO_RADIANS,
  1910.                                              -0.0000000951 * Constants.ARC_SECONDS_TO_RADIANS);
  1911.             final PolynomialNutation omegaA =
  1912.                     new PolynomialNutation(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
  1913.                                            -0.025754     * Constants.ARC_SECONDS_TO_RADIANS,
  1914.                                             0.0512623    * Constants.ARC_SECONDS_TO_RADIANS,
  1915.                                            -0.00772503   * Constants.ARC_SECONDS_TO_RADIANS,
  1916.                                            -0.000000467  * Constants.ARC_SECONDS_TO_RADIANS,
  1917.                                             0.0000003337 * Constants.ARC_SECONDS_TO_RADIANS);
  1918.             final PolynomialNutation chiA =
  1919.                     new PolynomialNutation( 0.0,
  1920.                                            10.556403     * Constants.ARC_SECONDS_TO_RADIANS,
  1921.                                            -2.3814292    * Constants.ARC_SECONDS_TO_RADIANS,
  1922.                                            -0.00121197   * Constants.ARC_SECONDS_TO_RADIANS,
  1923.                                             0.000170663  * Constants.ARC_SECONDS_TO_RADIANS,
  1924.                                            -0.0000000560 * Constants.ARC_SECONDS_TO_RADIANS);

  1926.             return new PrecessionFunction(psiA, omegaA, chiA);

  1928.         }

  1930.          /** {@inheritDoc} */
  1931.         @Override
  1932.         public TimeVectorFunction getNutationFunction()
  1933.             throws OrekitException {

  1935.             // set up nutation arguments
  1936.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  1938.             // set up Poisson series
  1939.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1940.             final PoissonSeriesParser parser =
  1941.                     new PoissonSeriesParser(17).
  1942.                         withFirstDelaunay(4).
  1943.                         withFirstPlanetary(9).
  1944.                         withSinCos(0, 2, microAS, 3, microAS);
  1945.             final PoissonSeries psiSeries     = parser.parse(getStream(PSI_SERIES), PSI_SERIES);
  1946.             final PoissonSeries epsilonSeries = parser.parse(getStream(EPSILON_SERIES), EPSILON_SERIES);
  1947.             final PoissonSeries.CompiledSeries psiEpsilonSeries =
  1948.                     PoissonSeries.compile(psiSeries, epsilonSeries);

  1950.             return new TimeVectorFunction() {

  1952.                 /** {@inheritDoc} */
  1953.                 @Override
  1954.                 public double[] value(final AbsoluteDate date) {
  1955.                     final BodiesElements elements = arguments.evaluateAll(date);
  1956.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  1957.                     return new double[] {
  1958.                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements)
  1959.                     };
  1960.                 }

  1962.                 /** {@inheritDoc} */
  1963.                 @Override
  1964.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1965.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1966.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  1967.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  1968.                     result[0] = psiEpsilon[0];
  1969.                     result[1] = psiEpsilon[1];
  1970.                     result[2] = IAU1994ResolutionC7.value(elements);
  1971.                     return result;
  1972.                 }

  1974.             };

  1976.         }

  1978.         /** {@inheritDoc} */
  1979.         @Override
  1980.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1) throws OrekitException {

  1982.             // Earth Rotation Angle
  1983.             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);

  1985.             // Polynomial part of the apparent sidereal time series
  1986.             // which is the opposite of Equation of Origins (EO)
  1987.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1988.             final PoissonSeriesParser parser =
  1989.                     new PoissonSeriesParser(17).
  1990.                         withFirstDelaunay(4).
  1991.                         withFirstPlanetary(9).
  1992.                         withSinCos(0, 2, microAS, 3, microAS).
  1993.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1994.             final PolynomialNutation minusEO =
  1995.                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();

  1997.             // create a function evaluating the series
  1998.             return new TimeScalarFunction() {

  2000.                 /** {@inheritDoc} */
  2001.                 @Override
  2002.                 public double value(final AbsoluteDate date) {
  2003.                     return era.value(date) + minusEO.value(evaluateTC(date));
  2004.                 }

  2006.                 /** {@inheritDoc} */
  2007.                 @Override
  2008.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  2009.                     return era.value(date).add(minusEO.value(evaluateTC(date)));
  2010.                 }

  2012.             };

  2014.         }

  2016.         /** {@inheritDoc} */
  2017.         @Override
  2018.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1) throws OrekitException {

  2020.             // Earth Rotation Angle
  2021.             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);

  2023.             // Polynomial part of the apparent sidereal time series
  2024.             // which is the opposite of Equation of Origins (EO)
  2025.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2026.             final PoissonSeriesParser parser =
  2027.                     new PoissonSeriesParser(17).
  2028.                         withFirstDelaunay(4).
  2029.                         withFirstPlanetary(9).
  2030.                         withSinCos(0, 2, microAS, 3, microAS).
  2031.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  2032.             final PolynomialNutation minusEO =
  2033.                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();

  2035.             // create a function evaluating the series
  2036.             return new TimeScalarFunction() {

  2038.                 /** {@inheritDoc} */
  2039.                 @Override
  2040.                 public double value(final AbsoluteDate date) {
  2041.                     return era.getRate() + minusEO.derivative(evaluateTC(date));
  2042.                 }

  2044.                 /** {@inheritDoc} */
  2045.                 @Override
  2046.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  2047.                     return minusEO.derivative(evaluateTC(date)).add(era.getRate());
  2048.                 }

  2050.             };

  2052.         }

  2054.         /** {@inheritDoc} */
  2055.         @Override
  2056.         public TimeScalarFunction getGASTFunction(final TimeScale ut1, final EOPHistory eopHistory)
  2057.             throws OrekitException {

  2059.             // set up nutation arguments
  2060.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  2062.             // mean obliquity function
  2063.             final TimeScalarFunction epsilon = getMeanObliquityFunction();

  2065.             // set up Poisson series
  2066.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2067.             final PoissonSeriesParser baseParser =
  2068.                     new PoissonSeriesParser(17).
  2069.                         withFirstDelaunay(4).
  2070.                         withFirstPlanetary(9).
  2071.                         withSinCos(0, 2, microAS, 3, microAS);
  2072.             final PoissonSeriesParser gstParser  = baseParser.withPolynomialPart('t', Unit.ARC_SECONDS);
  2073.             final PoissonSeries psiSeries        = baseParser.parse(getStream(PSI_SERIES), PSI_SERIES);
  2074.             final PoissonSeries gstSeries        = gstParser.parse(getStream(GST_SERIES), GST_SERIES);
  2075.             final PoissonSeries.CompiledSeries psiGstSeries =
  2076.                     PoissonSeries.compile(psiSeries, gstSeries);

  2078.             // ERA function
  2079.             final TimeScalarFunction era = getEarthOrientationAngleFunction(ut1);

  2081.             return new TimeScalarFunction() {

  2083.                 /** {@inheritDoc} */
  2084.                 @Override
  2085.                 public double value(final AbsoluteDate date) {

  2087.                     // evaluate equation of origins
  2088.                     final BodiesElements elements = arguments.evaluateAll(date);
  2089.                     final double[] angles = psiGstSeries.value(elements);
  2090.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  2091.                     final double deltaPsi = angles[0] + ddPsi;
  2092.                     final double epsilonA = epsilon.value(date);

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

  2098.                 }

  2100.                 /** {@inheritDoc} */
  2101.                 @Override
  2102.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  2104.                     // evaluate equation of origins
  2105.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  2106.                     final T[] angles = psiGstSeries.value(elements);
  2107.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  2108.                     final T deltaPsi = angles[0].add(ddPsi);
  2109.                     final T epsilonA = epsilon.value(date);

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

  2115.                 }

  2117.             };

  2119.         }

  2121.         /** {@inheritDoc} */
  2122.         @Override
  2123.         public TimeVectorFunction getEOPTidalCorrection()
  2124.             throws OrekitException {

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

  2133.             // set up Poisson series
  2134.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2135.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
  2136.                     withOptionalColumn(1).
  2137.                     withGamma(2).
  2138.                     withFirstDelaunay(3);
  2139.             final PoissonSeries xSeries =
  2140.                     xyParser.
  2141.                     withSinCos(0, 10, microAS, 11, microAS).
  2142.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
  2143.             final PoissonSeries ySeries =
  2144.                     xyParser.
  2145.                     withSinCos(0, 12, microAS, 13, microAS).
  2146.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);

  2148.             final double microS = 1.0e-6;
  2149.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  2150.                     withOptionalColumn(1).
  2151.                     withGamma(2).
  2152.                     withFirstDelaunay(3).
  2153.                     withSinCos(0, 10, microS, 11, microS);
  2154.             final PoissonSeries ut1Series =
  2155.                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);

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

  2159.         }

  2161.         /** {@inheritDoc} */
  2162.         @Override
  2163.         public double[] getNominalTidalDisplacement() {
  2164.             return new double[] {
  2165.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  2166.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  2167.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  2168.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  2169.                 // H₀
  2170.                 -0.31460
  2171.             };
  2172.         }

  2174.         /** {@inheritDoc} */
  2175.         @Override
  2176.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal()
  2177.             throws OrekitException {
  2178.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  2179.                                                                   18, 15, 16, 17, 18);
  2180.         }

  2182.         /** {@inheritDoc} */
  2183.         @Override
  2184.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal()
  2185.             throws OrekitException {
  2186.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  2187.                                                                 18, 15, 16, 17, 18);
  2188.         }

  2190.     };

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

  2195.     /** Get the reference epoch for fundamental nutation arguments.
  2196.      * @return reference epoch for fundamental nutation arguments
  2197.      * @since 6.1
  2198.      */
  2199.     public AbsoluteDate getNutationReferenceEpoch() {
  2200.         // IERS 1996, IERS 2003 and IERS 2010 use the same J2000.0 reference date
  2201.         return AbsoluteDate.J2000_EPOCH;
  2202.     }

  2204.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2205.      * @param date current date
  2206.      * @return date offset in Julian centuries
  2207.      * @since 6.1
  2208.      */
  2209.     public double evaluateTC(final AbsoluteDate date) {
  2210.         return date.durationFrom(getNutationReferenceEpoch()) / Constants.JULIAN_CENTURY;
  2211.     }

  2213.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2214.      * @param date current date
  2215.      * @param <T> type of the field elements
  2216.      * @return date offset in Julian centuries
  2217.      * @since 9.0
  2218.      */
  2219.     public <T extends RealFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date) {
  2220.         return date.durationFrom(getNutationReferenceEpoch()).divide(Constants.JULIAN_CENTURY);
  2221.     }

  2223.     /** Get the fundamental nutation arguments.
  2224.      * @param timeScale time scale for computing Greenwich Mean Sidereal Time
  2225.      * (typically {@link TimeScalesFactory#getUT1(IERSConventions, boolean) UT1})
  2226.      * @return fundamental nutation arguments
  2227.      * @exception OrekitException if fundamental nutation arguments cannot be loaded
  2228.      * @since 6.1
  2229.      */
  2230.     public abstract FundamentalNutationArguments getNutationArguments(TimeScale timeScale)
  2231.         throws OrekitException;

  2233.     /** Get the function computing mean obliquity of the ecliptic.
  2234.      * @return function computing mean obliquity of the ecliptic
  2235.      * @exception OrekitException if table cannot be loaded
  2236.      * @since 6.1
  2237.      */
  2238.     public abstract TimeScalarFunction getMeanObliquityFunction() throws OrekitException;

  2240.     /** Get the function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components.
  2241.      * <p>
  2242.      * The returned function computes the two X, Y components of CIP and the S+XY/2 component of the non-rotating CIO.
  2243.      * </p>
  2244.      * @return function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components
  2245.      * @exception OrekitException if table cannot be loaded
  2246.      * @since 6.1
  2247.      */
  2248.     public abstract TimeVectorFunction getXYSpXY2Function()
  2249.         throws OrekitException;

  2251.     /** Get the function computing the raw Earth Orientation Angle.
  2252.      * <p>
  2253.      * The raw angle does not contain any correction. If for example dTU1 correction
  2254.      * due to tidal effect is desired, it must be added afterward by the caller.
  2255.      * The returned value contain the angle as the value and the angular rate as
  2256.      * the first derivative.
  2257.      * </p>
  2258.      * @param ut1 UT1 time scale
  2259.      * @return function computing the rawEarth Orientation Angle, in the non-rotating origin paradigm
  2260.      * @since 6.1
  2261.      */
  2262.     public TimeScalarFunction getEarthOrientationAngleFunction(final TimeScale ut1) {
  2263.         return new StellarAngleCapitaine(ut1);
  2264.     }


  2267.     /** Get the function computing the precession angles.
  2268.      * <p>
  2269.      * The function returned computes the three precession angles
  2270.      * ψ<sub>A</sub> (around Z axis), ω<sub>A</sub> (around X axis)
  2271.      * and χ<sub>A</sub> (around Z axis). The constant angle ε₀
  2272.      * for the fourth rotation (around X axis) can be retrieved by evaluating the
  2273.      * function returned by {@link #getMeanObliquityFunction()} at {@link
  2274.      * #getNutationReferenceEpoch() nutation reference epoch}.
  2275.      * </p>
  2276.      * @return function computing the precession angle
  2277.      * @exception OrekitException if table cannot be loaded
  2278.      * @since 6.1
  2279.      */
  2280.     public abstract TimeVectorFunction getPrecessionFunction() throws OrekitException;

  2282.     /** Get the function computing the nutation angles.
  2283.      * <p>
  2284.      * The function returned computes the two classical angles ΔΨ and Δε,
  2285.      * and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
  2286.      * resolution C7 (the correction is forced to 0 before this date)
  2287.      * </p>
  2288.      * @return function computing the nutation in longitude ΔΨ and Δε
  2289.      * and the correction of equation of equinoxes
  2290.      * @exception OrekitException if table cannot be loaded
  2291.      * @since 6.1
  2292.      */
  2293.     public abstract TimeVectorFunction getNutationFunction()
  2294.         throws OrekitException;

  2296.     /** Get the function computing Greenwich mean sidereal time, in radians.
  2297.      * @param ut1 UT1 time scale
  2298.      * @return function computing Greenwich mean sidereal time
  2299.      * @exception OrekitException if table cannot be loaded
  2300.      * @since 6.1
  2301.      */
  2302.     public abstract TimeScalarFunction getGMSTFunction(TimeScale ut1)
  2303.         throws OrekitException;

  2305.     /** Get the function computing Greenwich mean sidereal time rate, in radians per second.
  2306.      * @param ut1 UT1 time scale
  2307.      * @return function computing Greenwich mean sidereal time rate
  2308.      * @exception OrekitException if table cannot be loaded
  2309.      * @since 9.0
  2310.      */
  2311.     public abstract TimeScalarFunction getGMSTRateFunction(TimeScale ut1)
  2312.         throws OrekitException;

  2314.     /** Get the function computing Greenwich apparent sidereal time, in radians.
  2315.      * @param ut1 UT1 time scale
  2316.      * @param eopHistory EOP history
  2317.      * @return function computing Greenwich apparent sidereal time
  2318.      * @exception OrekitException if table cannot be loaded
  2319.      * @since 6.1
  2320.      */
  2321.     public abstract TimeScalarFunction getGASTFunction(TimeScale ut1, EOPHistory eopHistory)
  2322.         throws OrekitException;

  2324.     /** Get the function computing tidal corrections for Earth Orientation Parameters.
  2325.      * @return function computing tidal corrections for Earth Orientation Parameters,
  2326.      * for xp, yp, ut1 and lod respectively
  2327.      * @exception OrekitException if table cannot be loaded
  2328.      * @since 6.1
  2329.      */
  2330.     public abstract TimeVectorFunction getEOPTidalCorrection()
  2331.         throws OrekitException;

  2333.     /** Get the Love numbers.
  2334.      * @return Love numbers
  2335.      * @exception OrekitException if table cannot be loaded
  2336.      * @since 6.1
  2337.      */
  2338.     public abstract LoveNumbers getLoveNumbers()
  2339.         throws OrekitException;

  2341.     /** Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2342.      * @param ut1 UT1 time scale
  2343.      * @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2344.      * @exception OrekitException if table cannot be loaded
  2345.      * @since 6.1
  2346.      */
  2347.     public abstract TimeVectorFunction getTideFrequencyDependenceFunction(TimeScale ut1)
  2348.         throws OrekitException;

  2350.     /** Get the permanent tide to be <em>removed</em> from ΔC₂₀ when zero-tide potentials are used.
  2351.      * @return permanent tide to remove
  2352.      * @exception OrekitException if table cannot be loaded
  2353.      */
  2354.     public abstract double getPermanentTide() throws OrekitException;

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

  2365.     /** Get the function computing ocean pole tide (ΔC₂₁, ΔS₂₁).
  2366.      * @param eopHistory EOP history
  2367.      * @return model for ocean pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2368.      * @exception OrekitException if table cannot be loaded
  2369.      * @since 6.1
  2370.      */
  2371.     public abstract TimeVectorFunction getOceanPoleTide(EOPHistory eopHistory)
  2372.         throws OrekitException;

  2374.     /** Get the nominal values of the displacement numbers.
  2375.      * @return an array containing h⁽⁰⁾, h⁽²⁾, h₃, hI diurnal, hI semi-diurnal,
  2376.      * l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾, l₃, lI diurnal, lI semi-diurnal,
  2377.      * H₀ permanent deformation amplitude
  2378.      * @since 9.1
  2379.      */
  2380.     public abstract double[] getNominalTidalDisplacement();

  2382.     /** Get the correction function for tidal displacement for diurnal tides.
  2383.      * <ul>
  2384.      *  <li>f[0]: radial correction, longitude cosine part</li>
  2385.      *  <li>f[1]: radial correction, longitude sine part</li>
  2386.      *  <li>f[2]: North correction, longitude cosine part</li>
  2387.      *  <li>f[3]: North correction, longitude sine part</li>
  2388.      *  <li>f[4]: East correction, longitude cosine part</li>
  2389.      *  <li>f[5]: East correction, longitude sine part</li>
  2390.      * </ul>
  2391.      * @return correction function for tidal displacement
  2392.      * @throws OrekitException if Poisson series cannot be loaded
  2393.      * @since 9.1
  2394.      */
  2395.     public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal()
  2396.         throws OrekitException;

  2398.     /** Get the correction function for tidal displacement for diurnal tides.
  2399.      * <ul>
  2400.      *  <li>f[0]: radial correction, longitude cosine part</li>
  2401.      *  <li>f[1]: radial correction, longitude sine part</li>
  2402.      *  <li>f[2]: North correction, longitude cosine part</li>
  2403.      *  <li>f[3]: North correction, longitude sine part</li>
  2404.      *  <li>f[4]: East correction, longitude cosine part</li>
  2405.      *  <li>f[5]: East correction, longitude sine part</li>
  2406.      * </ul>
  2407.      * @param tableName name for the diurnal tides table
  2408.      * @param cols total number of columns of the diurnal tides table
  2409.      * @param rIp column holding ∆Rf(ip) in the diurnal tides table, counting from 1
  2410.      * @param rOp column holding ∆Rf(op) in the diurnal tides table, counting from 1
  2411.      * @param tIp column holding ∆Tf(ip) in the diurnal tides table, counting from 1
  2412.      * @param tOp column holding ∆Tf(op) in the diurnal tides table, counting from 1
  2413.      * @return correction function for tidal displacement for diurnal tides
  2414.      * @exception OrekitException if Poisson series cannot be loaded
  2415.      * @since 9.1
  2416.      */
  2417.     protected static CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal(final String tableName, final int cols,
  2418.                                                                                    final int rIp, final int rOp,
  2419.                                                                                    final int tIp, final int tOp)
  2420.         throws OrekitException {

  2422.         // radial component, missing the sin 2φ factor; this corresponds to:
  2423.         //  - equation 15a in IERS conventions 1996, chapter 7
  2424.         //  - equation 16a in IERS conventions 2003, chapter 7
  2425.         //  - equation 7.12a in IERS conventions 2010, chapter 7
  2426.         final PoissonSeries drCos = new PoissonSeriesParser(cols).
  2427.                                     withOptionalColumn(1).
  2428.                                     withDoodson(4, 3).
  2429.                                     withFirstDelaunay(10).
  2430.                                     withSinCos(0, rIp, +1.0e-3, rOp, +1.0e-3).
  2431.                                     parse(getStream(tableName), tableName);
  2432.         final PoissonSeries drSin = new PoissonSeriesParser(cols).
  2433.                                     withOptionalColumn(1).
  2434.                                     withDoodson(4, 3).
  2435.                                     withFirstDelaunay(10).
  2436.                                     withSinCos(0, rOp, -1.0e-3, rIp, +1.0e-3).
  2437.                                     parse(getStream(tableName), tableName);

  2439.         // North component, missing the cos 2φ factor; this corresponds to:
  2440.         //  - equation 15b in IERS conventions 1996, chapter 7
  2441.         //  - equation 16b in IERS conventions 2003, chapter 7
  2442.         //  - equation 7.12b in IERS conventions 2010, chapter 7
  2443.         final PoissonSeries dnCos = new PoissonSeriesParser(cols).
  2444.                                     withOptionalColumn(1).
  2445.                                     withDoodson(4, 3).
  2446.                                     withFirstDelaunay(10).
  2447.                                     withSinCos(0, tIp, +1.0e-3, tOp, +1.0e-3).
  2448.                                     parse(getStream(tableName), tableName);
  2449.         final PoissonSeries dnSin = new PoissonSeriesParser(cols).
  2450.                                     withOptionalColumn(1).
  2451.                                     withDoodson(4, 3).
  2452.                                     withFirstDelaunay(10).
  2453.                                     withSinCos(0, tOp, -1.0e-3, tIp, +1.0e-3).
  2454.                                     parse(getStream(tableName), tableName);

  2456.         // East component, missing the sin φ factor; this corresponds to:
  2457.         //  - equation 15b in IERS conventions 1996, chapter 7
  2458.         //  - equation 16b in IERS conventions 2003, chapter 7
  2459.         //  - equation 7.12b in IERS conventions 2010, chapter 7
  2460.         final PoissonSeries deCos = new PoissonSeriesParser(cols).
  2461.                                     withOptionalColumn(1).
  2462.                                     withDoodson(4, 3).
  2463.                                     withFirstDelaunay(10).
  2464.                                     withSinCos(0, tOp, -1.0e-3, tIp, +1.0e-3).
  2465.                                     parse(getStream(tableName), tableName);
  2466.         final PoissonSeries deSin = new PoissonSeriesParser(cols).
  2467.                                     withOptionalColumn(1).
  2468.                                     withDoodson(4, 3).
  2469.                                     withFirstDelaunay(10).
  2470.                                     withSinCos(0, tIp, -1.0e-3, tOp, -1.0e-3).
  2471.                                     parse(getStream(tableName), tableName);

  2473.         return PoissonSeries.compile(drCos, drSin,
  2474.                                      dnCos, dnSin,
  2475.                                      deCos, deSin);

  2477.     }

  2479.     /** Get the correction function for tidal displacement for zonal tides.
  2480.      * <ul>
  2481.      *  <li>f[0]: radial correction</li>
  2482.      *  <li>f[1]: North correction</li>
  2483.      * </ul>
  2484.      * @return correction function for tidal displacement
  2485.      * @throws OrekitException if Poisson series cannot be loaded
  2486.      * @since 9.1
  2487.      */
  2488.     public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionZonal()
  2489.         throws OrekitException;

  2491.     /** Get the correction function for tidal displacement for zonal tides.
  2492.      * <ul>
  2493.      *  <li>f[0]: radial correction</li>
  2494.      *  <li>f[1]: North correction</li>
  2495.      * </ul>
  2496.      * @param tableName name for the zonal tides table
  2497.      * @param cols total number of columns of the table
  2498.      * @param rIp column holding ∆Rf(ip) in the table, counting from 1
  2499.      * @param rOp column holding ∆Rf(op) in the table, counting from 1
  2500.      * @param tIp column holding ∆Tf(ip) in the table, counting from 1
  2501.      * @param tOp column holding ∆Tf(op) in the table, counting from 1
  2502.      * @return correction function for tidal displacement for zonal tides
  2503.      * @exception OrekitException if Poisson series cannot be loaded
  2504.      * @since 9.1
  2505.      */
  2506.     protected static CompiledSeries getTidalDisplacementFrequencyCorrectionZonal(final String tableName, final int cols,
  2507.                                                                                  final int rIp, final int rOp,
  2508.                                                                                  final int tIp, final int tOp)
  2509.         throws OrekitException {

  2511.         // radial component, missing the 3⁄2 sin² φ - 1⁄2 factor; this corresponds to:
  2512.         //  - equation 16a in IERS conventions 1996, chapter 7
  2513.         //  - equation 17a in IERS conventions 2003, chapter 7
  2514.         //  - equation 7.13a in IERS conventions 2010, chapter 7
  2515.         final PoissonSeries dr = new PoissonSeriesParser(cols).
  2516.                                  withOptionalColumn(1).
  2517.                                  withDoodson(4, 3).
  2518.                                  withFirstDelaunay(10).
  2519.                                  withSinCos(0, rOp, +1.0e-3, rIp, +1.0e-3).
  2520.                                  parse(getStream(tableName), tableName);

  2522.         // North component, missing the sin 2φ factor; this corresponds to:
  2523.         //  - equation 16b in IERS conventions 1996, chapter 7
  2524.         //  - equation 17b in IERS conventions 2003, chapter 7
  2525.         //  - equation 7.13b in IERS conventions 2010, chapter 7
  2526.         final PoissonSeries dn = new PoissonSeriesParser(cols).
  2527.                                  withOptionalColumn(1).
  2528.                                  withDoodson(4, 3).
  2529.                                  withFirstDelaunay(10).
  2530.                                  withSinCos(0, tOp, +1.0e-3, tIp, +1.0e-3).
  2531.                                  parse(getStream(tableName), tableName);

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

  2535.     }

  2537.     /** Interface for functions converting nutation corrections between
  2538.      * δΔψ/δΔε to δX/δY.
  2539.      * <ul>
  2540.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  2541.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  2542.      * </ul>
  2543.      * @since 6.1
  2544.      */
  2545.     public interface NutationCorrectionConverter {

  2547.         /** Convert nutation corrections.
  2548.          * @param date current date
  2549.          * @param ddPsi δΔψ part of the nutation correction
  2550.          * @param ddEpsilon δΔε part of the nutation correction
  2551.          * @return array containing δX and δY
  2552.          * @exception OrekitException if correction cannot be converted
  2553.          */
  2554.         double[] toNonRotating(AbsoluteDate date, double ddPsi, double ddEpsilon)
  2555.             throws OrekitException;

  2557.         /** Convert nutation corrections.
  2558.          * @param date current date
  2559.          * @param dX δX part of the nutation correction
  2560.          * @param dY δY part of the nutation correction
  2561.          * @return array containing δΔψ and δΔε
  2562.          * @exception OrekitException if correction cannot be converted
  2563.          */
  2564.         double[] toEquinox(AbsoluteDate date, double dX, double dY)
  2565.             throws OrekitException;

  2567.     }

  2569.     /** Create a function converting nutation corrections between
  2570.      * δX/δY and δΔψ/δΔε.
  2571.      * <ul>
  2572.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  2573.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  2574.      * </ul>
  2575.      * @return a new converter
  2576.      * @exception OrekitException if some convention table cannot be loaded
  2577.      * @since 6.1
  2578.      */
  2579.     public NutationCorrectionConverter getNutationCorrectionConverter()
  2580.         throws OrekitException {

  2582.         // get models parameters
  2583.         final TimeVectorFunction precessionFunction = getPrecessionFunction();
  2584.         final TimeScalarFunction epsilonAFunction = getMeanObliquityFunction();
  2585.         final AbsoluteDate date0 = getNutationReferenceEpoch();
  2586.         final double cosE0 = FastMath.cos(epsilonAFunction.value(date0));

  2588.         return new NutationCorrectionConverter() {

  2590.             /** {@inheritDoc} */
  2591.             @Override
  2592.             public double[] toNonRotating(final AbsoluteDate date,
  2593.                                           final double ddPsi, final double ddEpsilon)
  2594.                 throws OrekitException {
  2595.                 // compute precession angles psiA, omegaA and chiA
  2596.                 final double[] angles = precessionFunction.value(date);

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

  2602.                 // convert nutation corrections (equation 23/IERS-2003 or 5.25/IERS-2010)
  2603.                 return new double[] {
  2604.                     sinEA * ddPsi + c * ddEpsilon,
  2605.                     ddEpsilon - c * sinEA * ddPsi
  2606.                 };

  2608.             }

  2610.             /** {@inheritDoc} */
  2611.             @Override
  2612.             public double[] toEquinox(final AbsoluteDate date,
  2613.                                       final double dX, final double dY)
  2614.                 throws OrekitException {
  2615.                 // compute precession angles psiA, omegaA and chiA
  2616.                 final double[] angles   = precessionFunction.value(date);

  2618.                 // conversion coefficients
  2619.                 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
  2620.                 final double c     = angles[0] * cosE0 - angles[2];
  2621.                 final double opC2  = 1 + c * c;

  2623.                 // convert nutation corrections (inverse of equation 23/IERS-2003 or 5.25/IERS-2010)
  2624.                 return new double[] {
  2625.                     (dX - c * dY) / (sinEA * opC2),
  2626.                     (dY + c * dX) / opC2
  2627.                 };
  2628.             }

  2630.         };

  2632.     }

  2634.     /** Load the Love numbers.
  2635.      * @param nameLove name of the Love number resource
  2636.      * @return Love numbers
  2637.      * @exception OrekitException if the Love numbers embedded in the
  2638.      * library cannot be read
  2639.      */
  2640.     protected LoveNumbers loadLoveNumbers(final String nameLove) throws OrekitException {
  2641.         try {

  2643.             // allocate the three triangular arrays for real, imaginary and time-dependent numbers
  2644.             final double[][] real      = new double[4][];
  2645.             final double[][] imaginary = new double[4][];
  2646.             final double[][] plus      = new double[4][];
  2647.             for (int i = 0; i < real.length; ++i) {
  2648.                 real[i]      = new double[i + 1];
  2649.                 imaginary[i] = new double[i + 1];
  2650.                 plus[i]      = new double[i + 1];
  2651.             }

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

  2655.                 if (stream == null) {
  2656.                     // this should never happen with files embedded within Orekit
  2657.                     throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, nameLove);
  2658.                 }

  2660.                 // setup the reader
  2661.                 try (BufferedReader reader = new BufferedReader(new InputStreamReader(stream, "UTF-8"))) {

  2663.                     String line = reader.readLine();
  2664.                     int lineNumber = 1;

  2666.                     // look for the Love numbers
  2667.                     while (line != null) {

  2669.                         line = line.trim();
  2670.                         if (!(line.isEmpty() || line.startsWith("#"))) {
  2671.                             final String[] fields = line.split("\\p{Space}+");
  2672.                             if (fields.length != 5) {
  2673.                                 // this should never happen with files embedded within Orekit
  2674.                                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  2675.                                                           lineNumber, nameLove, line);
  2676.                             }
  2677.                             final int n = Integer.parseInt(fields[0]);
  2678.                             final int m = Integer.parseInt(fields[1]);
  2679.                             if ((n < 2) || (n > 3) || (m < 0) || (m > n)) {
  2680.                                 // this should never happen with files embedded within Orekit
  2681.                                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  2682.                                                           lineNumber, nameLove, line);

  2684.                             }
  2685.                             real[n][m]      = Double.parseDouble(fields[2]);
  2686.                             imaginary[n][m] = Double.parseDouble(fields[3]);
  2687.                             plus[n][m]      = Double.parseDouble(fields[4]);
  2688.                         }

  2690.                         // next line
  2691.                         lineNumber++;
  2692.                         line = reader.readLine();

  2694.                     }
  2695.                 }
  2696.             }

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

  2700.         } catch (IOException ioe) {
  2701.             // this should never happen with files embedded within Orekit
  2702.             throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, nameLove);
  2703.         }
  2704.     }

  2706.     /** Get a data stream.
  2707.      * @param name file name of the resource stream
  2708.      * @return stream
  2709.      */
  2710.     private static InputStream getStream(final String name) {
  2711.         return IERSConventions.class.getResourceAsStream(name);
  2712.     }

  2714.     /** Correction to equation of equinoxes.
  2715.      * <p>IAU 1994 resolution C7 added two terms to the equation of equinoxes
  2716.      * taking effect since 1997-02-27 for continuity
  2717.      * </p>
  2718.      */
  2719.     private static class IAU1994ResolutionC7 {

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

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

  2727.         /** Start date for applying Moon corrections to the equation of the equinoxes.
  2728.          * This date corresponds to 1997-02-27T00:00:00 UTC, hence the 30s offset from TAI.
  2729.          */
  2730.         private static final AbsoluteDate MODEL_START =
  2731.             new AbsoluteDate(1997, 2, 27, 0, 0, 30, TimeScalesFactory.getTAI());

  2733.         /** Evaluate the correction.
  2734.          * @param arguments Delaunay for nutation
  2735.          * @return correction value (0 before 1997-02-27)
  2736.          */
  2737.         public static double value(final DelaunayArguments arguments) {
  2738.             if (arguments.getDate().compareTo(MODEL_START) >= 0) {

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

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

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

  2749.             } else {
  2750.                 return 0.0;
  2751.             }
  2752.         }

  2754.         /** Evaluate the correction.
  2755.          * @param arguments Delaunay for nutation
  2756.          * @param <T> type of the field elements
  2757.          * @return correction value (0 before 1997-02-27)
  2758.          */
  2759.         public static <T extends RealFieldElement<T>> T value(final FieldDelaunayArguments<T> arguments) {
  2760.             if (arguments.getDate().toAbsoluteDate().compareTo(MODEL_START) >= 0) {

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

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

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

  2771.             } else {
  2772.                 return arguments.getDate().getField().getZero();
  2773.             }
  2774.         }

  2776.     };

  2778.     /** Stellar angle model.
  2779.      * <p>
  2780.      * The stellar angle computed here has been defined in the paper "A non-rotating origin on the
  2781.      * instantaneous equator: Definition, properties and use", N. Capitaine, Guinot B., and Souchay J.,
  2782.      * Celestial Mechanics, Volume 39, Issue 3, pp 283-307. It has been proposed as a conventional
  2783.      * conventional relationship between UT1 and Earth rotation in the paper "Definition of the Celestial
  2784.      * Ephemeris origin and of UT1 in the International Celestial Reference Frame", Capitaine, N.,
  2785.      * Guinot, B., and McCarthy, D. D., 2000, Astronomy and Astrophysics, 355(1), pp. 398–405.
  2786.      * </p>
  2787.      * <p>
  2788.      * It is presented simply as stellar angle in IERS conventions 1996 but as since been adopted as
  2789.      * the conventional relationship defining UT1 from ICRF and is called Earth Rotation Angle in
  2790.      * IERS conventions 2003 and 2010.
  2791.      * </p>
  2792.      */
  2793.     private static class StellarAngleCapitaine implements TimeScalarFunction {

  2795.         /** Reference date of Capitaine's Earth Rotation Angle model. */
  2796.         private static final AbsoluteDate REFERENCE_DATE = new AbsoluteDate(DateComponents.J2000_EPOCH,
  2797.                                                                             TimeComponents.H12,
  2798.                                                                             TimeScalesFactory.getTAI());

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

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

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

  2811.         /** UT1 time scale. */
  2812.         private final TimeScale ut1;

  2814.         /** Simple constructor.
  2815.          * @param ut1 UT1 time scale
  2816.          */
  2817.         StellarAngleCapitaine(final TimeScale ut1) {
  2818.             this.ut1 = ut1;
  2819.         }

  2821.         /** Get the rotation rate.
  2822.          * @return rotation rate
  2823.          */
  2824.         public double getRate() {
  2825.             return ERA_1A + ERA_1B;
  2826.         }

  2828.         /** {@inheritDoc} */
  2829.         @Override
  2830.         public double value(final AbsoluteDate date) {

  2832.             // split the date offset as a full number of days plus a smaller part
  2833.             final int secondsInDay = 86400;
  2834.             final double dt  = date.durationFrom(REFERENCE_DATE);
  2835.             final long days  = ((long) dt) / secondsInDay;
  2836.             final double dtA = secondsInDay * days;
  2837.             final double dtB = (dt - dtA) + ut1.offsetFromTAI(date);

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

  2841.         }

  2843.         /** {@inheritDoc} */
  2844.         @Override
  2845.         public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  2847.             // split the date offset as a full number of days plus a smaller part
  2848.             final int secondsInDay = 86400;
  2849.             final T dt  = date.durationFrom(REFERENCE_DATE);
  2850.             final long days  = ((long) dt.getReal()) / secondsInDay;
  2851.             final double dtA = secondsInDay * days;
  2852.             final T dtB = dt.subtract(dtA).add(ut1.offsetFromTAI(date.toAbsoluteDate()));

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

  2856.         }

  2858.     }

  2860.     /** Mean pole. */
  2861.     private static class MeanPole implements TimeStamped, Serializable {

  2863.         /** Serializable UID. */
  2864.         private static final long serialVersionUID = 20131028l;

  2866.         /** Date. */
  2867.         private final AbsoluteDate date;

  2869.         /** X coordinate. */
  2870.         private double x;

  2872.         /** Y coordinate. */
  2873.         private double y;

  2875.         /** Simple constructor.
  2876.          * @param date date
  2877.          * @param x x coordinate
  2878.          * @param y y coordinate
  2879.          */
  2880.         MeanPole(final AbsoluteDate date, final double x, final double y) {
  2881.             this.date = date;
  2882.             this.x    = x;
  2883.             this.y    = y;
  2884.         }

  2886.         /** {@inheritDoc} */
  2887.         @Override
  2888.         public AbsoluteDate getDate() {
  2889.             return date;
  2890.         }

  2892.         /** Get x coordinate.
  2893.          * @return x coordinate
  2894.          */
  2895.         public double getX() {
  2896.             return x;
  2897.         }

  2899.         /** Get y coordinate.
  2900.          * @return y coordinate
  2901.          */
  2902.         public double getY() {
  2903.             return y;
  2904.         }

  2906.     }

  2908.     /** Local class for precession function. */
  2909.     private class PrecessionFunction implements TimeVectorFunction {

  2911.         /** Polynomial nutation for psiA. */
  2912.         private final PolynomialNutation psiA;

  2914.         /** Polynomial nutation for omegaA. */
  2915.         private final PolynomialNutation omegaA;

  2917.         /** Polynomial nutation for chiA. */
  2918.         private final PolynomialNutation chiA;

  2920.         /** Simple constructor.
  2921.          * @param psiA polynomial nutation for psiA
  2922.          * @param omegaA polynomial nutation for omegaA
  2923.          * @param chiA polynomial nutation for chiA
  2924.          */
  2925.         PrecessionFunction(final PolynomialNutation psiA,
  2926.                            final PolynomialNutation omegaA,
  2927.                            final PolynomialNutation chiA) {
  2928.             this.psiA   = psiA;
  2929.             this.omegaA = omegaA;
  2930.             this.chiA   = chiA;
  2931.         }


  2934.         /** {@inheritDoc} */
  2935.         @Override
  2936.         public double[] value(final AbsoluteDate date) {
  2937.             final double tc = evaluateTC(date);
  2938.             return new double[] {
  2939.                 psiA.value(tc), omegaA.value(tc), chiA.value(tc)
  2940.             };
  2941.         }

  2943.         /** {@inheritDoc} */
  2944.         @Override
  2945.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  2946.             final T[] a = MathArrays.buildArray(date.getField(), 3);
  2947.             final T tc = evaluateTC(date);
  2948.             a[0] = psiA.value(tc);
  2949.             a[1] = omegaA.value(tc);
  2950.             a[2] = chiA.value(tc);
  2951.             return a;
  2952.         }

  2954.     }

  2956.     /** Local class for tides frequency function. */
  2957.     private static class TideFrequencyDependenceFunction implements TimeVectorFunction {

  2959.         /** Nutation arguments. */
  2960.         private final FundamentalNutationArguments arguments;

  2962.         /** Correction series. */
  2963.         private final PoissonSeries.CompiledSeries kSeries;

  2965.         /** Simple constructor.
  2966.          * @param arguments nutation arguments
  2967.          * @param c20Series correction series for the C20 term
  2968.          * @param c21Series correction series for the C21 term
  2969.          * @param s21Series correction series for the S21 term
  2970.          * @param c22Series correction series for the C22 term
  2971.          * @param s22Series correction series for the S22 term
  2972.          */
  2973.         TideFrequencyDependenceFunction(final FundamentalNutationArguments arguments,
  2974.                                         final PoissonSeries c20Series,
  2975.                                         final PoissonSeries c21Series,
  2976.                                         final PoissonSeries s21Series,
  2977.                                         final PoissonSeries c22Series,
  2978.                                         final PoissonSeries s22Series) {
  2979.             this.arguments = arguments;
  2980.             this.kSeries   = PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);
  2981.         }


  2984.         /** {@inheritDoc} */
  2985.         @Override
  2986.         public double[] value(final AbsoluteDate date) {
  2987.             return kSeries.value(arguments.evaluateAll(date));
  2988.         }

  2990.         /** {@inheritDoc} */
  2991.         @Override
  2992.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  2993.             return kSeries.value(arguments.evaluateAll(date));
  2994.         }

  2996.     }

  2998.     /** Local class for EOP tidal corrections. */
  2999.     private static class EOPTidalCorrection implements TimeVectorFunction {

  3001.         /** Nutation arguments. */
  3002.         private final FundamentalNutationArguments arguments;

  3004.         /** Correction series. */
  3005.         private final PoissonSeries.CompiledSeries correctionSeries;

  3007.         /** Simple constructor.
  3008.          * @param arguments nutation arguments
  3009.          * @param xSeries correction series for the x coordinate
  3010.          * @param ySeries correction series for the y coordinate
  3011.          * @param ut1Series correction series for the UT1
  3012.          */
  3013.         EOPTidalCorrection(final FundamentalNutationArguments arguments,
  3014.                            final PoissonSeries xSeries,
  3015.                            final PoissonSeries ySeries,
  3016.                            final PoissonSeries ut1Series) {
  3017.             this.arguments        = arguments;
  3018.             this.correctionSeries = PoissonSeries.compile(xSeries, ySeries, ut1Series);
  3019.         }

  3021.         /** {@inheritDoc} */
  3022.         @Override
  3023.         public double[] value(final AbsoluteDate date) {
  3024.             final BodiesElements elements = arguments.evaluateAll(date);
  3025.             final double[] correction    = correctionSeries.value(elements);
  3026.             final double[] correctionDot = correctionSeries.derivative(elements);
  3027.             return new double[] {
  3028.                 correction[0],
  3029.                 correction[1],
  3030.                 correction[2],
  3031.                 correctionDot[2] * (-Constants.JULIAN_DAY)
  3032.             };
  3033.         }

  3035.         /** {@inheritDoc} */
  3036.         @Override
  3037.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

  3039.             final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  3040.             final T[] correction    = correctionSeries.value(elements);
  3041.             final T[] correctionDot = correctionSeries.derivative(elements);
  3042.             final T[] a = MathArrays.buildArray(date.getField(), 4);
  3043.             a[0] = correction[0];
  3044.             a[1] = correction[1];
  3045.             a[2] = correction[2];
  3046.             a[3] = correctionDot[2].multiply(-Constants.JULIAN_DAY);
  3047.             return a;
  3048.         }

  3050.     }

  3052. }
Created with JaCoCo 0.7.9.201702052155