IERSConventions.java

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  133.             return new TimeScalarFunction() {

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

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

  147.             };

  149.         }

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

  155.             // set up nutation arguments
  156.             final FundamentalNutationArguments arguments =
  157.                     getNutationArguments(timeScales);

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

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

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

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

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

  200.             final double fYCosOm    = -0.00231 * Constants.ARC_SECONDS_TO_RADIANS;
  201.             final double fYCos2FDOm = -0.00014 * Constants.ARC_SECONDS_TO_RADIANS;

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

  214.             final PoissonSeries.CompiledSeries xySum =
  215.                     PoissonSeries.compile(xSum, ySum);

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

  226.             return new TimeVectorFunction() {

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

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

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

  240.                     final SinCos scOmega   = FastMath.sinCos(omega);
  241.                     final SinCos sc2omega  = SinCos.sum(scOmega, scOmega);
  242.                     final SinCos sc2FD0m   = FastMath.sinCos(2 * (f - d + omega));
  243.                     final double cosOmega  = scOmega.cos();
  244.                     final double sinOmega  = scOmega.sin();
  245.                     final double sin2Omega = sc2omega.sin();
  246.                     final double cos2FDOm  = sc2FD0m.cos();
  247.                     final double sin2FDOm  = sc2FD0m.sin();

  249.                     final double x = xPolynomial.value(t) + sinEps0 * xy[0] +
  250.                             t * t * (fXCosOm * cosOmega + fXSinOm * sinOmega + fXSin2FDOm * cos2FDOm);
  251.                     final double y = yPolynomial.value(t) + xy[1] +
  252.                             t * t * (fYCosOm * cosOmega + fYCos2FDOm * cos2FDOm);
  253.                     final double sPxy2 = fSSinOm * sinOmega + fSSin2Om * sin2Omega +
  254.                             t * (fST + t * (fST2SinOm * sinOmega + fST2Sin2FDOm * sin2FDOm + t * fST3));

  256.                     return new double[] {
  257.                         x, y, sPxy2
  258.                     };

  260.                 }

  262.                 /** {@inheritDoc} */
  263.                 @Override
  264.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

  269.                     final T omega     = elements.getOmega();
  270.                     final T f         = elements.getF();
  271.                     final T d         = elements.getD();
  272.                     final T t         = elements.getTC();
  273.                     final T t2        = t.multiply(t);

  275.                     final FieldSinCos<T> scOmega  = FastMath.sinCos(omega);
  276.                     final FieldSinCos<T> sc2omega = FieldSinCos.sum(scOmega, scOmega);
  277.                     final FieldSinCos<T> sc2FD0m  = FastMath.sinCos(f.subtract(d).add(omega).multiply(2));
  278.                     final T cosOmega  = scOmega.cos();
  279.                     final T sinOmega  = scOmega.sin();
  280.                     final T sin2Omega = sc2omega.sin();
  281.                     final T cos2FDOm  = sc2FD0m.cos();
  282.                     final T sin2FDOm  = sc2FD0m.sin();

  284.                     final T x = xPolynomial.value(t).
  285.                                 add(xy[0].multiply(sinEps0)).
  286.                                 add(t2.multiply(cosOmega.multiply(fXCosOm).add(sinOmega.multiply(fXSinOm)).add(cos2FDOm.multiply(fXSin2FDOm))));
  287.                     final T y = yPolynomial.value(t).
  288.                                 add(xy[1]).
  289.                                 add(t2.multiply(cosOmega.multiply(fYCosOm).add(cos2FDOm.multiply(fYCos2FDOm))));
  290.                     final T sPxy2 = sinOmega.multiply(fSSinOm).
  291.                                     add(sin2Omega.multiply(fSSin2Om)).
  292.                                     add(t.multiply(fST3).add(sinOmega.multiply(fST2SinOm)).add(sin2FDOm.multiply(fST2Sin2FDOm)).multiply(t).add(fST).multiply(t));

  294.                     final T[] a = MathArrays.buildArray(date.getField(), 3);
  295.                     a[0] = x;
  296.                     a[1] = y;
  297.                     a[2] = sPxy2;
  298.                     return a;

  300.                 }

  302.             };

  304.         }

  306.         /** {@inheritDoc} */
  307.         @Override
  308.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

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

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

  334.         }

  336.         /** {@inheritDoc} */
  337.         @Override
  338.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  340.             // set up nutation arguments
  341.             final FundamentalNutationArguments arguments =
  342.                     getNutationArguments(timeScales);

  344.             // set up Poisson series
  345.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  346.             final PoissonSeriesParser baseParser =
  347.                     new PoissonSeriesParser(10).withFirstDelaunay(1);

  349.             final PoissonSeriesParser psiParser =
  350.                     baseParser.
  351.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  352.                     withSinCos(1, 8, deciMilliAS, -1, deciMilliAS);
  353.             final PoissonSeries psiSeries;
  354.             try (InputStream in = getStream(PSI_EPSILON_SERIES)) {
  355.                 psiSeries = psiParser.parse(in, PSI_EPSILON_SERIES);
  356.             } catch (IOException e) {
  357.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  358.             }

  360.             final PoissonSeriesParser epsilonParser =
  361.                     baseParser.
  362.                     withSinCos(0, -1, deciMilliAS, 9, deciMilliAS).
  363.                     withSinCos(1, -1, deciMilliAS, 10, deciMilliAS);
  364.             final PoissonSeries epsilonSeries;
  365.             try (InputStream in = getStream(PSI_EPSILON_SERIES)) {
  366.                 epsilonSeries = epsilonParser.parse(in, PSI_EPSILON_SERIES);
  367.             } catch (IOException e) {
  368.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  369.             }

  371.             final PoissonSeries.CompiledSeries psiEpsilonSeries =
  372.                     PoissonSeries.compile(psiSeries, epsilonSeries);

  374.             return new TimeVectorFunction() {

  376.                 /** {@inheritDoc} */
  377.                 @Override
  378.                 public double[] value(final AbsoluteDate date) {
  379.                     final BodiesElements elements = arguments.evaluateAll(date);
  380.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  381.                     return new double[] {
  382.                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  383.                     };
  384.                 }

  386.                 /** {@inheritDoc} */
  387.                 @Override
  388.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  389.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  390.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  391.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  392.                     result[0] = psiEpsilon[0];
  393.                     result[1] = psiEpsilon[1];
  394.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  395.                     return result;
  396.                 }

  398.             };

  400.         }

  402.         /** {@inheritDoc} */
  403.         @Override
  404.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  405.                                                   final TimeScales timeScales) {

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

  410.             // constants from IERS 1996 page 21
  411.             // the underlying model is IAU 1982 GMST-UT1
  412.             final AbsoluteDate gmstReference = new AbsoluteDate(
  413.                     DateComponents.J2000_EPOCH, TimeComponents.H12, timeScales.getTAI());
  414.             final double gmst0 = 24110.54841;
  415.             final double gmst1 = 8640184.812866;
  416.             final double gmst2 = 0.093104;
  417.             final double gmst3 = -6.2e-6;

  419.             return new TimeScalarFunction() {

  421.                 /** {@inheritDoc} */
  422.                 @Override
  423.                 public double value(final AbsoluteDate date) {

  425.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  426.                     final double dtai = date.durationFrom(gmstReference);
  427.                     final double tut1 = dtai + ut1.offsetFromTAI(date);
  428.                     final double tt   = tut1 / Constants.JULIAN_CENTURY;

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

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

  437.                 }

  439.                 /** {@inheritDoc} */
  440.                 @Override
  441.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  443.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  444.                     final T dtai = date.durationFrom(gmstReference);
  445.                     final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
  446.                     final T tt   = tut1.divide(Constants.JULIAN_CENTURY);

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

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

  455.                 }

  457.             };

  459.         }

  461.         /** {@inheritDoc} */
  462.         @Override
  463.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  464.                                                       final TimeScales timeScales) {

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

  469.             // constants from IERS 1996 page 21
  470.             // the underlying model is IAU 1982 GMST-UT1
  471.             final AbsoluteDate gmstReference = new AbsoluteDate(
  472.                     DateComponents.J2000_EPOCH, TimeComponents.H12, timeScales.getTAI());
  473.             final double gmst1 = 8640184.812866;
  474.             final double gmst2 = 0.093104;
  475.             final double gmst3 = -6.2e-6;

  477.             return new TimeScalarFunction() {

  479.                 /** {@inheritDoc} */
  480.                 @Override
  481.                 public double value(final AbsoluteDate date) {

  483.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  484.                     final double dtai = date.durationFrom(gmstReference);
  485.                     final double tut1 = dtai + ut1.offsetFromTAI(date);
  486.                     final double tt   = tut1 / Constants.JULIAN_CENTURY;

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

  491.                 }

  493.                 /** {@inheritDoc} */
  494.                 @Override
  495.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  497.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  498.                     final T dtai = date.durationFrom(gmstReference);
  499.                     final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
  500.                     final T tt   = tut1.divide(Constants.JULIAN_CENTURY);

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

  505.                 }

  507.             };

  509.         }

  511.         /** {@inheritDoc} */
  512.         @Override
  513.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  514.                                                   final EOPHistory eopHistory,
  515.                                                   final TimeScales timeScales) {

  517.             // obliquity
  518.             final TimeScalarFunction epsilonA = getMeanObliquityFunction(timeScales);

  520.             // GMST function
  521.             final TimeScalarFunction gmst = getGMSTFunction(ut1, timeScales);

  523.             // nutation function
  524.             final TimeVectorFunction nutation = getNutationFunction(timeScales);

  526.             return new TimeScalarFunction() {

  528.                 /** {@inheritDoc} */
  529.                 @Override
  530.                 public double value(final AbsoluteDate date) {

  532.                     // compute equation of equinoxes
  533.                     final double[] angles = nutation.value(date);
  534.                     double deltaPsi = angles[0];
  535.                     if (eopHistory != null) {
  536.                         deltaPsi += eopHistory.getEquinoxNutationCorrection(date)[0];
  537.                     }
  538.                     final double eqe = deltaPsi  * FastMath.cos(epsilonA.value(date)) + angles[2];

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

  543.                 }

  545.                 /** {@inheritDoc} */
  546.                 @Override
  547.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  549.                     // compute equation of equinoxes
  550.                     final T[] angles = nutation.value(date);
  551.                     T deltaPsi = angles[0];
  552.                     if (eopHistory != null) {
  553.                         deltaPsi = deltaPsi.add(eopHistory.getEquinoxNutationCorrection(date)[0]);
  554.                     }
  555.                     final T eqe = deltaPsi.multiply(epsilonA.value(date).cos()).add(angles[2]);

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

  560.                 }

  562.             };

  564.         }

  566.         /** {@inheritDoc} */
  567.         @Override
  568.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

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

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

  601.             final double deciMilliS = 1.0e-4;
  602.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(17).
  603.                     withOptionalColumn(1).
  604.                     withGamma(7).
  605.                     withFirstDelaunay(2).
  606.                     withSinCos(0, 16, deciMilliS, 17, deciMilliS);
  607.             final PoissonSeries ut1Series;
  608.             try (InputStream in = getStream(TIDAL_CORRECTION_UT1_SERIES)) {
  609.                 ut1Series = ut1Parser.parse(in, TIDAL_CORRECTION_UT1_SERIES);
  610.             } catch (IOException e) {
  611.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  612.             }

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

  616.         }

  618.         /** {@inheritDoc} */
  619.         public LoveNumbers getLoveNumbers() {
  620.             return loadLoveNumbers(LOVE_NUMBERS);
  621.         }

  623.         /** {@inheritDoc} */
  624.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  625.                                                                      final TimeScales timeScales) {

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

  630.             // set up Poisson series
  631.             final PoissonSeriesParser k20Parser =
  632.                     new PoissonSeriesParser(18).
  633.                         withOptionalColumn(1).
  634.                         withDoodson(4, 2).
  635.                         withFirstDelaunay(10);
  636.             final PoissonSeriesParser k21Parser =
  637.                     new PoissonSeriesParser(18).
  638.                         withOptionalColumn(1).
  639.                         withDoodson(4, 3).
  640.                         withFirstDelaunay(10);
  641.             final PoissonSeriesParser k22Parser =
  642.                     new PoissonSeriesParser(16).
  643.                         withOptionalColumn(1).
  644.                         withDoodson(4, 2).
  645.                         withFirstDelaunay(10);

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

  674.                 return new TideFrequencyDependenceFunction(arguments,
  675.                                                            c20Series,
  676.                                                            c21Series, s21Series,
  677.                                                            c22Series, s22Series);
  678.             } catch (IOException e) {
  679.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  680.             }

  682.         }

  684.         /** {@inheritDoc} */
  685.         @Override
  686.         public double getPermanentTide() {
  687.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  688.         }

  690.         /** {@inheritDoc} */
  691.         @Override
  692.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  694.             // constants from IERS 1996 page 47
  695.             final double globalFactor = -1.348e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  696.             final double coupling     =  0.00112;

  698.             return new TimeVectorFunction() {

  700.                 /** {@inheritDoc} */
  701.                 @Override
  702.                 public double[] value(final AbsoluteDate date) {
  703.                     final PoleCorrection pole = eopHistory.getPoleCorrection(date);
  704.                     return new double[] {
  705.                         globalFactor * (pole.getXp() + coupling * pole.getYp()),
  706.                         globalFactor * (coupling * pole.getXp() - pole.getYp()),
  707.                     };
  708.                 }

  710.                 /** {@inheritDoc} */
  711.                 @Override
  712.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  713.                     final FieldPoleCorrection<T> pole = eopHistory.getPoleCorrection(date);
  714.                     final T[] a = MathArrays.buildArray(date.getField(), 2);
  715.                     a[0] = pole.getXp().add(pole.getYp().multiply(coupling)).multiply(globalFactor);
  716.                     a[1] = pole.getXp().multiply(coupling).subtract(pole.getYp()).multiply(globalFactor);
  717.                     return a;
  718.                 }

  720.             };
  721.         }

  723.         /** {@inheritDoc} */
  724.         @Override
  725.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {

  727.             return new TimeVectorFunction() {

  729.                 /** {@inheritDoc} */
  730.                 @Override
  731.                 public double[] value(final AbsoluteDate date) {
  732.                     // there are no model for ocean pole tide prior to conventions 2010
  733.                     return new double[] {
  734.                         0.0, 0.0
  735.                     };
  736.                 }

  738.                 /** {@inheritDoc} */
  739.                 @Override
  740.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  741.                     // there are no model for ocean pole tide prior to conventions 2010
  742.                     return MathArrays.buildArray(date.getField(), 2);
  743.                 }

  745.             };
  746.         }

  748.         /** {@inheritDoc} */
  749.         @Override
  750.         public double[] getNominalTidalDisplacement() {

  752.             //  // elastic Earth values
  753.             //  return new double[] {
  754.             //      // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  755.             //      0.6026,                                  -0.0006, 0.292, -0.0025, -0.0022,
  756.             //      // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  757.             //      0.0831,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  758.             //      // H₀
  759.             //      -0.31460
  760.             //  };

  762.             // anelastic Earth values
  763.             return new double[] {
  764.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  765.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  766.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  767.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  768.                 // H₀
  769.                 -0.31460
  770.             };

  772.         }

  774.         /** {@inheritDoc} */
  775.         @Override
  776.         public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
  777.             return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
  778.                                                                   18, 17, -1, 18, -1);
  779.         }

  781.         /** {@inheritDoc} */
  782.         @Override
  783.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  784.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  785.                                                                 20, 17, 19, 18, 20);
  786.         }

  788.     },

  790.     /** Constant for IERS 2003 conventions. */
  791.     IERS_2003 {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  841.         /** {@inheritDoc} */
  842.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
  843.                                                                  final TimeScales timeScales) {
  844.             try (InputStream in = getStream(NUTATION_ARGUMENTS)) {
  845.                 return new FundamentalNutationArguments(this, timeScale,
  846.                         in, NUTATION_ARGUMENTS, timeScales);
  847.             } catch (IOException e) {
  848.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  849.             }
  850.         }

  852.         /** {@inheritDoc} */
  853.         @Override
  854.         public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {

  856.             // epsilon 0 value from chapter 5, page 41, other terms from equation 32 page 45
  857.             final PolynomialNutation epsilonA =
  858.                     new PolynomialNutation(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  859.                                              -46.84024  * Constants.ARC_SECONDS_TO_RADIANS,
  860.                                               -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  861.                                                0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  863.             return new TimeScalarFunction() {

  865.                 /** {@inheritDoc} */
  866.                 @Override
  867.                 public double value(final AbsoluteDate date) {
  868.                     return epsilonA.value(evaluateTC(date, timeScales));
  869.                 }

  871.                 /** {@inheritDoc} */
  872.                 @Override
  873.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  874.                     return epsilonA.value(evaluateTC(date, timeScales));
  875.                 }

  877.             };

  879.         }

  881.         /** {@inheritDoc} */
  882.         @Override
  883.         public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {

  885.             // set up nutation arguments
  886.             final FundamentalNutationArguments arguments =
  887.                     getNutationArguments(timeScales);

  889.             // set up Poisson series
  890.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  891.             final PoissonSeriesParser parser =
  892.                     new PoissonSeriesParser(17).
  893.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  894.                         withFirstDelaunay(4).
  895.                         withFirstPlanetary(9).
  896.                         withSinCos(0, 2, microAS, 3, microAS);

  898.             final PoissonSeries.CompiledSeries xys;
  899.             try (InputStream xIn = getStream(X_SERIES);
  900.                  InputStream yIn = getStream(Y_SERIES);
  901.                  InputStream sIn = getStream(S_SERIES)) {
  902.                 final PoissonSeries xSeries = parser.parse(xIn, X_SERIES);
  903.                 final PoissonSeries ySeries = parser.parse(yIn, Y_SERIES);
  904.                 final PoissonSeries sSeries = parser.parse(sIn, S_SERIES);
  905.                 xys = PoissonSeries.compile(xSeries, ySeries, sSeries);
  906.             } catch (IOException e) {
  907.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  908.             }

  910.             // create a function evaluating the series
  911.             return new TimeVectorFunction() {

  913.                 /** {@inheritDoc} */
  914.                 @Override
  915.                 public double[] value(final AbsoluteDate date) {
  916.                     return xys.value(arguments.evaluateAll(date));
  917.                 }

  919.                 /** {@inheritDoc} */
  920.                 @Override
  921.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  922.                     return xys.value(arguments.evaluateAll(date));
  923.                 }

  925.             };

  927.         }


  930.         /** {@inheritDoc} */
  931.         @Override
  932.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

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

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

  955.         }

  957.         /** {@inheritDoc} */
  958.         @Override
  959.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  961.             // set up nutation arguments
  962.             final FundamentalNutationArguments arguments =
  963.                     getNutationArguments(timeScales);

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

  990.             final PoissonSeriesParser planetaryParser =
  991.                     new PoissonSeriesParser(21).
  992.                         withFirstDelaunay(2).
  993.                         withFirstPlanetary(7);
  994.             final PoissonSeriesParser planetaryPsiParser =
  995.                     planetaryParser.withSinCos(0, 17, milliAS, 18, milliAS);
  996.             final PoissonSeries psiPlanetarySeries;
  997.             try (InputStream in = getStream(PLANETARY_SERIES)) {
  998.                 psiPlanetarySeries = planetaryPsiParser.parse(in, PLANETARY_SERIES);
  999.             } catch (IOException e) {
  1000.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1001.             }
  1002.             final PoissonSeriesParser planetaryEpsilonParser =
  1003.                     planetaryParser.withSinCos(0, 19, milliAS, 20, milliAS);
  1004.             final PoissonSeries epsilonPlanetarySeries;
  1005.             try (InputStream in = getStream(PLANETARY_SERIES)) {
  1006.                 epsilonPlanetarySeries = planetaryEpsilonParser.parse(in, PLANETARY_SERIES);
  1007.             } catch (IOException e) {
  1008.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1009.             }

  1011.             final PoissonSeries.CompiledSeries luniSolarSeries =
  1012.                     PoissonSeries.compile(psiLuniSolarSeries, epsilonLuniSolarSeries);
  1013.             final PoissonSeries.CompiledSeries planetarySeries =
  1014.                     PoissonSeries.compile(psiPlanetarySeries, epsilonPlanetarySeries);

  1016.             return new TimeVectorFunction() {

  1018.                 /** {@inheritDoc} */
  1019.                 @Override
  1020.                 public double[] value(final AbsoluteDate date) {
  1021.                     final BodiesElements elements = arguments.evaluateAll(date);
  1022.                     final double[] luniSolar = luniSolarSeries.value(elements);
  1023.                     final double[] planetary = planetarySeries.value(elements);
  1024.                     return new double[] {
  1025.                         luniSolar[0] + planetary[0], luniSolar[1] + planetary[1],
  1026.                         IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  1027.                     };
  1028.                 }

  1030.                 /** {@inheritDoc} */
  1031.                 @Override
  1032.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1033.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1034.                     final T[] luniSolar = luniSolarSeries.value(elements);
  1035.                     final T[] planetary = planetarySeries.value(elements);
  1036.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  1037.                     result[0] = luniSolar[0].add(planetary[0]);
  1038.                     result[1] = luniSolar[1].add(planetary[1]);
  1039.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  1040.                     return result;
  1041.                 }

  1043.             };

  1045.         }

  1047.         /** {@inheritDoc} */
  1048.         @Override
  1049.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  1050.                                                   final TimeScales timeScales) {

  1052.             // Earth Rotation Angle
  1053.             final StellarAngleCapitaine era =
  1054.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

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

  1072.             // create a function evaluating the series
  1073.             return new TimeScalarFunction() {

  1075.                 /** {@inheritDoc} */
  1076.                 @Override
  1077.                 public double value(final AbsoluteDate date) {
  1078.                     return era.value(date) + minusEO.value(evaluateTC(date, timeScales));
  1079.                 }

  1081.                 /** {@inheritDoc} */
  1082.                 @Override
  1083.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1084.                     return era.value(date).add(minusEO.value(evaluateTC(date, timeScales)));
  1085.                 }

  1087.             };

  1089.         }

  1091.         /** {@inheritDoc} */
  1092.         @Override
  1093.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  1094.                                                       final TimeScales timeScales) {

  1096.             // Earth Rotation Angle
  1097.             final StellarAngleCapitaine era =
  1098.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

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

  1116.             // create a function evaluating the series
  1117.             return new TimeScalarFunction() {

  1119.                 /** {@inheritDoc} */
  1120.                 @Override
  1121.                 public double value(final AbsoluteDate date) {
  1122.                     return era.getRate() + minusEO.derivative(evaluateTC(date, timeScales));
  1123.                 }

  1125.                 /** {@inheritDoc} */
  1126.                 @Override
  1127.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1128.                     return minusEO.derivative(evaluateTC(date, timeScales)).add(era.getRate());
  1129.                 }

  1131.             };

  1133.         }

  1135.         /** {@inheritDoc} */
  1136.         @Override
  1137.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  1138.                                                   final EOPHistory eopHistory,
  1139.                                                   final TimeScales timeScales) {

  1141.             // set up nutation arguments
  1142.             final FundamentalNutationArguments arguments =
  1143.                     getNutationArguments(timeScales);

  1145.             // mean obliquity function
  1146.             final TimeScalarFunction epsilon = getMeanObliquityFunction(timeScales);

  1148.             // set up Poisson series
  1149.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  1150.             final PoissonSeriesParser luniSolarPsiParser =
  1151.                     new PoissonSeriesParser(14).
  1152.                         withFirstDelaunay(1).
  1153.                         withSinCos(0, 7, milliAS, 11, milliAS).
  1154.                         withSinCos(1, 8, milliAS, 12, milliAS);
  1155.             final PoissonSeries psiLuniSolarSeries;
  1156.             try (InputStream in = getStream(LUNI_SOLAR_SERIES)) {
  1157.                 psiLuniSolarSeries = luniSolarPsiParser.parse(in, LUNI_SOLAR_SERIES);
  1158.             } catch (IOException e) {
  1159.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1160.             }

  1162.             final PoissonSeriesParser planetaryPsiParser =
  1163.                     new PoissonSeriesParser(21).
  1164.                         withFirstDelaunay(2).
  1165.                         withFirstPlanetary(7).
  1166.                         withSinCos(0, 17, milliAS, 18, milliAS);
  1167.             final PoissonSeries psiPlanetarySeries;
  1168.             try (InputStream in = getStream(PLANETARY_SERIES)) {
  1169.                 psiPlanetarySeries = planetaryPsiParser.parse(in, PLANETARY_SERIES);
  1170.             } catch (IOException e) {
  1171.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1172.             }


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

  1191.             // ERA function
  1192.             final TimeScalarFunction era =
  1193.                     getEarthOrientationAngleFunction(ut1, timeScales.getTAI());

  1195.             return new TimeScalarFunction() {

  1197.                 /** {@inheritDoc} */
  1198.                 @Override
  1199.                 public double value(final AbsoluteDate date) {

  1201.                     // evaluate equation of origins
  1202.                     final BodiesElements elements = arguments.evaluateAll(date);
  1203.                     final double[] angles = psiGstSeries.value(elements);
  1204.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  1205.                     final double deltaPsi = angles[0] + angles[1] + ddPsi;
  1206.                     final double epsilonA = epsilon.value(date);

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

  1212.                 }

  1214.                 /** {@inheritDoc} */
  1215.                 @Override
  1216.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  1218.                     // evaluate equation of origins
  1219.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1220.                     final T[] angles = psiGstSeries.value(elements);
  1221.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  1222.                     final T deltaPsi = angles[0].add(angles[1]).add(ddPsi);
  1223.                     final T epsilonA = epsilon.value(date);

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

  1229.                 }

  1231.             };

  1233.         }

  1235.         /** {@inheritDoc} */
  1236.         @Override
  1237.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

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

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

  1265.                 final double microS = 1.0e-6;
  1266.                 final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  1267.                         withOptionalColumn(1).
  1268.                         withGamma(2).
  1269.                         withFirstDelaunay(3).
  1270.                         withSinCos(0, 10, microS, 11, microS);
  1271.                 final PoissonSeries ut1Series =
  1272.                         ut1Parser.parse(tidalUt1, TIDAL_CORRECTION_UT1_SERIES);

  1274.                 return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
  1275.             } catch (IOException e) {
  1276.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1277.             }

  1279.         }

  1281.         /** {@inheritDoc} */
  1282.         public LoveNumbers getLoveNumbers() {
  1283.             return loadLoveNumbers(LOVE_NUMBERS);
  1284.         }

  1286.         /** {@inheritDoc} */
  1287.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  1288.                                                                      final TimeScales timeScales) {

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

  1293.             // set up Poisson series
  1294.             final PoissonSeriesParser k20Parser =
  1295.                     new PoissonSeriesParser(18).
  1296.                         withOptionalColumn(1).
  1297.                         withDoodson(4, 2).
  1298.                         withFirstDelaunay(10);
  1299.             final PoissonSeriesParser k21Parser =
  1300.                     new PoissonSeriesParser(18).
  1301.                         withOptionalColumn(1).
  1302.                         withDoodson(4, 3).
  1303.                         withFirstDelaunay(10);
  1304.             final PoissonSeriesParser k22Parser =
  1305.                     new PoissonSeriesParser(16).
  1306.                         withOptionalColumn(1).
  1307.                         withDoodson(4, 2).
  1308.                         withFirstDelaunay(10);

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

  1337.                 return new TideFrequencyDependenceFunction(arguments,
  1338.                                                            c20Series,
  1339.                                                            c21Series, s21Series,
  1340.                                                            c22Series, s22Series);
  1341.             } catch (IOException e) {
  1342.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1343.             }

  1345.         }

  1347.         /** {@inheritDoc} */
  1348.         @Override
  1349.         public double getPermanentTide() {
  1350.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1351.         }

  1353.         /** {@inheritDoc} */
  1354.         @Override
  1355.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

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

  1382.             // polynomial extension from IERS 2003, section 7.1.4, equations 23a and 23b
  1383.             final double xp0    = 0.054   * Constants.ARC_SECONDS_TO_RADIANS;
  1384.             final double xp0Dot = 0.00083 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
  1385.             final double yp0    = 0.357   * Constants.ARC_SECONDS_TO_RADIANS;
  1386.             final double yp0Dot = 0.00395 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;

  1388.             // constants from IERS 2003, section 6.2
  1389.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1390.             final double ratio        =  0.00115;

  1392.             return new TimeVectorFunction() {

  1394.                 /** {@inheritDoc} */
  1395.                 @Override
  1396.                 public double[] value(final AbsoluteDate date) {

  1398.                     // we can't compute anything before the range covered by the annual pole file
  1399.                     if (date.compareTo(firstAnnualPoleDate) <= 0) {
  1400.                         return new double[] {
  1401.                             0.0, 0.0
  1402.                         };
  1403.                     }

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

  1427.                         // we are after the range covered by the annual pole file,
  1428.                         // we use the polynomial extension
  1429.                         final double t = date.durationFrom(
  1430.                                 eopHistory.getTimeScales().getJ2000Epoch());
  1431.                         meanPoleX = xp0 + t * xp0Dot;
  1432.                         meanPoleY = yp0 + t * yp0Dot;

  1434.                     }

  1436.                     // evaluate wobble variables
  1437.                     final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1438.                     final double m1 = correction.getXp() - meanPoleX;
  1439.                     final double m2 = meanPoleY - correction.getYp();

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

  1472.                 }

  1474.                 /** {@inheritDoc} */
  1475.                 @Override
  1476.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

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

  1512.                         // we are after the range covered by the annual pole file,
  1513.                         // we use the polynomial extension
  1514.                         final T t = date.durationFrom(
  1515.                                 eopHistory.getTimeScales().getJ2000Epoch());
  1516.                         meanPoleX = t.multiply(xp0Dot).add(xp0);
  1517.                         meanPoleY = t.multiply(yp0Dot).add(yp0);

  1519.                     }

  1521.                     // evaluate wobble variables
  1522.                     final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1523.                     final T m1 = correction.getXp().subtract(meanPoleX);
  1524.                     final T m2 = meanPoleY.subtract(correction.getYp());

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

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

  1557.                     return a;

  1559.                 }

  1561.             };

  1563.         }

  1565.         /** {@inheritDoc} */
  1566.         @Override
  1567.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {

  1569.             return new TimeVectorFunction() {

  1571.                 /** {@inheritDoc} */
  1572.                 @Override
  1573.                 public double[] value(final AbsoluteDate date) {
  1574.                     // there are no model for ocean pole tide prior to conventions 2010
  1575.                     return new double[] {
  1576.                         0.0, 0.0
  1577.                     };
  1578.                 }

  1580.                 /** {@inheritDoc} */
  1581.                 @Override
  1582.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1583.                     // there are no model for ocean pole tide prior to conventions 2010
  1584.                     return MathArrays.buildArray(date.getField(), 2);
  1585.                 }

  1587.             };
  1588.         }

  1590.         /** {@inheritDoc} */
  1591.         @Override
  1592.         public double[] getNominalTidalDisplacement() {
  1593.             return new double[] {
  1594.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  1595.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  1596.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  1597.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  1598.                 // H₀
  1599.                 -0.31460
  1600.             };
  1601.         }

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

  1610.         /** {@inheritDoc} */
  1611.         @Override
  1612.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  1613.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  1614.                                                                 18, 15, 16, 17, 18);
  1615.         }

  1617.     },

  1619.     /** Constant for IERS 2010 conventions. */
  1620.     IERS_2010 {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  1667.         /** {@inheritDoc} */
  1668.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
  1669.                                                                  final TimeScales timeScales) {
  1670.             try (InputStream in = getStream(NUTATION_ARGUMENTS)) {
  1671.                 return new FundamentalNutationArguments(this, timeScale,
  1672.                         in, NUTATION_ARGUMENTS, timeScales);
  1673.             } catch (IOException e) {
  1674.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1675.             }
  1676.         }

  1678.         /** {@inheritDoc} */
  1679.         @Override
  1680.         public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {

  1682.             // epsilon 0 value from chapter 5, page 56, other terms from equation 5.40 page 65
  1683.             final PolynomialNutation epsilonA =
  1684.                     new PolynomialNutation(84381.406        * Constants.ARC_SECONDS_TO_RADIANS,
  1685.                                              -46.836769     * Constants.ARC_SECONDS_TO_RADIANS,
  1686.                                               -0.0001831    * Constants.ARC_SECONDS_TO_RADIANS,
  1687.                                                0.00200340   * Constants.ARC_SECONDS_TO_RADIANS,
  1688.                                               -0.000000576  * Constants.ARC_SECONDS_TO_RADIANS,
  1689.                                               -0.0000000434 * Constants.ARC_SECONDS_TO_RADIANS);

  1691.             return new TimeScalarFunction() {

  1693.                 /** {@inheritDoc} */
  1694.                 @Override
  1695.                 public double value(final AbsoluteDate date) {
  1696.                     return epsilonA.value(evaluateTC(date, timeScales));
  1697.                 }

  1699.                 /** {@inheritDoc} */
  1700.                 @Override
  1701.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1702.                     return epsilonA.value(evaluateTC(date, timeScales));
  1703.                 }

  1705.             };

  1707.         }

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

  1713.             // set up nutation arguments
  1714.             final FundamentalNutationArguments arguments =
  1715.                     getNutationArguments(timeScales);

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

  1737.             // create a function evaluating the series
  1738.             return new TimeVectorFunction() {

  1740.                 /** {@inheritDoc} */
  1741.                 @Override
  1742.                 public double[] value(final AbsoluteDate date) {
  1743.                     return xys.value(arguments.evaluateAll(date));
  1744.                 }

  1746.                 /** {@inheritDoc} */
  1747.                 @Override
  1748.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1749.                     return xys.value(arguments.evaluateAll(date));
  1750.                 }

  1752.             };

  1754.         }

  1756.         /** {@inheritDoc} */
  1757.         public LoveNumbers getLoveNumbers() {
  1758.             return loadLoveNumbers(LOVE_NUMBERS);
  1759.         }

  1761.         /** {@inheritDoc} */
  1762.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
  1763.                                                                      final TimeScales timeScales) {

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

  1768.             // set up Poisson series
  1769.             final PoissonSeriesParser k20Parser =
  1770.                     new PoissonSeriesParser(18).
  1771.                         withOptionalColumn(1).
  1772.                         withDoodson(4, 2).
  1773.                         withFirstDelaunay(10);
  1774.             final PoissonSeriesParser k21Parser =
  1775.                     new PoissonSeriesParser(18).
  1776.                         withOptionalColumn(1).
  1777.                         withDoodson(4, 3).
  1778.                         withFirstDelaunay(10);
  1779.             final PoissonSeriesParser k22Parser =
  1780.                     new PoissonSeriesParser(16).
  1781.                         withOptionalColumn(1).
  1782.                         withDoodson(4, 2).
  1783.                         withFirstDelaunay(10);

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

  1812.                 return new TideFrequencyDependenceFunction(arguments,
  1813.                                                            c20Series,
  1814.                                                            c21Series, s21Series,
  1815.                                                            c22Series, s22Series);
  1816.             } catch (IOException e) {
  1817.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  1818.             }

  1820.         }

  1822.         /** {@inheritDoc} */
  1823.         @Override
  1824.         public double getPermanentTide() {
  1825.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1826.         }

  1828.         /** Compute pole wobble variables m₁ and m₂.
  1829.          * @param date current date
  1830.          * @param eopHistory EOP history
  1831.          * @return array containing m₁ and m₂
  1832.          */
  1833.         private double[] computePoleWobble(final AbsoluteDate date, final EOPHistory eopHistory) {

  1835.             // polynomial model from IERS 2010, table 7.7
  1836.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1837.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1838.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1839.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1840.             final AbsoluteDate changeDate =
  1841.                     new AbsoluteDate(2010, 1, 1, eopHistory.getTimeScales().getTT());

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

  1872.             // evaluate wobble variables
  1873.             final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1874.             final double m1 = correction.getXp() - meanPoleX;
  1875.             final double m2 = meanPoleY - correction.getYp();

  1877.             return new double[] {
  1878.                 m1, m2
  1879.             };

  1881.         }

  1883.         /** Compute pole wobble variables m₁ and m₂.
  1884.          * @param date current date
  1885.          * @param <T> type of the field elements
  1886.          * @param eopHistory EOP history
  1887.          * @return array containing m₁ and m₂
  1888.          */
  1889.         private <T extends RealFieldElement<T>> T[] computePoleWobble(final FieldAbsoluteDate<T> date, final EOPHistory eopHistory) {

  1891.             // polynomial model from IERS 2010, table 7.7
  1892.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1893.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1894.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1895.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1896.             final AbsoluteDate changeDate =
  1897.                     new AbsoluteDate(2010, 1, 1, eopHistory.getTimeScales().getTT());

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

  1928.             // evaluate wobble variables
  1929.             final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1930.             final T[] m = MathArrays.buildArray(date.getField(), 2);
  1931.             m[0] = correction.getXp().subtract(meanPoleX);
  1932.             m[1] = meanPoleY.subtract(correction.getYp());

  1934.             return m;

  1936.         }

  1938.         /** {@inheritDoc} */
  1939.         @Override
  1940.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  1942.             // constants from IERS 2010, section 6.4
  1943.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1944.             final double ratio        =  0.00115;

  1946.             return new TimeVectorFunction() {

  1948.                 /** {@inheritDoc} */
  1949.                 @Override
  1950.                 public double[] value(final AbsoluteDate date) {

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

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

  1968.                 }

  1970.                 /** {@inheritDoc} */
  1971.                 @Override
  1972.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

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

  1990.                     return a;

  1992.                 }

  1994.             };

  1996.         }

  1998.         /** {@inheritDoc} */
  1999.         @Override
  2000.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {

  2002.             return new TimeVectorFunction() {

  2004.                 /** {@inheritDoc} */
  2005.                 @Override
  2006.                 public double[] value(final AbsoluteDate date) {

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

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

  2020.                 }

  2022.                 /** {@inheritDoc} */
  2023.                 @Override
  2024.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

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

  2038.                     return v;

  2040.                 }

  2042.             };

  2044.         }

  2046.         /** {@inheritDoc} */
  2047.         @Override
  2048.         public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {

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

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

  2077.         }

  2079.         /** {@inheritDoc} */
  2080.         @Override
  2081.         public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {

  2083.             // set up nutation arguments
  2084.             final FundamentalNutationArguments arguments =
  2085.                     getNutationArguments(timeScales);

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

  2104.             return new TimeVectorFunction() {

  2106.                 /** {@inheritDoc} */
  2107.                 @Override
  2108.                 public double[] value(final AbsoluteDate date) {
  2109.                     final BodiesElements elements = arguments.evaluateAll(date);
  2110.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  2111.                     return new double[] {
  2112.                         psiEpsilon[0], psiEpsilon[1],
  2113.                         IAU1994ResolutionC7.value(elements, timeScales.getTAI())
  2114.                     };
  2115.                 }

  2117.                 /** {@inheritDoc} */
  2118.                 @Override
  2119.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  2120.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  2121.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  2122.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  2123.                     result[0] = psiEpsilon[0];
  2124.                     result[1] = psiEpsilon[1];
  2125.                     result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
  2126.                     return result;
  2127.                 }

  2129.             };

  2131.         }

  2133.         /** {@inheritDoc} */
  2134.         @Override
  2135.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
  2136.                                                   final TimeScales timeScales) {

  2138.             // Earth Rotation Angle
  2139.             final StellarAngleCapitaine era =
  2140.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

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

  2158.             // create a function evaluating the series
  2159.             return new TimeScalarFunction() {

  2161.                 /** {@inheritDoc} */
  2162.                 @Override
  2163.                 public double value(final AbsoluteDate date) {
  2164.                     return era.value(date) + minusEO.value(evaluateTC(date, timeScales));
  2165.                 }

  2167.                 /** {@inheritDoc} */
  2168.                 @Override
  2169.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  2170.                     return era.value(date).add(minusEO.value(evaluateTC(date, timeScales)));
  2171.                 }

  2173.             };

  2175.         }

  2177.         /** {@inheritDoc} */
  2178.         @Override
  2179.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
  2180.                                                       final TimeScales timeScales) {

  2182.             // Earth Rotation Angle
  2183.             final StellarAngleCapitaine era =
  2184.                     new StellarAngleCapitaine(ut1, timeScales.getTAI());

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

  2202.             // create a function evaluating the series
  2203.             return new TimeScalarFunction() {

  2205.                 /** {@inheritDoc} */
  2206.                 @Override
  2207.                 public double value(final AbsoluteDate date) {
  2208.                     return era.getRate() + minusEO.derivative(evaluateTC(date, timeScales));
  2209.                 }

  2211.                 /** {@inheritDoc} */
  2212.                 @Override
  2213.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  2214.                     return minusEO.derivative(evaluateTC(date, timeScales)).add(era.getRate());
  2215.                 }

  2217.             };

  2219.         }

  2221.         /** {@inheritDoc} */
  2222.         @Override
  2223.         public TimeScalarFunction getGASTFunction(final TimeScale ut1,
  2224.                                                   final EOPHistory eopHistory,
  2225.                                                   final TimeScales timeScales) {

  2227.             // set up nutation arguments
  2228.             final FundamentalNutationArguments arguments =
  2229.                     getNutationArguments(timeScales);

  2231.             // mean obliquity function
  2232.             final TimeScalarFunction epsilon = getMeanObliquityFunction(timeScales);

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

  2252.             // ERA function
  2253.             final TimeScalarFunction era =
  2254.                     getEarthOrientationAngleFunction(ut1, timeScales.getTAI());

  2256.             return new TimeScalarFunction() {

  2258.                 /** {@inheritDoc} */
  2259.                 @Override
  2260.                 public double value(final AbsoluteDate date) {

  2262.                     // evaluate equation of origins
  2263.                     final BodiesElements elements = arguments.evaluateAll(date);
  2264.                     final double[] angles = psiGstSeries.value(elements);
  2265.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  2266.                     final double deltaPsi = angles[0] + ddPsi;
  2267.                     final double epsilonA = epsilon.value(date);

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

  2273.                 }

  2275.                 /** {@inheritDoc} */
  2276.                 @Override
  2277.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  2279.                     // evaluate equation of origins
  2280.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  2281.                     final T[] angles = psiGstSeries.value(elements);
  2282.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  2283.                     final T deltaPsi = angles[0].add(ddPsi);
  2284.                     final T epsilonA = epsilon.value(date);

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

  2290.                 }

  2292.             };

  2294.         }

  2296.         /** {@inheritDoc} */
  2297.         @Override
  2298.         public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {

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

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

  2326.                 final double microS = 1.0e-6;
  2327.                 final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  2328.                         withOptionalColumn(1).
  2329.                         withGamma(2).
  2330.                         withFirstDelaunay(3).
  2331.                         withSinCos(0, 10, microS, 11, microS);
  2332.                 final PoissonSeries ut1Series =
  2333.                         ut1Parser.parse(ut1In, TIDAL_CORRECTION_UT1_SERIES);

  2335.                 return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
  2336.             } catch (IOException e) {
  2337.                 throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2338.             }

  2340.         }

  2342.         /** {@inheritDoc} */
  2343.         @Override
  2344.         public double[] getNominalTidalDisplacement() {
  2345.             return new double[] {
  2346.                 // h⁽⁰⁾,                                  h⁽²⁾,   h₃,     hI diurnal, hI semi-diurnal,
  2347.                 0.6078,                                  -0.0006, 0.292, -0.0025,    -0.0022,
  2348.                 // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾,   l₃,     lI diurnal, lI semi-diurnal
  2349.                 0.0847,  0.0012,       0.0024,            0.0002, 0.015, -0.0007,    -0.0007,
  2350.                 // H₀
  2351.                 -0.31460
  2352.             };
  2353.         }

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

  2362.         /** {@inheritDoc} */
  2363.         @Override
  2364.         public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
  2365.             return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
  2366.                                                                 18, 15, 16, 17, 18);
  2367.         }

  2369.     };

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

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

  2377.     /** Get the reference epoch for fundamental nutation arguments.
  2378.      *
  2379.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2380.      *
  2381.      * @return reference epoch for fundamental nutation arguments
  2382.      * @since 6.1
  2383.      * @see #getNutationReferenceEpoch(TimeScales)
  2384.      */
  2385.     @DefaultDataContext
  2386.     public AbsoluteDate getNutationReferenceEpoch() {
  2387.         return getNutationReferenceEpoch(DataContext.getDefault().getTimeScales());
  2388.     }

  2390.     /**
  2391.      * Get the reference epoch for fundamental nutation arguments.
  2392.      *
  2393.      * @param timeScales to use for the reference epoch.
  2394.      * @return reference epoch for fundamental nutation arguments
  2395.      * @since 10.1
  2396.      */
  2397.     public AbsoluteDate getNutationReferenceEpoch(final TimeScales timeScales) {
  2398.         // IERS 1996, IERS 2003 and IERS 2010 use the same J2000.0 reference date
  2399.         return timeScales.getJ2000Epoch();
  2400.     }

  2402.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2403.      *
  2404.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2405.      *
  2406.      * @param date current date
  2407.      * @return date offset in Julian centuries
  2408.      * @since 6.1
  2409.      * @see #evaluateTC(AbsoluteDate, TimeScales)
  2410.      */
  2411.     @DefaultDataContext
  2412.     public double evaluateTC(final AbsoluteDate date) {
  2413.         return evaluateTC(date, DataContext.getDefault().getTimeScales());
  2414.     }

  2416.     /**
  2417.      * Evaluate the date offset between the current date and the {@link
  2418.      * #getNutationReferenceEpoch() reference date}.
  2419.      *
  2420.      * @param date       current date
  2421.      * @param timeScales used in the evaluation.
  2422.      * @return date offset in Julian centuries
  2423.      * @since 10.1
  2424.      */
  2425.     public double evaluateTC(final AbsoluteDate date, final TimeScales timeScales) {
  2426.         return date.durationFrom(getNutationReferenceEpoch(timeScales)) /
  2427.                 Constants.JULIAN_CENTURY;
  2428.     }

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

  2445.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2446.      * @param <T> type of the field elements
  2447.      * @param date current date
  2448.      * @param timeScales used in the evaluation.
  2449.      * @return date offset in Julian centuries
  2450.      * @since 10.1
  2451.      */
  2452.     public <T extends RealFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date,
  2453.                                                         final TimeScales timeScales) {
  2454.         return date.durationFrom(getNutationReferenceEpoch(timeScales))
  2455.                 .divide(Constants.JULIAN_CENTURY);
  2456.     }

  2458.     /**
  2459.      * Get the fundamental nutation arguments. Does not compute GMST based values: gamma,
  2460.      * gammaDot.
  2461.      *
  2462.      * @param timeScales other time scales used in the computation including TAI and TT.
  2463.      * @return fundamental nutation arguments
  2464.      * @see #getNutationArguments(TimeScale, TimeScales)
  2465.      * @since 10.1
  2466.      */
  2467.     protected FundamentalNutationArguments getNutationArguments(
  2468.             final TimeScales timeScales) {

  2470.         return getNutationArguments(null, timeScales);
  2471.     }

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

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

  2502.     /** Get the function computing mean obliquity of the ecliptic.
  2503.      *
  2504.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2505.      *
  2506.      * @return function computing mean obliquity of the ecliptic
  2507.      * @since 6.1
  2508.      * @see #getMeanObliquityFunction(TimeScales)
  2509.      */
  2510.     @DefaultDataContext
  2511.     public TimeScalarFunction getMeanObliquityFunction() {
  2512.         return getMeanObliquityFunction(DataContext.getDefault().getTimeScales());
  2513.     }

  2515.     /**
  2516.      * Get the function computing mean obliquity of the ecliptic.
  2517.      *
  2518.      * @param timeScales used in computing the function.
  2519.      * @return function computing mean obliquity of the ecliptic
  2520.      * @since 10.1
  2521.      */
  2522.     public abstract TimeScalarFunction getMeanObliquityFunction(TimeScales timeScales);

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

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

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

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

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

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

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

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

  2661.     /** Get the function computing Greenwich mean sidereal time, in radians.
  2662.      *
  2663.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2664.      *
  2665.      * @param ut1 UT1 time scale
  2666.      * @return function computing Greenwich mean sidereal time
  2667.      * @since 6.1
  2668.      * @see #getGMSTFunction(TimeScale, TimeScales)
  2669.      */
  2670.     @DefaultDataContext
  2671.     public TimeScalarFunction getGMSTFunction(final TimeScale ut1) {
  2672.         return getGMSTFunction(ut1, DataContext.getDefault().getTimeScales());
  2673.     }

  2675.     /**
  2676.      * Get the function computing Greenwich mean sidereal time, in radians.
  2677.      *
  2678.      * @param ut1 UT1 time scale
  2679.      * @param timeScales other time scales used in the computation including TAI and TT.
  2680.      * @return function computing Greenwich mean sidereal time
  2681.      * @since 10.1
  2682.      */
  2683.     public abstract TimeScalarFunction getGMSTFunction(TimeScale ut1,
  2684.                                                        TimeScales timeScales);

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

  2701.     /**
  2702.      * Get the function computing Greenwich mean sidereal time rate, in radians per
  2703.      * second.
  2704.      *
  2705.      * @param ut1 UT1 time scale
  2706.      * @param timeScales other time scales used in the computation including TAI and TT.
  2707.      * @return function computing Greenwich mean sidereal time rate
  2708.      * @since 10.1
  2709.      */
  2710.     public abstract TimeScalarFunction getGMSTRateFunction(TimeScale ut1,
  2711.                                                            TimeScales timeScales);

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

  2735.     /**
  2736.      * Get the function computing Greenwich apparent sidereal time, in radians.
  2737.      *
  2738.      * @param ut1        UT1 time scale
  2739.      * @param eopHistory EOP history. If {@code null} then no nutation correction is
  2740.      *                   applied for EOP.
  2741.      * @param timeScales        TAI time scale.
  2742.      * @return function computing Greenwich apparent sidereal time
  2743.      * @since 10.1
  2744.      */
  2745.     public abstract TimeScalarFunction getGASTFunction(TimeScale ut1,
  2746.                                                        EOPHistory eopHistory,
  2747.                                                        TimeScales timeScales);

  2749.     /** Get the function computing tidal corrections for Earth Orientation Parameters.
  2750.      *
  2751.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  2752.      *
  2753.      * @return function computing tidal corrections for Earth Orientation Parameters,
  2754.      * for xp, yp, ut1 and lod respectively
  2755.      * @since 6.1
  2756.      * @see #getEOPTidalCorrection(TimeScales)
  2757.      */
  2758.     @DefaultDataContext
  2759.     public TimeVectorFunction getEOPTidalCorrection() {
  2760.         return getEOPTidalCorrection(DataContext.getDefault().getTimeScales());
  2761.     }

  2763.     /**
  2764.      * Get the function computing tidal corrections for Earth Orientation Parameters.
  2765.      *
  2766.      * @param timeScales used in the computation. The TT and TAI scales are used.
  2767.      * @return function computing tidal corrections for Earth Orientation Parameters, for
  2768.      * xp, yp, ut1 and lod respectively
  2769.      * @since 10.1
  2770.      */
  2771.     public abstract TimeVectorFunction getEOPTidalCorrection(TimeScales timeScales);

  2773.     /** Get the Love numbers.
  2774.      * @return Love numbers
  2775.           * @since 6.1
  2776.      */
  2777.     public abstract LoveNumbers getLoveNumbers();

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

  2794.     /**
  2795.      * Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂,
  2796.      * ΔS₂₂).
  2797.      *
  2798.      * @param ut1 UT1 time scale
  2799.      * @param timeScales other time scales used in the computation including TAI and TT.
  2800.      * @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂,
  2801.      * ΔS₂₂).
  2802.      * @since 10.1
  2803.      */
  2804.     public abstract TimeVectorFunction getTideFrequencyDependenceFunction(
  2805.             TimeScale ut1,
  2806.             TimeScales timeScales);

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

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

  2820.     /** Get the function computing ocean pole tide (ΔC₂₁, ΔS₂₁).
  2821.      * @param eopHistory EOP history
  2822.      * @return model for ocean pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2823.           * @since 6.1
  2824.      */
  2825.     public abstract TimeVectorFunction getOceanPoleTide(EOPHistory eopHistory);

  2827.     /** Get the nominal values of the displacement numbers.
  2828.      * @return an array containing h⁽⁰⁾, h⁽²⁾, h₃, hI diurnal, hI semi-diurnal,
  2829.      * l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾, l₃, lI diurnal, lI semi-diurnal,
  2830.      * H₀ permanent deformation amplitude
  2831.      * @since 9.1
  2832.      */
  2833.     public abstract double[] getNominalTidalDisplacement();

  2835.     /** Get the correction function for tidal displacement for diurnal tides.
  2836.      * <ul>
  2837.      *  <li>f[0]: radial correction, longitude cosine part</li>
  2838.      *  <li>f[1]: radial correction, longitude sine part</li>
  2839.      *  <li>f[2]: North correction, longitude cosine part</li>
  2840.      *  <li>f[3]: North correction, longitude sine part</li>
  2841.      *  <li>f[4]: East correction, longitude cosine part</li>
  2842.      *  <li>f[5]: East correction, longitude sine part</li>
  2843.      * </ul>
  2844.      * @return correction function for tidal displacement
  2845.      * @since 9.1
  2846.      */
  2847.     public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal();

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

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

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

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

  2928.             return PoissonSeries.compile(drCos, drSin,
  2929.                                          dnCos, dnSin,
  2930.                                          deCos, deSin);
  2931.         } catch (IOException e) {
  2932.             throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2933.         }

  2935.     }

  2937.     /** Get the correction function for tidal displacement for zonal tides.
  2938.      * <ul>
  2939.      *  <li>f[0]: radial correction</li>
  2940.      *  <li>f[1]: North correction</li>
  2941.      * </ul>
  2942.      * @return correction function for tidal displacement
  2943.      * @since 9.1
  2944.      */
  2945.     public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionZonal();

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

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

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

  2989.             return PoissonSeries.compile(dr, dn);
  2990.         } catch (IOException e) {
  2991.             throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
  2992.         }

  2994.     }

  2996.     /** Interface for functions converting nutation corrections between
  2997.      * δΔψ/δΔε to δX/δY.
  2998.      * <ul>
  2999.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  3000.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  3001.      * </ul>
  3002.      * @since 6.1
  3003.      */
  3004.     public interface NutationCorrectionConverter {

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

  3014.         /** Convert nutation corrections.
  3015.          * @param date current date
  3016.          * @param dX δX part of the nutation correction
  3017.          * @param dY δY part of the nutation correction
  3018.          * @return array containing δΔψ and δΔε
  3019.          */
  3020.         double[] toEquinox(AbsoluteDate date, double dX, double dY);

  3022.     }

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

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

  3055.         // get models parameters
  3056.         final TimeVectorFunction precessionFunction = getPrecessionFunction(timeScales);
  3057.         final TimeScalarFunction epsilonAFunction = getMeanObliquityFunction(timeScales);
  3058.         final AbsoluteDate date0 = getNutationReferenceEpoch(timeScales);
  3059.         final double cosE0 = FastMath.cos(epsilonAFunction.value(date0));

  3061.         return new NutationCorrectionConverter() {

  3063.             /** {@inheritDoc} */
  3064.             @Override
  3065.             public double[] toNonRotating(final AbsoluteDate date,
  3066.                                           final double ddPsi, final double ddEpsilon) {
  3067.                 // compute precession angles psiA, omegaA and chiA
  3068.                 final double[] angles = precessionFunction.value(date);

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

  3074.                 // convert nutation corrections (equation 23/IERS-2003 or 5.25/IERS-2010)
  3075.                 return new double[] {
  3076.                     sinEA * ddPsi + c * ddEpsilon,
  3077.                     ddEpsilon - c * sinEA * ddPsi
  3078.                 };

  3080.             }

  3082.             /** {@inheritDoc} */
  3083.             @Override
  3084.             public double[] toEquinox(final AbsoluteDate date,
  3085.                                       final double dX, final double dY) {
  3086.                 // compute precession angles psiA, omegaA and chiA
  3087.                 final double[] angles   = precessionFunction.value(date);

  3089.                 // conversion coefficients
  3090.                 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
  3091.                 final double c     = angles[0] * cosE0 - angles[2];
  3092.                 final double opC2  = 1 + c * c;

  3094.                 // convert nutation corrections (inverse of equation 23/IERS-2003 or 5.25/IERS-2010)
  3095.                 return new double[] {
  3096.                     (dX - c * dY) / (sinEA * opC2),
  3097.                     (dY + c * dX) / opC2
  3098.                 };
  3099.             }

  3101.         };

  3103.     }

  3105.     /** Load the Love numbers.
  3106.      * @param nameLove name of the Love number resource
  3107.      * @return Love numbers
  3108.      */
  3109.     protected LoveNumbers loadLoveNumbers(final String nameLove) {
  3110.         try {

  3112.             // allocate the three triangular arrays for real, imaginary and time-dependent numbers
  3113.             final double[][] real      = new double[4][];
  3114.             final double[][] imaginary = new double[4][];
  3115.             final double[][] plus      = new double[4][];
  3116.             for (int i = 0; i < real.length; ++i) {
  3117.                 real[i]      = new double[i + 1];
  3118.                 imaginary[i] = new double[i + 1];
  3119.                 plus[i]      = new double[i + 1];
  3120.             }

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

  3124.                 if (stream == null) {
  3125.                     // this should never happen with files embedded within Orekit
  3126.                     throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, nameLove);
  3127.                 }

  3129.                 int lineNumber = 1;
  3130.                 String line = null;
  3131.                 // setup the reader
  3132.                 try (BufferedReader reader = new BufferedReader(new InputStreamReader(stream, StandardCharsets.UTF_8))) {

  3134.                     line = reader.readLine();

  3136.                     // look for the Love numbers
  3137.                     while (line != null) {

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

  3154.                             }
  3155.                             real[n][m]      = Double.parseDouble(fields[2]);
  3156.                             imaginary[n][m] = Double.parseDouble(fields[3]);
  3157.                             plus[n][m]      = Double.parseDouble(fields[4]);
  3158.                         }

  3160.                         // next line
  3161.                         lineNumber++;
  3162.                         line = reader.readLine();

  3164.                     }
  3165.                 } catch (NumberFormatException nfe) {
  3166.                     // this should never happen with files embedded within Orekit
  3167.                     throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  3168.                                               lineNumber, nameLove, line);
  3169.                 }
  3170.             }

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

  3174.         } catch (IOException ioe) {
  3175.             // this should never happen with files embedded within Orekit
  3176.             throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, nameLove);
  3177.         }
  3178.     }

  3180.     /** Get a data stream.
  3181.      * @param name file name of the resource stream
  3182.      * @return stream
  3183.      */
  3184.     private static InputStream getStream(final String name) {
  3185.         return IERSConventions.class.getResourceAsStream(name);
  3186.     }

  3188.     /** Correction to equation of equinoxes.
  3189.      * <p>IAU 1994 resolution C7 added two terms to the equation of equinoxes
  3190.      * taking effect since 1997-02-27 for continuity
  3191.      * </p>
  3192.      */
  3193.     private static class IAU1994ResolutionC7 {

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

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


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

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

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

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

  3224.             } else {
  3225.                 return 0.0;
  3226.             }
  3227.         }

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

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

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

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

  3253.             } else {
  3254.                 return arguments.getDate().getField().getZero();
  3255.             }
  3256.         }

  3258.     };

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

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

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

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

  3288.         /** UT1 time scale. */
  3289.         private final TimeScale ut1;

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

  3294.         /** Simple constructor.
  3295.          * @param ut1 UT1 time scale
  3296.          * @param tai TAI time scale
  3297.          */
  3298.         StellarAngleCapitaine(final TimeScale ut1, final TimeScale tai) {
  3299.             this.ut1 = ut1;
  3300.             referenceDate = new AbsoluteDate(
  3301.                     DateComponents.J2000_EPOCH,
  3302.                     TimeComponents.H12,
  3303.                     tai);
  3304.         }

  3306.         /** Get the rotation rate.
  3307.          * @return rotation rate
  3308.          */
  3309.         public double getRate() {
  3310.             return ERA_1A + ERA_1B;
  3311.         }

  3313.         /** {@inheritDoc} */
  3314.         @Override
  3315.         public double value(final AbsoluteDate date) {

  3317.             // split the date offset as a full number of days plus a smaller part
  3318.             final int secondsInDay = 86400;
  3319.             final double dt  = date.durationFrom(referenceDate);
  3320.             final long days  = ((long) dt) / secondsInDay;
  3321.             final double dtA = secondsInDay * days;
  3322.             final double dtB = (dt - dtA) + ut1.offsetFromTAI(date);

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

  3326.         }

  3328.         /** {@inheritDoc} */
  3329.         @Override
  3330.         public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  3332.             // split the date offset as a full number of days plus a smaller part
  3333.             final int secondsInDay = 86400;
  3334.             final T dt  = date.durationFrom(referenceDate);
  3335.             final long days  = ((long) dt.getReal()) / secondsInDay;
  3336.             final double dtA = secondsInDay * days;
  3337.             final T dtB = dt.subtract(dtA).add(ut1.offsetFromTAI(date.toAbsoluteDate()));

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

  3341.         }

  3343.     }

  3345.     /** Mean pole. */
  3346.     private static class MeanPole implements TimeStamped, Serializable {

  3348.         /** Serializable UID. */
  3349.         private static final long serialVersionUID = 20131028l;

  3351.         /** Date. */
  3352.         private final AbsoluteDate date;

  3354.         /** X coordinate. */
  3355.         private double x;

  3357.         /** Y coordinate. */
  3358.         private double y;

  3360.         /** Simple constructor.
  3361.          * @param date date
  3362.          * @param x x coordinate
  3363.          * @param y y coordinate
  3364.          */
  3365.         MeanPole(final AbsoluteDate date, final double x, final double y) {
  3366.             this.date = date;
  3367.             this.x    = x;
  3368.             this.y    = y;
  3369.         }

  3371.         /** {@inheritDoc} */
  3372.         @Override
  3373.         public AbsoluteDate getDate() {
  3374.             return date;
  3375.         }

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

  3384.         /** Get y coordinate.
  3385.          * @return y coordinate
  3386.          */
  3387.         public double getY() {
  3388.             return y;
  3389.         }

  3391.     }

  3393.     /** Local class for precession function. */
  3394.     private class PrecessionFunction implements TimeVectorFunction {

  3396.         /** Polynomial nutation for psiA. */
  3397.         private final PolynomialNutation psiA;

  3399.         /** Polynomial nutation for omegaA. */
  3400.         private final PolynomialNutation omegaA;

  3402.         /** Polynomial nutation for chiA. */
  3403.         private final PolynomialNutation chiA;

  3405.         /** Time scales to use. */
  3406.         private final TimeScales timeScales;

  3408.         /** Simple constructor.
  3409.          * @param psiA polynomial nutation for psiA
  3410.          * @param omegaA polynomial nutation for omegaA
  3411.          * @param chiA polynomial nutation for chiA
  3412.          * @param timeScales used in the computation.
  3413.          */
  3414.         PrecessionFunction(final PolynomialNutation psiA,
  3415.                            final PolynomialNutation omegaA,
  3416.                            final PolynomialNutation chiA,
  3417.                            final TimeScales timeScales) {
  3418.             this.psiA   = psiA;
  3419.             this.omegaA = omegaA;
  3420.             this.chiA   = chiA;
  3421.             this.timeScales = timeScales;
  3422.         }


  3425.         /** {@inheritDoc} */
  3426.         @Override
  3427.         public double[] value(final AbsoluteDate date) {
  3428.             final double tc = evaluateTC(date, timeScales);
  3429.             return new double[] {
  3430.                 psiA.value(tc), omegaA.value(tc), chiA.value(tc)
  3431.             };
  3432.         }

  3434.         /** {@inheritDoc} */
  3435.         @Override
  3436.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  3437.             final T[] a = MathArrays.buildArray(date.getField(), 3);
  3438.             final T tc = evaluateTC(date, timeScales);
  3439.             a[0] = psiA.value(tc);
  3440.             a[1] = omegaA.value(tc);
  3441.             a[2] = chiA.value(tc);
  3442.             return a;
  3443.         }

  3445.     }

  3447.     /** Local class for tides frequency function. */
  3448.     private static class TideFrequencyDependenceFunction implements TimeVectorFunction {

  3450.         /** Nutation arguments. */
  3451.         private final FundamentalNutationArguments arguments;

  3453.         /** Correction series. */
  3454.         private final PoissonSeries.CompiledSeries kSeries;

  3456.         /** Simple constructor.
  3457.          * @param arguments nutation arguments
  3458.          * @param c20Series correction series for the C20 term
  3459.          * @param c21Series correction series for the C21 term
  3460.          * @param s21Series correction series for the S21 term
  3461.          * @param c22Series correction series for the C22 term
  3462.          * @param s22Series correction series for the S22 term
  3463.          */
  3464.         TideFrequencyDependenceFunction(final FundamentalNutationArguments arguments,
  3465.                                         final PoissonSeries c20Series,
  3466.                                         final PoissonSeries c21Series,
  3467.                                         final PoissonSeries s21Series,
  3468.                                         final PoissonSeries c22Series,
  3469.                                         final PoissonSeries s22Series) {
  3470.             this.arguments = arguments;
  3471.             this.kSeries   = PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);
  3472.         }


  3475.         /** {@inheritDoc} */
  3476.         @Override
  3477.         public double[] value(final AbsoluteDate date) {
  3478.             return kSeries.value(arguments.evaluateAll(date));
  3479.         }

  3481.         /** {@inheritDoc} */
  3482.         @Override
  3483.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  3484.             return kSeries.value(arguments.evaluateAll(date));
  3485.         }

  3487.     }

  3489.     /** Local class for EOP tidal corrections. */
  3490.     private static class EOPTidalCorrection implements TimeVectorFunction {

  3492.         /** Nutation arguments. */
  3493.         private final FundamentalNutationArguments arguments;

  3495.         /** Correction series. */
  3496.         private final PoissonSeries.CompiledSeries correctionSeries;

  3498.         /** Simple constructor.
  3499.          * @param arguments nutation arguments
  3500.          * @param xSeries correction series for the x coordinate
  3501.          * @param ySeries correction series for the y coordinate
  3502.          * @param ut1Series correction series for the UT1
  3503.          */
  3504.         EOPTidalCorrection(final FundamentalNutationArguments arguments,
  3505.                            final PoissonSeries xSeries,
  3506.                            final PoissonSeries ySeries,
  3507.                            final PoissonSeries ut1Series) {
  3508.             this.arguments        = arguments;
  3509.             this.correctionSeries = PoissonSeries.compile(xSeries, ySeries, ut1Series);
  3510.         }

  3512.         /** {@inheritDoc} */
  3513.         @Override
  3514.         public double[] value(final AbsoluteDate date) {
  3515.             final BodiesElements elements = arguments.evaluateAll(date);
  3516.             final double[] correction    = correctionSeries.value(elements);
  3517.             final double[] correctionDot = correctionSeries.derivative(elements);
  3518.             return new double[] {
  3519.                 correction[0],
  3520.                 correction[1],
  3521.                 correction[2],
  3522.                 correctionDot[2] * (-Constants.JULIAN_DAY)
  3523.             };
  3524.         }

  3526.         /** {@inheritDoc} */
  3527.         @Override
  3528.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

  3530.             final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  3531.             final T[] correction    = correctionSeries.value(elements);
  3532.             final T[] correctionDot = correctionSeries.derivative(elements);
  3533.             final T[] a = MathArrays.buildArray(date.getField(), 4);
  3534.             a[0] = correction[0];
  3535.             a[1] = correction[1];
  3536.             a[2] = correction[2];
  3537.             a[3] = correctionDot[2].multiply(-Constants.JULIAN_DAY);
  3538.             return a;
  3539.         }

  3541.     }

  3543. }
Created with JaCoCo 0.8.6.202009150832