IERSConventions.java

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

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

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


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

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

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

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

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

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

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

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

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

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

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

  97.         /** {@inheritDoc} */
  98.         @Override
  99.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
  100.             throws OrekitException {
  101.             return new FundamentalNutationArguments(this, timeScale,
  102.                                                     getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
  103.         }

  105.         /** {@inheritDoc} */
  106.         @Override
  107.         public TimeScalarFunction getMeanObliquityFunction() throws OrekitException {

  109.             // value from chapter 5, page 22
  110.             final PolynomialNutation epsilonA =
  111.                     new PolynomialNutation(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  112.                                              -46.8150   * Constants.ARC_SECONDS_TO_RADIANS,
  113.                                               -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  114.                                                0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  116.             return new TimeScalarFunction() {

  118.                 /** {@inheritDoc} */
  119.                 @Override
  120.                 public double value(final AbsoluteDate date) {
  121.                     return epsilonA.value(evaluateTC(date));
  122.                 }

  124.                 /** {@inheritDoc} */
  125.                 @Override
  126.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  127.                     return epsilonA.value(evaluateTC(date));
  128.                 }

  130.             };

  132.         }

  134.         /** {@inheritDoc} */
  135.         @Override
  136.         public TimeVectorFunction getXYSpXY2Function()
  137.             throws OrekitException {

  139.             // set up nutation arguments
  140.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  142.             // X = 2004.3109″t - 0.42665″t² - 0.198656″t³ + 0.0000140″t⁴
  143.             //     + 0.00006″t² cos Ω + sin ε0 { Σ [(Ai + Ai' t) sin(ARGUMENT) + Ai'' t cos(ARGUMENT)]}
  144.             //     + 0.00204″t² sin Ω + 0.00016″t² sin 2(F - D + Ω),
  145.             final PolynomialNutation xPolynomial =
  146.                     new PolynomialNutation(0,
  147.                                            2004.3109 * Constants.ARC_SECONDS_TO_RADIANS,
  148.                                            -0.42665  * Constants.ARC_SECONDS_TO_RADIANS,
  149.                                            -0.198656 * Constants.ARC_SECONDS_TO_RADIANS,
  150.                                            0.0000140 * Constants.ARC_SECONDS_TO_RADIANS);

  152.             final double fXCosOm    = 0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
  153.             final double fXSinOm    = 0.00204 * Constants.ARC_SECONDS_TO_RADIANS;
  154.             final double fXSin2FDOm = 0.00016 * Constants.ARC_SECONDS_TO_RADIANS;
  155.             final double sinEps0   = FastMath.sin(getMeanObliquityFunction().value(getNutationReferenceEpoch()));

  157.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  158.             final PoissonSeriesParser baseParser =
  159.                     new PoissonSeriesParser(12).withFirstDelaunay(1);

  161.             final PoissonSeriesParser xParser =
  162.                     baseParser.
  163.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  164.                     withSinCos(1, 8, deciMilliAS,  9, deciMilliAS);
  165.             final PoissonSeries xSum = xParser.parse(getStream(X_Y_SERIES), X_Y_SERIES);

  167.             // Y = -0.00013″ - 22.40992″t² + 0.001836″t³ + 0.0011130″t⁴
  168.             //     + Σ [(Bi + Bi' t) cos(ARGUMENT) + Bi'' t sin(ARGUMENT)]
  169.             //    - 0.00231″t² cos Ω − 0.00014″t² cos 2(F - D + Ω)
  170.             final PolynomialNutation yPolynomial =
  171.                     new PolynomialNutation(-0.00013  * Constants.ARC_SECONDS_TO_RADIANS,
  172.                                            0.0,
  173.                                            -22.40992 * Constants.ARC_SECONDS_TO_RADIANS,
  174.                                            0.001836  * Constants.ARC_SECONDS_TO_RADIANS,
  175.                                            0.0011130 * Constants.ARC_SECONDS_TO_RADIANS);

  177.             final double fYCosOm    = -0.00231 * Constants.ARC_SECONDS_TO_RADIANS;
  178.             final double fYCos2FDOm = -0.00014 * Constants.ARC_SECONDS_TO_RADIANS;

  180.             final PoissonSeriesParser yParser =
  181.                     baseParser.
  182.                     withSinCos(0, -1, deciMilliAS, 10, deciMilliAS).
  183.                     withSinCos(1, 12, deciMilliAS, 11, deciMilliAS);
  184.             final PoissonSeries ySum = yParser.parse(getStream(X_Y_SERIES), X_Y_SERIES);

  186.             final PoissonSeries.CompiledSeries xySum =
  187.                     PoissonSeries.compile(xSum, ySum);

  189.             // s = -XY/2 + 0.00385″t - 0.07259″t³ - 0.00264″ sin Ω - 0.00006″ sin 2Ω
  190.             //     + 0.00074″t² sin Ω + 0.00006″t² sin 2(F - D + Ω)
  191.             final double fST          =  0.00385 * Constants.ARC_SECONDS_TO_RADIANS;
  192.             final double fST3         = -0.07259 * Constants.ARC_SECONDS_TO_RADIANS;
  193.             final double fSSinOm      = -0.00264 * Constants.ARC_SECONDS_TO_RADIANS;
  194.             final double fSSin2Om     = -0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
  195.             final double fST2SinOm    =  0.00074 * Constants.ARC_SECONDS_TO_RADIANS;
  196.             final double fST2Sin2FDOm =  0.00006 * Constants.ARC_SECONDS_TO_RADIANS;

  198.             return new TimeVectorFunction() {

  200.                 /** {@inheritDoc} */
  201.                 @Override
  202.                 public double[] value(final AbsoluteDate date) {

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

  207.                     final double omega     = elements.getOmega();
  208.                     final double f         = elements.getF();
  209.                     final double d         = elements.getD();
  210.                     final double t         = elements.getTC();

  212.                     final double cosOmega  = FastMath.cos(omega);
  213.                     final double sinOmega  = FastMath.sin(omega);
  214.                     final double sin2Omega = FastMath.sin(2 * omega);
  215.                     final double cos2FDOm  = FastMath.cos(2 * (f - d + omega));
  216.                     final double sin2FDOm  = FastMath.sin(2 * (f - d + omega));

  218.                     final double x = xPolynomial.value(t) + sinEps0 * xy[0] +
  219.                             t * t * (fXCosOm * cosOmega + fXSinOm * sinOmega + fXSin2FDOm * cos2FDOm);
  220.                     final double y = yPolynomial.value(t) + xy[1] +
  221.                             t * t * (fYCosOm * cosOmega + fYCos2FDOm * cos2FDOm);
  222.                     final double sPxy2 = fSSinOm * sinOmega + fSSin2Om * sin2Omega +
  223.                             t * (fST + t * (fST2SinOm * sinOmega + fST2Sin2FDOm * sin2FDOm + t * fST3));

  225.                     return new double[] {
  226.                         x, y, sPxy2
  227.                     };

  229.                 }

  231.                 /** {@inheritDoc} */
  232.                 @Override
  233.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

  238.                     final T omega     = elements.getOmega();
  239.                     final T f         = elements.getF();
  240.                     final T d         = elements.getD();
  241.                     final T t         = elements.getTC();
  242.                     final T t2        = t.multiply(t);

  244.                     final T cosOmega  = omega.cos();
  245.                     final T sinOmega  = omega.sin();
  246.                     final T sin2Omega = omega.multiply(2).sin();
  247.                     final T fMDPO2 = f.subtract(d).add(omega).multiply(2);
  248.                     final T cos2FDOm  = fMDPO2.cos();
  249.                     final T sin2FDOm  = fMDPO2.sin();

  251.                     final T x = xPolynomial.value(t).
  252.                                 add(xy[0].multiply(sinEps0)).
  253.                                 add(t2.multiply(cosOmega.multiply(fXCosOm).add(sinOmega.multiply(fXSinOm)).add(cos2FDOm.multiply(fXSin2FDOm))));
  254.                     final T y = yPolynomial.value(t).
  255.                                 add(xy[1]).
  256.                                 add(t2.multiply(cosOmega.multiply(fYCosOm).add(cos2FDOm.multiply(fYCos2FDOm))));
  257.                     final T sPxy2 = sinOmega.multiply(fSSinOm).
  258.                                     add(sin2Omega.multiply(fSSin2Om)).
  259.                                     add(t.multiply(fST3).add(sinOmega.multiply(fST2SinOm)).add(sin2FDOm.multiply(fST2Sin2FDOm)).multiply(t).add(fST).multiply(t));

  261.                     final T[] a = MathArrays.buildArray(date.getField(), 3);
  262.                     a[0] = x;
  263.                     a[1] = y;
  264.                     a[2] = sPxy2;
  265.                     return a;

  267.                 }

  269.             };

  271.         }

  273.         /** {@inheritDoc} */
  274.         @Override
  275.         public TimeVectorFunction getPrecessionFunction() throws OrekitException {

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

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

  300.         }

  302.         /** {@inheritDoc} */
  303.         @Override
  304.         public TimeVectorFunction getNutationFunction()
  305.             throws OrekitException {

  307.             // set up nutation arguments
  308.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  310.             // set up Poisson series
  311.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  312.             final PoissonSeriesParser baseParser =
  313.                     new PoissonSeriesParser(10).withFirstDelaunay(1);

  315.             final PoissonSeriesParser psiParser =
  316.                     baseParser.
  317.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  318.                     withSinCos(1, 8, deciMilliAS, -1, deciMilliAS);
  319.             final PoissonSeries psiSeries = psiParser.parse(getStream(PSI_EPSILON_SERIES), PSI_EPSILON_SERIES);

  321.             final PoissonSeriesParser epsilonParser =
  322.                     baseParser.
  323.                     withSinCos(0, -1, deciMilliAS, 9, deciMilliAS).
  324.                     withSinCos(1, -1, deciMilliAS, 10, deciMilliAS);
  325.             final PoissonSeries epsilonSeries = epsilonParser.parse(getStream(PSI_EPSILON_SERIES), PSI_EPSILON_SERIES);

  327.             final PoissonSeries.CompiledSeries psiEpsilonSeries =
  328.                     PoissonSeries.compile(psiSeries, epsilonSeries);

  330.             return new TimeVectorFunction() {

  332.                 /** {@inheritDoc} */
  333.                 @Override
  334.                 public double[] value(final AbsoluteDate date) {
  335.                     final BodiesElements elements = arguments.evaluateAll(date);
  336.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  337.                     return new double[] {
  338.                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements)
  339.                     };
  340.                 }

  342.                 /** {@inheritDoc} */
  343.                 @Override
  344.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  345.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  346.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  347.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  348.                     result[0] = psiEpsilon[0];
  349.                     result[1] = psiEpsilon[1];
  350.                     result[2] = IAU1994ResolutionC7.value(elements);
  351.                     return result;
  352.                 }

  354.             };

  356.         }

  358.         /** {@inheritDoc} */
  359.         @Override
  360.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1)
  361.             throws OrekitException {

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

  366.             // constants from IERS 1996 page 21
  367.             // the underlying model is IAU 1982 GMST-UT1
  368.             final AbsoluteDate gmstReference =
  369.                 new AbsoluteDate(DateComponents.J2000_EPOCH, TimeComponents.H12, TimeScalesFactory.getTAI());
  370.             final double gmst0 = 24110.54841;
  371.             final double gmst1 = 8640184.812866;
  372.             final double gmst2 = 0.093104;
  373.             final double gmst3 = -6.2e-6;

  375.             return new TimeScalarFunction() {

  377.                 /** {@inheritDoc} */
  378.                 @Override
  379.                 public double value(final AbsoluteDate date) {

  381.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  382.                     final double dtai = date.durationFrom(gmstReference);
  383.                     final double tut1 = dtai + ut1.offsetFromTAI(date);
  384.                     final double tt   = tut1 / Constants.JULIAN_CENTURY;

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

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

  393.                 }

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

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

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

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

  411.                 }

  413.             };

  415.         }

  417.         /** {@inheritDoc} */
  418.         @Override
  419.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1)
  420.             throws OrekitException {

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

  425.             // constants from IERS 1996 page 21
  426.             // the underlying model is IAU 1982 GMST-UT1
  427.             final AbsoluteDate gmstReference =
  428.                 new AbsoluteDate(DateComponents.J2000_EPOCH, TimeComponents.H12, TimeScalesFactory.getTAI());
  429.             final double gmst1 = 8640184.812866;
  430.             final double gmst2 = 0.093104;
  431.             final double gmst3 = -6.2e-6;

  433.             return new TimeScalarFunction() {

  435.                 /** {@inheritDoc} */
  436.                 @Override
  437.                 public double value(final AbsoluteDate date) {

  439.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  440.                     final double dtai = date.durationFrom(gmstReference);
  441.                     final double tut1 = dtai + ut1.offsetFromTAI(date);
  442.                     final double tt   = tut1 / Constants.JULIAN_CENTURY;

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

  447.                 }

  449.                 /** {@inheritDoc} */
  450.                 @Override
  451.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  453.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  454.                     final T dtai = date.durationFrom(gmstReference);
  455.                     final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
  456.                     final T tt   = tut1.divide(Constants.JULIAN_CENTURY);

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

  461.                 }

  463.             };

  465.         }

  467.         /** {@inheritDoc} */
  468.         @Override
  469.         public TimeScalarFunction getGASTFunction(final TimeScale ut1, final EOPHistory eopHistory)
  470.             throws OrekitException {

  472.             // obliquity
  473.             final TimeScalarFunction epsilonA = getMeanObliquityFunction();

  475.             // GMST function
  476.             final TimeScalarFunction gmst = getGMSTFunction(ut1);

  478.             // nutation function
  479.             final TimeVectorFunction nutation = getNutationFunction();

  481.             return new TimeScalarFunction() {

  483.                 /** {@inheritDoc} */
  484.                 @Override
  485.                 public double value(final AbsoluteDate date) {

  487.                     // compute equation of equinoxes
  488.                     final double[] angles = nutation.value(date);
  489.                     double deltaPsi = angles[0];
  490.                     if (eopHistory != null) {
  491.                         deltaPsi += eopHistory.getEquinoxNutationCorrection(date)[0];
  492.                     }
  493.                     final double eqe = deltaPsi  * FastMath.cos(epsilonA.value(date)) + angles[2];

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

  498.                 }

  500.                 /** {@inheritDoc} */
  501.                 @Override
  502.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  504.                     // compute equation of equinoxes
  505.                     final T[] angles = nutation.value(date);
  506.                     T deltaPsi = angles[0];
  507.                     if (eopHistory != null) {
  508.                         deltaPsi = deltaPsi.add(eopHistory.getEquinoxNutationCorrection(date)[0]);
  509.                     }
  510.                     final T eqe = deltaPsi.multiply(epsilonA.value(date).cos()).add(angles[2]);

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

  515.                 }

  517.             };

  519.         }

  521.         /** {@inheritDoc} */
  522.         @Override
  523.         public TimeVectorFunction getEOPTidalCorrection()
  524.             throws OrekitException {

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

  533.             // set up Poisson series
  534.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  535.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(17).
  536.                     withOptionalColumn(1).
  537.                     withGamma(7).
  538.                     withFirstDelaunay(2);
  539.             final PoissonSeries xSeries =
  540.                     xyParser.
  541.                     withSinCos(0, 14, milliAS, 15, milliAS).
  542.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
  543.             final PoissonSeries ySeries =
  544.                     xyParser.
  545.                     withSinCos(0, 16, milliAS, 17, milliAS).
  546.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES),
  547.                                                          TIDAL_CORRECTION_XP_YP_SERIES);

  549.             final double deciMilliS = 1.0e-4;
  550.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(17).
  551.                     withOptionalColumn(1).
  552.                     withGamma(7).
  553.                     withFirstDelaunay(2).
  554.                     withSinCos(0, 16, deciMilliS, 17, deciMilliS);
  555.             final PoissonSeries ut1Series =
  556.                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);

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

  560.         }

  562.         /** {@inheritDoc} */
  563.         public LoveNumbers getLoveNumbers() throws OrekitException {
  564.             return loadLoveNumbers(LOVE_NUMBERS);
  565.         }

  567.         /** {@inheritDoc} */
  568.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1)
  569.             throws OrekitException {

  571.             // set up nutation arguments
  572.             final FundamentalNutationArguments arguments = getNutationArguments(ut1);

  574.             // set up Poisson series
  575.             final PoissonSeriesParser k20Parser =
  576.                     new PoissonSeriesParser(18).
  577.                         withOptionalColumn(1).
  578.                         withDoodson(4, 2).
  579.                         withFirstDelaunay(10);
  580.             final PoissonSeriesParser k21Parser =
  581.                     new PoissonSeriesParser(18).
  582.                         withOptionalColumn(1).
  583.                         withDoodson(4, 3).
  584.                         withFirstDelaunay(10);
  585.             final PoissonSeriesParser k22Parser =
  586.                     new PoissonSeriesParser(16).
  587.                         withOptionalColumn(1).
  588.                         withDoodson(4, 2).
  589.                         withFirstDelaunay(10);

  591.             final double pico = 1.0e-12;
  592.             final PoissonSeries c20Series =
  593.                     k20Parser.
  594.                   withSinCos(0, 18, -pico, 16, pico).
  595.                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
  596.             final PoissonSeries c21Series =
  597.                     k21Parser.
  598.                     withSinCos(0, 17, pico, 18, pico).
  599.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  600.             final PoissonSeries s21Series =
  601.                     k21Parser.
  602.                     withSinCos(0, 18, -pico, 17, pico).
  603.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  604.             final PoissonSeries c22Series =
  605.                     k22Parser.
  606.                     withSinCos(0, -1, pico, 16, pico).
  607.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
  608.             final PoissonSeries s22Series =
  609.                     k22Parser.
  610.                     withSinCos(0, 16, -pico, -1, pico).
  611.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);

  613.             return new TideFrequencyDependenceFunction(arguments,
  614.                                                        c20Series,
  615.                                                        c21Series, s21Series,
  616.                                                        c22Series, s22Series);

  618.         }

  620.         /** {@inheritDoc} */
  621.         @Override
  622.         public double getPermanentTide() throws OrekitException {
  623.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  624.         }

  626.         /** {@inheritDoc} */
  627.         @Override
  628.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {

  630.             // constants from IERS 1996 page 47
  631.             final double globalFactor = -1.348e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  632.             final double coupling     =  0.00112;

  634.             return new TimeVectorFunction() {

  636.                 /** {@inheritDoc} */
  637.                 @Override
  638.                 public double[] value(final AbsoluteDate date) {
  639.                     final PoleCorrection pole = eopHistory.getPoleCorrection(date);
  640.                     return new double[] {
  641.                         globalFactor * (pole.getXp() + coupling * pole.getYp()),
  642.                         globalFactor * (coupling * pole.getXp() - pole.getYp()),
  643.                     };
  644.                 }

  646.                 /** {@inheritDoc} */
  647.                 @Override
  648.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  649.                     final FieldPoleCorrection<T> pole = eopHistory.getPoleCorrection(date);
  650.                     final T[] a = MathArrays.buildArray(date.getField(), 2);
  651.                     a[0] = pole.getXp().add(pole.getYp().multiply(coupling)).multiply(globalFactor);
  652.                     a[1] = pole.getXp().multiply(coupling).subtract(pole.getYp()).multiply(globalFactor);
  653.                     return a;
  654.                 }

  656.             };
  657.         }

  659.         /** {@inheritDoc} */
  660.         @Override
  661.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory)
  662.             throws OrekitException {

  664.             return new TimeVectorFunction() {

  666.                 /** {@inheritDoc} */
  667.                 @Override
  668.                 public double[] value(final AbsoluteDate date) {
  669.                     // there are no model for ocean pole tide prior to conventions 2010
  670.                     return new double[] {
  671.                         0.0, 0.0
  672.                     };
  673.                 }

  675.                 /** {@inheritDoc} */
  676.                 @Override
  677.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  678.                     // there are no model for ocean pole tide prior to conventions 2010
  679.                     return MathArrays.buildArray(date.getField(), 2);
  680.                 }

  682.             };
  683.         }

  685.     },

  687.     /** Constant for IERS 2003 conventions. */
  688.     IERS_2003 {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  732.         /** {@inheritDoc} */
  733.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
  734.             throws OrekitException {
  735.             return new FundamentalNutationArguments(this, timeScale,
  736.                                                     getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
  737.         }

  739.         /** {@inheritDoc} */
  740.         @Override
  741.         public TimeScalarFunction getMeanObliquityFunction() throws OrekitException {

  743.             // epsilon 0 value from chapter 5, page 41, other terms from equation 32 page 45
  744.             final PolynomialNutation epsilonA =
  745.                     new PolynomialNutation(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  746.                                              -46.84024  * Constants.ARC_SECONDS_TO_RADIANS,
  747.                                               -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  748.                                                0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  750.             return new TimeScalarFunction() {

  752.                 /** {@inheritDoc} */
  753.                 @Override
  754.                 public double value(final AbsoluteDate date) {
  755.                     return epsilonA.value(evaluateTC(date));
  756.                 }

  758.                 /** {@inheritDoc} */
  759.                 @Override
  760.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  761.                     return epsilonA.value(evaluateTC(date));
  762.                 }

  764.             };

  766.         }

  768.         /** {@inheritDoc} */
  769.         @Override
  770.         public TimeVectorFunction getXYSpXY2Function()
  771.             throws OrekitException {

  773.             // set up nutation arguments
  774.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  776.             // set up Poisson series
  777.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  778.             final PoissonSeriesParser parser =
  779.                     new PoissonSeriesParser(17).
  780.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  781.                         withFirstDelaunay(4).
  782.                         withFirstPlanetary(9).
  783.                         withSinCos(0, 2, microAS, 3, microAS);

  785.             final PoissonSeries xSeries = parser.parse(getStream(X_SERIES), X_SERIES);
  786.             final PoissonSeries ySeries = parser.parse(getStream(Y_SERIES), Y_SERIES);
  787.             final PoissonSeries sSeries = parser.parse(getStream(S_SERIES), S_SERIES);
  788.             final PoissonSeries.CompiledSeries xys = PoissonSeries.compile(xSeries, ySeries, sSeries);

  790.             // create a function evaluating the series
  791.             return new TimeVectorFunction() {

  793.                 /** {@inheritDoc} */
  794.                 @Override
  795.                 public double[] value(final AbsoluteDate date) {
  796.                     return xys.value(arguments.evaluateAll(date));
  797.                 }

  799.                 /** {@inheritDoc} */
  800.                 @Override
  801.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  802.                     return xys.value(arguments.evaluateAll(date));
  803.                 }

  805.             };

  807.         }


  810.         /** {@inheritDoc} */
  811.         @Override
  812.         public TimeVectorFunction getPrecessionFunction() throws OrekitException {

  814.             // set up the conventional polynomials
  815.             // the following values are from equation 32 in IERS 2003 conventions
  816.             final PolynomialNutation psiA =
  817.                     new PolynomialNutation(    0.0,
  818.                                             5038.47875   * Constants.ARC_SECONDS_TO_RADIANS,
  819.                                               -1.07259   * Constants.ARC_SECONDS_TO_RADIANS,
  820.                                               -0.001147  * Constants.ARC_SECONDS_TO_RADIANS);
  821.             final PolynomialNutation omegaA =
  822.                     new PolynomialNutation(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
  823.                                            -0.02524   * Constants.ARC_SECONDS_TO_RADIANS,
  824.                                             0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
  825.                                            -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
  826.             final PolynomialNutation chiA =
  827.                     new PolynomialNutation( 0.0,
  828.                                            10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
  829.                                            -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
  830.                                            -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);

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

  834.         }

  836.         /** {@inheritDoc} */
  837.         @Override
  838.         public TimeVectorFunction getNutationFunction()
  839.             throws OrekitException {

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

  844.             // set up Poisson series
  845.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  846.             final PoissonSeriesParser luniSolarParser =
  847.                     new PoissonSeriesParser(14).withFirstDelaunay(1);
  848.             final PoissonSeriesParser luniSolarPsiParser =
  849.                     luniSolarParser.
  850.                     withSinCos(0, 7, milliAS, 11, milliAS).
  851.                     withSinCos(1, 8, milliAS, 12, milliAS);
  852.             final PoissonSeries psiLuniSolarSeries =
  853.                     luniSolarPsiParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);
  854.             final PoissonSeriesParser luniSolarEpsilonParser =
  855.                     luniSolarParser.
  856.                     withSinCos(0, 13, milliAS, 9, milliAS).
  857.                     withSinCos(1, 14, milliAS, 10, milliAS);
  858.             final PoissonSeries epsilonLuniSolarSeries =
  859.                     luniSolarEpsilonParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);

  861.             final PoissonSeriesParser planetaryParser =
  862.                     new PoissonSeriesParser(21).
  863.                         withFirstDelaunay(2).
  864.                         withFirstPlanetary(7);
  865.             final PoissonSeriesParser planetaryPsiParser =
  866.                     planetaryParser.withSinCos(0, 17, milliAS, 18, milliAS);
  867.             final PoissonSeries psiPlanetarySeries =
  868.                     planetaryPsiParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);
  869.             final PoissonSeriesParser planetaryEpsilonParser =
  870.                     planetaryParser.withSinCos(0, 19, milliAS, 20, milliAS);
  871.             final PoissonSeries epsilonPlanetarySeries =
  872.                     planetaryEpsilonParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);

  874.             final PoissonSeries.CompiledSeries luniSolarSeries =
  875.                     PoissonSeries.compile(psiLuniSolarSeries, epsilonLuniSolarSeries);
  876.             final PoissonSeries.CompiledSeries planetarySeries =
  877.                     PoissonSeries.compile(psiPlanetarySeries, epsilonPlanetarySeries);

  879.             return new TimeVectorFunction() {

  881.                 /** {@inheritDoc} */
  882.                 @Override
  883.                 public double[] value(final AbsoluteDate date) {
  884.                     final BodiesElements elements = arguments.evaluateAll(date);
  885.                     final double[] luniSolar = luniSolarSeries.value(elements);
  886.                     final double[] planetary = planetarySeries.value(elements);
  887.                     return new double[] {
  888.                         luniSolar[0] + planetary[0], luniSolar[1] + planetary[1],
  889.                         IAU1994ResolutionC7.value(elements)
  890.                     };
  891.                 }

  893.                 /** {@inheritDoc} */
  894.                 @Override
  895.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  896.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  897.                     final T[] luniSolar = luniSolarSeries.value(elements);
  898.                     final T[] planetary = planetarySeries.value(elements);
  899.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  900.                     result[0] = luniSolar[0].add(planetary[0]);
  901.                     result[1] = luniSolar[1].add(planetary[1]);
  902.                     result[2] = IAU1994ResolutionC7.value(elements);
  903.                     return result;
  904.                 }

  906.             };

  908.         }

  910.         /** {@inheritDoc} */
  911.         @Override
  912.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1)
  913.             throws OrekitException {

  915.             // Earth Rotation Angle
  916.             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);

  918.             // Polynomial part of the apparent sidereal time series
  919.             // which is the opposite of Equation of Origins (EO)
  920.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  921.             final PoissonSeriesParser parser =
  922.                     new PoissonSeriesParser(17).
  923.                         withFirstDelaunay(4).
  924.                         withFirstPlanetary(9).
  925.                         withSinCos(0, 2, microAS, 3, microAS).
  926.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  927.             final PolynomialNutation minusEO =
  928.                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();

  930.             // create a function evaluating the series
  931.             return new TimeScalarFunction() {

  933.                 /** {@inheritDoc} */
  934.                 @Override
  935.                 public double value(final AbsoluteDate date) {
  936.                     return era.value(date) + minusEO.value(evaluateTC(date));
  937.                 }

  939.                 /** {@inheritDoc} */
  940.                 @Override
  941.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  942.                     return era.value(date).add(minusEO.value(evaluateTC(date)));
  943.                 }

  945.             };

  947.         }

  949.         /** {@inheritDoc} */
  950.         @Override
  951.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1)
  952.             throws OrekitException {

  954.             // Earth Rotation Angle
  955.             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);

  957.             // Polynomial part of the apparent sidereal time series
  958.             // which is the opposite of Equation of Origins (EO)
  959.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  960.             final PoissonSeriesParser parser =
  961.                     new PoissonSeriesParser(17).
  962.                         withFirstDelaunay(4).
  963.                         withFirstPlanetary(9).
  964.                         withSinCos(0, 2, microAS, 3, microAS).
  965.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  966.             final PolynomialNutation minusEO =
  967.                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();

  969.             // create a function evaluating the series
  970.             return new TimeScalarFunction() {

  972.                 /** {@inheritDoc} */
  973.                 @Override
  974.                 public double value(final AbsoluteDate date) {
  975.                     return era.getRate() + minusEO.derivative(evaluateTC(date));
  976.                 }

  978.                 /** {@inheritDoc} */
  979.                 @Override
  980.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  981.                     return minusEO.derivative(evaluateTC(date)).add(era.getRate());
  982.                 }

  984.             };

  986.         }

  988.         /** {@inheritDoc} */
  989.         @Override
  990.         public TimeScalarFunction getGASTFunction(final TimeScale ut1, final EOPHistory eopHistory)
  991.             throws OrekitException {

  993.             // set up nutation arguments
  994.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  996.             // mean obliquity function
  997.             final TimeScalarFunction epsilon = getMeanObliquityFunction();

  999.             // set up Poisson series
  1000.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  1001.             final PoissonSeriesParser luniSolarPsiParser =
  1002.                     new PoissonSeriesParser(14).
  1003.                         withFirstDelaunay(1).
  1004.                         withSinCos(0, 7, milliAS, 11, milliAS).
  1005.                         withSinCos(1, 8, milliAS, 12, milliAS);
  1006.             final PoissonSeries psiLuniSolarSeries =
  1007.                     luniSolarPsiParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);

  1009.             final PoissonSeriesParser planetaryPsiParser =
  1010.                     new PoissonSeriesParser(21).
  1011.                         withFirstDelaunay(2).
  1012.                         withFirstPlanetary(7).
  1013.                         withSinCos(0, 17, milliAS, 18, milliAS);
  1014.             final PoissonSeries psiPlanetarySeries =
  1015.                     planetaryPsiParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);

  1017.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1018.             final PoissonSeriesParser gstParser =
  1019.                     new PoissonSeriesParser(17).
  1020.                         withFirstDelaunay(4).
  1021.                         withFirstPlanetary(9).
  1022.                         withSinCos(0, 2, microAS, 3, microAS).
  1023.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1024.             final PoissonSeries gstSeries = gstParser.parse(getStream(GST_SERIES), GST_SERIES);
  1025.             final PoissonSeries.CompiledSeries psiGstSeries =
  1026.                     PoissonSeries.compile(psiLuniSolarSeries, psiPlanetarySeries, gstSeries);

  1028.             // ERA function
  1029.             final TimeScalarFunction era = getEarthOrientationAngleFunction(ut1);

  1031.             return new TimeScalarFunction() {

  1033.                 /** {@inheritDoc} */
  1034.                 @Override
  1035.                 public double value(final AbsoluteDate date) {

  1037.                     // evaluate equation of origins
  1038.                     final BodiesElements elements = arguments.evaluateAll(date);
  1039.                     final double[] angles = psiGstSeries.value(elements);
  1040.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  1041.                     final double deltaPsi = angles[0] + angles[1] + ddPsi;
  1042.                     final double epsilonA = epsilon.value(date);

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

  1048.                 }

  1050.                 /** {@inheritDoc} */
  1051.                 @Override
  1052.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  1054.                     // evaluate equation of origins
  1055.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1056.                     final T[] angles = psiGstSeries.value(elements);
  1057.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  1058.                     final T deltaPsi = angles[0].add(angles[1]).add(ddPsi);
  1059.                     final T epsilonA = epsilon.value(date);

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

  1065.                 }

  1067.             };

  1069.         }

  1071.         /** {@inheritDoc} */
  1072.         @Override
  1073.         public TimeVectorFunction getEOPTidalCorrection()
  1074.             throws OrekitException {

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

  1083.             // set up Poisson series
  1084.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1085.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
  1086.                     withOptionalColumn(1).
  1087.                     withGamma(2).
  1088.                     withFirstDelaunay(3);
  1089.             final PoissonSeries xSeries =
  1090.                     xyParser.
  1091.                     withSinCos(0, 10, microAS, 11, microAS).
  1092.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
  1093.             final PoissonSeries ySeries =
  1094.                     xyParser.
  1095.                     withSinCos(0, 12, microAS, 13, microAS).
  1096.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);

  1098.             final double microS = 1.0e-6;
  1099.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  1100.                     withOptionalColumn(1).
  1101.                     withGamma(2).
  1102.                     withFirstDelaunay(3).
  1103.                     withSinCos(0, 10, microS, 11, microS);
  1104.             final PoissonSeries ut1Series =
  1105.                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);

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

  1109.         }

  1111.         /** {@inheritDoc} */
  1112.         public LoveNumbers getLoveNumbers() throws OrekitException {
  1113.             return loadLoveNumbers(LOVE_NUMBERS);
  1114.         }

  1116.         /** {@inheritDoc} */
  1117.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1)
  1118.             throws OrekitException {

  1120.             // set up nutation arguments
  1121.             final FundamentalNutationArguments arguments = getNutationArguments(ut1);

  1123.             // set up Poisson series
  1124.             final PoissonSeriesParser k20Parser =
  1125.                     new PoissonSeriesParser(18).
  1126.                         withOptionalColumn(1).
  1127.                         withDoodson(4, 2).
  1128.                         withFirstDelaunay(10);
  1129.             final PoissonSeriesParser k21Parser =
  1130.                     new PoissonSeriesParser(18).
  1131.                         withOptionalColumn(1).
  1132.                         withDoodson(4, 3).
  1133.                         withFirstDelaunay(10);
  1134.             final PoissonSeriesParser k22Parser =
  1135.                     new PoissonSeriesParser(16).
  1136.                         withOptionalColumn(1).
  1137.                         withDoodson(4, 2).
  1138.                         withFirstDelaunay(10);

  1140.             final double pico = 1.0e-12;
  1141.             final PoissonSeries c20Series =
  1142.                     k20Parser.
  1143.                   withSinCos(0, 18, -pico, 16, pico).
  1144.                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
  1145.             final PoissonSeries c21Series =
  1146.                     k21Parser.
  1147.                     withSinCos(0, 17, pico, 18, pico).
  1148.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  1149.             final PoissonSeries s21Series =
  1150.                     k21Parser.
  1151.                     withSinCos(0, 18, -pico, 17, pico).
  1152.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  1153.             final PoissonSeries c22Series =
  1154.                     k22Parser.
  1155.                     withSinCos(0, -1, pico, 16, pico).
  1156.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
  1157.             final PoissonSeries s22Series =
  1158.                     k22Parser.
  1159.                     withSinCos(0, 16, -pico, -1, pico).
  1160.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);

  1162.             return new TideFrequencyDependenceFunction(arguments,
  1163.                                                        c20Series,
  1164.                                                        c21Series, s21Series,
  1165.                                                        c22Series, s22Series);

  1167.         }

  1169.         /** {@inheritDoc} */
  1170.         @Override
  1171.         public double getPermanentTide() throws OrekitException {
  1172.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1173.         }

  1175.         /** {@inheritDoc} */
  1176.         @Override
  1177.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory)
  1178.             throws OrekitException {

  1180.             // annual pole from ftp://tai.bipm.org/iers/conv2003/chapter7/annual.pole
  1181.             final TimeScale utc = TimeScalesFactory.getUTC();
  1182.             final SimpleTimeStampedTableParser.RowConverter<MeanPole> converter =
  1183.                 new SimpleTimeStampedTableParser.RowConverter<MeanPole>() {
  1184.                     /** {@inheritDoc} */
  1185.                     @Override
  1186.                     public MeanPole convert(final double[] rawFields) throws OrekitException {
  1187.                         return new MeanPole(new AbsoluteDate((int) rawFields[0], 1, 1, utc),
  1188.                                             rawFields[1] * Constants.ARC_SECONDS_TO_RADIANS,
  1189.                                             rawFields[2] * Constants.ARC_SECONDS_TO_RADIANS);
  1190.                     }
  1191.                 };
  1192.             final SimpleTimeStampedTableParser<MeanPole> parser =
  1193.                     new SimpleTimeStampedTableParser<MeanPole>(3, converter);
  1194.             final List<MeanPole> annualPoleList = parser.parse(getStream(ANNUAL_POLE), ANNUAL_POLE);
  1195.             final AbsoluteDate firstAnnualPoleDate = annualPoleList.get(0).getDate();
  1196.             final AbsoluteDate lastAnnualPoleDate  = annualPoleList.get(annualPoleList.size() - 1).getDate();
  1197.             final ImmutableTimeStampedCache<MeanPole> annualCache =
  1198.                     new ImmutableTimeStampedCache<MeanPole>(2, annualPoleList);

  1200.             // polynomial extension from IERS 2003, section 7.1.4, equations 23a and 23b
  1201.             final double xp0    = 0.054   * Constants.ARC_SECONDS_TO_RADIANS;
  1202.             final double xp0Dot = 0.00083 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
  1203.             final double yp0    = 0.357   * Constants.ARC_SECONDS_TO_RADIANS;
  1204.             final double yp0Dot = 0.00395 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;

  1206.             // constants from IERS 2003, section 6.2
  1207.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1208.             final double ratio        =  0.00115;

  1210.             return new TimeVectorFunction() {

  1212.                 /** {@inheritDoc} */
  1213.                 @Override
  1214.                 public double[] value(final AbsoluteDate date) {

  1216.                     // we can't compute anything before the range covered by the annual pole file
  1217.                     if (date.compareTo(firstAnnualPoleDate) <= 0) {
  1218.                         return new double[] {
  1219.                             0.0, 0.0
  1220.                         };
  1221.                     }

  1223.                     // evaluate mean pole
  1224.                     double meanPoleX = 0;
  1225.                     double meanPoleY = 0;
  1226.                     if (date.compareTo(lastAnnualPoleDate) <= 0) {
  1227.                         // we are within the range covered by the annual pole file,
  1228.                         // we interpolate within it
  1229.                         try {
  1230.                             final HermiteInterpolator interpolator = new HermiteInterpolator();
  1231.                             annualCache.getNeighbors(date).forEach(neighbor ->
  1232.                                 interpolator.addSamplePoint(neighbor.getDate().durationFrom(date),
  1233.                                                             new double[] {
  1234.                                                                 neighbor.getX(), neighbor.getY()
  1235.                                                             }));
  1236.                             final double[] interpolated = interpolator.value(0);
  1237.                             meanPoleX = interpolated[0];
  1238.                             meanPoleY = interpolated[1];
  1239.                         } catch (TimeStampedCacheException tsce) {
  1240.                             // this should never happen
  1241.                             throw new OrekitInternalError(tsce);
  1242.                         }
  1243.                     } else {

  1245.                         // we are after the range covered by the annual pole file,
  1246.                         // we use the polynomial extension
  1247.                         final double t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
  1248.                         meanPoleX = xp0 + t * xp0Dot;
  1249.                         meanPoleY = yp0 + t * yp0Dot;

  1251.                     }

  1253.                     // evaluate wobble variables
  1254.                     final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1255.                     final double m1 = correction.getXp() - meanPoleX;
  1256.                     final double m2 = meanPoleY - correction.getYp();

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

  1289.                 }

  1291.                 /** {@inheritDoc} */
  1292.                 @Override
  1293.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

  1302.                     // evaluate mean pole
  1303.                     T meanPoleX = date.getField().getZero();
  1304.                     T meanPoleY = date.getField().getZero();
  1305.                     if (aDate.compareTo(lastAnnualPoleDate) <= 0) {
  1306.                         // we are within the range covered by the annual pole file,
  1307.                         // we interpolate within it
  1308.                         try {
  1309.                             final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
  1310.                             final T[] y = MathArrays.buildArray(date.getField(), 2);
  1311.                             final T zero = date.getField().getZero();
  1312.                             final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
  1313.                                                                                                        // for example removing derivatives
  1314.                                                                                                        // if T was DerivativeStructure
  1315.                             annualCache.getNeighbors(aDate).forEach(neighbor -> {
  1316.                                 y[0] = zero.add(neighbor.getX());
  1317.                                 y[1] = zero.add(neighbor.getY());
  1318.                                 interpolator.addSamplePoint(central.durationFrom(neighbor.getDate()).negate(), y);
  1319.                             });
  1320.                             final T[] interpolated = interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
  1321.                             meanPoleX = interpolated[0];
  1322.                             meanPoleY = interpolated[1];
  1323.                         } catch (TimeStampedCacheException tsce) {
  1324.                             // this should never happen
  1325.                             throw new OrekitInternalError(tsce);
  1326.                         }
  1327.                     } else {

  1329.                         // we are after the range covered by the annual pole file,
  1330.                         // we use the polynomial extension
  1331.                         final T t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
  1332.                         meanPoleX = t.multiply(xp0Dot).add(xp0);
  1333.                         meanPoleY = t.multiply(yp0Dot).add(yp0);

  1335.                     }

  1337.                     // evaluate wobble variables
  1338.                     final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1339.                     final T m1 = correction.getXp().subtract(meanPoleX);
  1340.                     final T m2 = meanPoleY.subtract(correction.getYp());

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

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

  1373.                     return a;

  1375.                 }

  1377.             };

  1379.         }

  1381.         /** {@inheritDoc} */
  1382.         @Override
  1383.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory)
  1384.             throws OrekitException {

  1386.             return new TimeVectorFunction() {

  1388.                 /** {@inheritDoc} */
  1389.                 @Override
  1390.                 public double[] value(final AbsoluteDate date) {
  1391.                     // there are no model for ocean pole tide prior to conventions 2010
  1392.                     return new double[] {
  1393.                         0.0, 0.0
  1394.                     };
  1395.                 }

  1397.                 /** {@inheritDoc} */
  1398.                 @Override
  1399.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1400.                     // there are no model for ocean pole tide prior to conventions 2010
  1401.                     return MathArrays.buildArray(date.getField(), 2);
  1402.                 }

  1404.             };
  1405.         }

  1407.     },

  1409.     /** Constant for IERS 2010 conventions. */
  1410.     IERS_2010 {

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

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

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

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

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

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

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

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

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

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

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

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

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

  1451.         /** {@inheritDoc} */
  1452.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
  1453.             throws OrekitException {
  1454.             return new FundamentalNutationArguments(this, timeScale,
  1455.                                                     getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
  1456.         }

  1458.         /** {@inheritDoc} */
  1459.         @Override
  1460.         public TimeScalarFunction getMeanObliquityFunction() throws OrekitException {

  1462.             // epsilon 0 value from chapter 5, page 56, other terms from equation 5.40 page 65
  1463.             final PolynomialNutation epsilonA =
  1464.                     new PolynomialNutation(84381.406        * Constants.ARC_SECONDS_TO_RADIANS,
  1465.                                              -46.836769     * Constants.ARC_SECONDS_TO_RADIANS,
  1466.                                               -0.0001831    * Constants.ARC_SECONDS_TO_RADIANS,
  1467.                                                0.00200340   * Constants.ARC_SECONDS_TO_RADIANS,
  1468.                                               -0.000000576  * Constants.ARC_SECONDS_TO_RADIANS,
  1469.                                               -0.0000000434 * Constants.ARC_SECONDS_TO_RADIANS);

  1471.             return new TimeScalarFunction() {

  1473.                 /** {@inheritDoc} */
  1474.                 @Override
  1475.                 public double value(final AbsoluteDate date) {
  1476.                     return epsilonA.value(evaluateTC(date));
  1477.                 }

  1479.                 /** {@inheritDoc} */
  1480.                 @Override
  1481.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1482.                     return epsilonA.value(evaluateTC(date));
  1483.                 }

  1485.             };

  1487.         }

  1489.         /** {@inheritDoc} */
  1490.         @Override
  1491.         public TimeVectorFunction getXYSpXY2Function() throws OrekitException {

  1493.             // set up nutation arguments
  1494.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  1496.             // set up Poisson series
  1497.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1498.             final PoissonSeriesParser parser =
  1499.                     new PoissonSeriesParser(17).
  1500.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  1501.                         withFirstDelaunay(4).
  1502.                         withFirstPlanetary(9).
  1503.                         withSinCos(0, 2, microAS, 3, microAS);
  1504.             final PoissonSeries xSeries = parser.parse(getStream(X_SERIES), X_SERIES);
  1505.             final PoissonSeries ySeries = parser.parse(getStream(Y_SERIES), Y_SERIES);
  1506.             final PoissonSeries sSeries = parser.parse(getStream(S_SERIES), S_SERIES);
  1507.             final PoissonSeries.CompiledSeries xys = PoissonSeries.compile(xSeries, ySeries, sSeries);

  1509.             // create a function evaluating the series
  1510.             return new TimeVectorFunction() {

  1512.                 /** {@inheritDoc} */
  1513.                 @Override
  1514.                 public double[] value(final AbsoluteDate date) {
  1515.                     return xys.value(arguments.evaluateAll(date));
  1516.                 }

  1518.                 /** {@inheritDoc} */
  1519.                 @Override
  1520.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1521.                     return xys.value(arguments.evaluateAll(date));
  1522.                 }

  1524.             };

  1526.         }

  1528.         /** {@inheritDoc} */
  1529.         public LoveNumbers getLoveNumbers() throws OrekitException {
  1530.             return loadLoveNumbers(LOVE_NUMBERS);
  1531.         }

  1533.         /** {@inheritDoc} */
  1534.         public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1)
  1535.             throws OrekitException {

  1537.             // set up nutation arguments
  1538.             final FundamentalNutationArguments arguments = getNutationArguments(ut1);

  1540.             // set up Poisson series
  1541.             final PoissonSeriesParser k20Parser =
  1542.                     new PoissonSeriesParser(18).
  1543.                         withOptionalColumn(1).
  1544.                         withDoodson(4, 2).
  1545.                         withFirstDelaunay(10);
  1546.             final PoissonSeriesParser k21Parser =
  1547.                     new PoissonSeriesParser(18).
  1548.                         withOptionalColumn(1).
  1549.                         withDoodson(4, 3).
  1550.                         withFirstDelaunay(10);
  1551.             final PoissonSeriesParser k22Parser =
  1552.                     new PoissonSeriesParser(16).
  1553.                         withOptionalColumn(1).
  1554.                         withDoodson(4, 2).
  1555.                         withFirstDelaunay(10);

  1557.             final double pico = 1.0e-12;
  1558.             final PoissonSeries c20Series =
  1559.                     k20Parser.
  1560.                     withSinCos(0, 18, -pico, 16, pico).
  1561.                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
  1562.             final PoissonSeries c21Series =
  1563.                     k21Parser.
  1564.                     withSinCos(0, 17, pico, 18, pico).
  1565.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  1566.             final PoissonSeries s21Series =
  1567.                     k21Parser.
  1568.                     withSinCos(0, 18, -pico, 17, pico).
  1569.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  1570.             final PoissonSeries c22Series =
  1571.                     k22Parser.
  1572.                     withSinCos(0, -1, pico, 16, pico).
  1573.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
  1574.             final PoissonSeries s22Series =
  1575.                     k22Parser.
  1576.                     withSinCos(0, 16, -pico, -1, pico).
  1577.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);

  1579.             return new TideFrequencyDependenceFunction(arguments,
  1580.                                                        c20Series,
  1581.                                                        c21Series, s21Series,
  1582.                                                        c22Series, s22Series);

  1584.         }

  1586.         /** {@inheritDoc} */
  1587.         @Override
  1588.         public double getPermanentTide() throws OrekitException {
  1589.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1590.         }

  1592.         /** Compute pole wobble variables m₁ and m₂.
  1593.          * @param date current date
  1594.          * @param eopHistory EOP history
  1595.          * @return array containing m₁ and m₂
  1596.          */
  1597.         private double[] computePoleWobble(final AbsoluteDate date, final EOPHistory eopHistory) {

  1599.             // polynomial model from IERS 2010, table 7.7
  1600.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1601.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1602.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1603.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1604.             final AbsoluteDate changeDate = new AbsoluteDate(2010, 1, 1, TimeScalesFactory.getTT());

  1606.             // evaluate mean pole
  1607.             final double[] xPolynomial;
  1608.             final double[] yPolynomial;
  1609.             if (date.compareTo(changeDate) <= 0) {
  1610.                 xPolynomial = new double[] {
  1611.                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
  1612.                 };
  1613.                 yPolynomial = new double[] {
  1614.                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
  1615.                 };
  1616.             } else {
  1617.                 xPolynomial = new double[] {
  1618.                     23.513 * f0, 7.6141 * f1
  1619.                 };
  1620.                 yPolynomial = new double[] {
  1621.                     358.891 * f0,  -0.6287 * f1
  1622.                 };
  1623.             }
  1624.             double meanPoleX = 0;
  1625.             double meanPoleY = 0;
  1626.             final double t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
  1627.             for (int i = xPolynomial.length - 1; i >= 0; --i) {
  1628.                 meanPoleX = meanPoleX * t + xPolynomial[i];
  1629.             }
  1630.             for (int i = yPolynomial.length - 1; i >= 0; --i) {
  1631.                 meanPoleY = meanPoleY * t + yPolynomial[i];
  1632.             }

  1634.             // evaluate wobble variables
  1635.             final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1636.             final double m1 = correction.getXp() - meanPoleX;
  1637.             final double m2 = meanPoleY - correction.getYp();

  1639.             return new double[] {
  1640.                 m1, m2
  1641.             };

  1643.         }

  1645.         /** Compute pole wobble variables m₁ and m₂.
  1646.          * @param date current date
  1647.          * @param <T> type of the field elements
  1648.          * @param eopHistory EOP history
  1649.          * @return array containing m₁ and m₂
  1650.          */
  1651.         private <T extends RealFieldElement<T>> T[] computePoleWobble(final FieldAbsoluteDate<T> date, final EOPHistory eopHistory) {

  1653.             // polynomial model from IERS 2010, table 7.7
  1654.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1655.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1656.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1657.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1658.             final AbsoluteDate changeDate = new AbsoluteDate(2010, 1, 1, TimeScalesFactory.getTT());

  1660.             // evaluate mean pole
  1661.             final double[] xPolynomial;
  1662.             final double[] yPolynomial;
  1663.             if (date.toAbsoluteDate().compareTo(changeDate) <= 0) {
  1664.                 xPolynomial = new double[] {
  1665.                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
  1666.                 };
  1667.                 yPolynomial = new double[] {
  1668.                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
  1669.                 };
  1670.             } else {
  1671.                 xPolynomial = new double[] {
  1672.                     23.513 * f0, 7.6141 * f1
  1673.                 };
  1674.                 yPolynomial = new double[] {
  1675.                     358.891 * f0,  -0.6287 * f1
  1676.                 };
  1677.             }
  1678.             T meanPoleX = date.getField().getZero();
  1679.             T meanPoleY = date.getField().getZero();
  1680.             final T t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
  1681.             for (int i = xPolynomial.length - 1; i >= 0; --i) {
  1682.                 meanPoleX = meanPoleX.multiply(t).add(xPolynomial[i]);
  1683.             }
  1684.             for (int i = yPolynomial.length - 1; i >= 0; --i) {
  1685.                 meanPoleY = meanPoleY.multiply(t).add(yPolynomial[i]);
  1686.             }

  1688.             // evaluate wobble variables
  1689.             final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
  1690.             final T[] m = MathArrays.buildArray(date.getField(), 2);
  1691.             m[0] = correction.getXp().subtract(meanPoleX);
  1692.             m[1] = meanPoleY.subtract(correction.getYp());

  1694.             return m;

  1696.         }

  1698.         /** {@inheritDoc} */
  1699.         @Override
  1700.         public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory)
  1701.             throws OrekitException {

  1703.             // constants from IERS 2010, section 6.4
  1704.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1705.             final double ratio        =  0.00115;

  1707.             return new TimeVectorFunction() {

  1709.                 /** {@inheritDoc} */
  1710.                 @Override
  1711.                 public double[] value(final AbsoluteDate date) {

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

  1716.                     return new double[] {
  1717.                         // the following correspond to the equations published in IERS 2010 conventions,
  1718.                         // section 6.4 page 94. The equations read:
  1719.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1720.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1721.                         // These equations seem to fix what was probably a sign error in IERS 2003
  1722.                         // conventions section 6.2 page 65. In this older publication, the equations read:
  1723.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1724.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1725.                         globalFactor * (wobbleM[0] + ratio * wobbleM[1]),
  1726.                         globalFactor * (wobbleM[1] - ratio * wobbleM[0])
  1727.                     };

  1729.                 }

  1731.                 /** {@inheritDoc} */
  1732.                 @Override
  1733.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

  1740.                     // the following correspond to the equations published in IERS 2010 conventions,
  1741.                     // section 6.4 page 94. The equations read:
  1742.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1743.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1744.                     // These equations seem to fix what was probably a sign error in IERS 2003
  1745.                     // conventions section 6.2 page 65. In this older publication, the equations read:
  1746.                     // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1747.                     // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1748.                     a[0] = wobbleM[0].add(wobbleM[1].multiply( ratio)).multiply(globalFactor);
  1749.                     a[1] = wobbleM[1].add(wobbleM[0].multiply(-ratio)).multiply(globalFactor);

  1751.                     return a;

  1753.                 }

  1755.             };

  1757.         }

  1759.         /** {@inheritDoc} */
  1760.         @Override
  1761.         public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory)
  1762.             throws OrekitException {

  1764.             return new TimeVectorFunction() {

  1766.                 /** {@inheritDoc} */
  1767.                 @Override
  1768.                 public double[] value(final AbsoluteDate date) {

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

  1773.                     return new double[] {
  1774.                         // the following correspond to the equations published in IERS 2010 conventions,
  1775.                         // section 6.4 page 94 equation 6.24:
  1776.                         // ∆C₂₁ = −2.1778 × 10⁻¹⁰ (m₁ − 0.01724m₂)
  1777.                         // ∆S₂₁ = −1.7232 × 10⁻¹⁰ (m₂ − 0.03365m₁)
  1778.                         -2.1778e-10 * (wobbleM[0] - 0.01724 * wobbleM[1]) / Constants.ARC_SECONDS_TO_RADIANS,
  1779.                         -1.7232e-10 * (wobbleM[1] - 0.03365 * wobbleM[0]) / Constants.ARC_SECONDS_TO_RADIANS
  1780.                     };

  1782.                 }

  1784.                 /** {@inheritDoc} */
  1785.                 @Override
  1786.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

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

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

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

  1800.                     return v;

  1802.                 }

  1804.             };

  1806.         }

  1808.         /** {@inheritDoc} */
  1809.         @Override
  1810.         public TimeVectorFunction getPrecessionFunction() throws OrekitException {

  1812.             // set up the conventional polynomials
  1813.             // the following values are from equation 5.40 in IERS 2010 conventions
  1814.             final PolynomialNutation psiA =
  1815.                     new PolynomialNutation(   0.0,
  1816.                                            5038.481507     * Constants.ARC_SECONDS_TO_RADIANS,
  1817.                                              -1.0790069    * Constants.ARC_SECONDS_TO_RADIANS,
  1818.                                              -0.00114045   * Constants.ARC_SECONDS_TO_RADIANS,
  1819.                                               0.000132851  * Constants.ARC_SECONDS_TO_RADIANS,
  1820.                                              -0.0000000951 * Constants.ARC_SECONDS_TO_RADIANS);
  1821.             final PolynomialNutation omegaA =
  1822.                     new PolynomialNutation(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
  1823.                                            -0.025754     * Constants.ARC_SECONDS_TO_RADIANS,
  1824.                                             0.0512623    * Constants.ARC_SECONDS_TO_RADIANS,
  1825.                                            -0.00772503   * Constants.ARC_SECONDS_TO_RADIANS,
  1826.                                            -0.000000467  * Constants.ARC_SECONDS_TO_RADIANS,
  1827.                                             0.0000003337 * Constants.ARC_SECONDS_TO_RADIANS);
  1828.             final PolynomialNutation chiA =
  1829.                     new PolynomialNutation( 0.0,
  1830.                                            10.556403     * Constants.ARC_SECONDS_TO_RADIANS,
  1831.                                            -2.3814292    * Constants.ARC_SECONDS_TO_RADIANS,
  1832.                                            -0.00121197   * Constants.ARC_SECONDS_TO_RADIANS,
  1833.                                             0.000170663  * Constants.ARC_SECONDS_TO_RADIANS,
  1834.                                            -0.0000000560 * Constants.ARC_SECONDS_TO_RADIANS);

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

  1838.         }

  1840.          /** {@inheritDoc} */
  1841.         @Override
  1842.         public TimeVectorFunction getNutationFunction()
  1843.             throws OrekitException {

  1845.             // set up nutation arguments
  1846.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  1848.             // set up Poisson series
  1849.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1850.             final PoissonSeriesParser parser =
  1851.                     new PoissonSeriesParser(17).
  1852.                         withFirstDelaunay(4).
  1853.                         withFirstPlanetary(9).
  1854.                         withSinCos(0, 2, microAS, 3, microAS);
  1855.             final PoissonSeries psiSeries     = parser.parse(getStream(PSI_SERIES), PSI_SERIES);
  1856.             final PoissonSeries epsilonSeries = parser.parse(getStream(EPSILON_SERIES), EPSILON_SERIES);
  1857.             final PoissonSeries.CompiledSeries psiEpsilonSeries =
  1858.                     PoissonSeries.compile(psiSeries, epsilonSeries);

  1860.             return new TimeVectorFunction() {

  1862.                 /** {@inheritDoc} */
  1863.                 @Override
  1864.                 public double[] value(final AbsoluteDate date) {
  1865.                     final BodiesElements elements = arguments.evaluateAll(date);
  1866.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  1867.                     return new double[] {
  1868.                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements)
  1869.                     };
  1870.                 }

  1872.                 /** {@inheritDoc} */
  1873.                 @Override
  1874.                 public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1875.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  1876.                     final T[] psiEpsilon = psiEpsilonSeries.value(elements);
  1877.                     final T[] result = MathArrays.buildArray(date.getField(), 3);
  1878.                     result[0] = psiEpsilon[0];
  1879.                     result[1] = psiEpsilon[1];
  1880.                     result[2] = IAU1994ResolutionC7.value(elements);
  1881.                     return result;
  1882.                 }

  1884.             };

  1886.         }

  1888.         /** {@inheritDoc} */
  1889.         @Override
  1890.         public TimeScalarFunction getGMSTFunction(final TimeScale ut1) throws OrekitException {

  1892.             // Earth Rotation Angle
  1893.             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);

  1895.             // Polynomial part of the apparent sidereal time series
  1896.             // which is the opposite of Equation of Origins (EO)
  1897.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1898.             final PoissonSeriesParser parser =
  1899.                     new PoissonSeriesParser(17).
  1900.                         withFirstDelaunay(4).
  1901.                         withFirstPlanetary(9).
  1902.                         withSinCos(0, 2, microAS, 3, microAS).
  1903.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1904.             final PolynomialNutation minusEO =
  1905.                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();

  1907.             // create a function evaluating the series
  1908.             return new TimeScalarFunction() {

  1910.                 /** {@inheritDoc} */
  1911.                 @Override
  1912.                 public double value(final AbsoluteDate date) {
  1913.                     return era.value(date) + minusEO.value(evaluateTC(date));
  1914.                 }

  1916.                 /** {@inheritDoc} */
  1917.                 @Override
  1918.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1919.                     return era.value(date).add(minusEO.value(evaluateTC(date)));
  1920.                 }

  1922.             };

  1924.         }

  1926.         /** {@inheritDoc} */
  1927.         @Override
  1928.         public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1) throws OrekitException {

  1930.             // Earth Rotation Angle
  1931.             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);

  1933.             // Polynomial part of the apparent sidereal time series
  1934.             // which is the opposite of Equation of Origins (EO)
  1935.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1936.             final PoissonSeriesParser parser =
  1937.                     new PoissonSeriesParser(17).
  1938.                         withFirstDelaunay(4).
  1939.                         withFirstPlanetary(9).
  1940.                         withSinCos(0, 2, microAS, 3, microAS).
  1941.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1942.             final PolynomialNutation minusEO =
  1943.                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();

  1945.             // create a function evaluating the series
  1946.             return new TimeScalarFunction() {

  1948.                 /** {@inheritDoc} */
  1949.                 @Override
  1950.                 public double value(final AbsoluteDate date) {
  1951.                     return era.getRate() + minusEO.derivative(evaluateTC(date));
  1952.                 }

  1954.                 /** {@inheritDoc} */
  1955.                 @Override
  1956.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
  1957.                     return minusEO.derivative(evaluateTC(date)).add(era.getRate());
  1958.                 }

  1960.             };

  1962.         }

  1964.         /** {@inheritDoc} */
  1965.         @Override
  1966.         public TimeScalarFunction getGASTFunction(final TimeScale ut1, final EOPHistory eopHistory)
  1967.             throws OrekitException {

  1969.             // set up nutation arguments
  1970.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  1972.             // mean obliquity function
  1973.             final TimeScalarFunction epsilon = getMeanObliquityFunction();

  1975.             // set up Poisson series
  1976.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1977.             final PoissonSeriesParser baseParser =
  1978.                     new PoissonSeriesParser(17).
  1979.                         withFirstDelaunay(4).
  1980.                         withFirstPlanetary(9).
  1981.                         withSinCos(0, 2, microAS, 3, microAS);
  1982.             final PoissonSeriesParser gstParser  = baseParser.withPolynomialPart('t', Unit.ARC_SECONDS);
  1983.             final PoissonSeries psiSeries        = baseParser.parse(getStream(PSI_SERIES), PSI_SERIES);
  1984.             final PoissonSeries gstSeries        = gstParser.parse(getStream(GST_SERIES), GST_SERIES);
  1985.             final PoissonSeries.CompiledSeries psiGstSeries =
  1986.                     PoissonSeries.compile(psiSeries, gstSeries);

  1988.             // ERA function
  1989.             final TimeScalarFunction era = getEarthOrientationAngleFunction(ut1);

  1991.             return new TimeScalarFunction() {

  1993.                 /** {@inheritDoc} */
  1994.                 @Override
  1995.                 public double value(final AbsoluteDate date) {

  1997.                     // evaluate equation of origins
  1998.                     final BodiesElements elements = arguments.evaluateAll(date);
  1999.                     final double[] angles = psiGstSeries.value(elements);
  2000.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  2001.                     final double deltaPsi = angles[0] + ddPsi;
  2002.                     final double epsilonA = epsilon.value(date);

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

  2008.                 }

  2010.                 /** {@inheritDoc} */
  2011.                 @Override
  2012.                 public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  2014.                     // evaluate equation of origins
  2015.                     final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  2016.                     final T[] angles = psiGstSeries.value(elements);
  2017.                     final T ddPsi    = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
  2018.                     final T deltaPsi = angles[0].add(ddPsi);
  2019.                     final T epsilonA = epsilon.value(date);

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

  2025.                 }

  2027.             };

  2029.         }

  2031.         /** {@inheritDoc} */
  2032.         @Override
  2033.         public TimeVectorFunction getEOPTidalCorrection()
  2034.             throws OrekitException {

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

  2043.             // set up Poisson series
  2044.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  2045.             final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
  2046.                     withOptionalColumn(1).
  2047.                     withGamma(2).
  2048.                     withFirstDelaunay(3);
  2049.             final PoissonSeries xSeries =
  2050.                     xyParser.
  2051.                     withSinCos(0, 10, microAS, 11, microAS).
  2052.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
  2053.             final PoissonSeries ySeries =
  2054.                     xyParser.
  2055.                     withSinCos(0, 12, microAS, 13, microAS).
  2056.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);

  2058.             final double microS = 1.0e-6;
  2059.             final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
  2060.                     withOptionalColumn(1).
  2061.                     withGamma(2).
  2062.                     withFirstDelaunay(3).
  2063.                     withSinCos(0, 10, microS, 11, microS);
  2064.             final PoissonSeries ut1Series =
  2065.                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);

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

  2069.         }

  2071.     };

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

  2076.     /** Get the reference epoch for fundamental nutation arguments.
  2077.      * @return reference epoch for fundamental nutation arguments
  2078.      * @since 6.1
  2079.      */
  2080.     public AbsoluteDate getNutationReferenceEpoch() {
  2081.         // IERS 1996, IERS 2003 and IERS 2010 use the same J2000.0 reference date
  2082.         return AbsoluteDate.J2000_EPOCH;
  2083.     }

  2085.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2086.      * @param date current date
  2087.      * @return date offset in Julian centuries
  2088.      * @since 6.1
  2089.      */
  2090.     public double evaluateTC(final AbsoluteDate date) {
  2091.         return date.durationFrom(getNutationReferenceEpoch()) / Constants.JULIAN_CENTURY;
  2092.     }

  2094.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  2095.      * @param date current date
  2096.      * @param <T> type of the field elements
  2097.      * @return date offset in Julian centuries
  2098.      * @since 9.0
  2099.      */
  2100.     public <T extends RealFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date) {
  2101.         return date.durationFrom(getNutationReferenceEpoch()).divide(Constants.JULIAN_CENTURY);
  2102.     }

  2104.     /** Get the fundamental nutation arguments.
  2105.      * @param timeScale time scale for computing Greenwich Mean Sidereal Time
  2106.      * (typically {@link TimeScalesFactory#getUT1(IERSConventions, boolean) UT1})
  2107.      * @return fundamental nutation arguments
  2108.      * @exception OrekitException if fundamental nutation arguments cannot be loaded
  2109.      * @since 6.1
  2110.      */
  2111.     public abstract FundamentalNutationArguments getNutationArguments(TimeScale timeScale)
  2112.         throws OrekitException;

  2114.     /** Get the function computing mean obliquity of the ecliptic.
  2115.      * @return function computing mean obliquity of the ecliptic
  2116.      * @exception OrekitException if table cannot be loaded
  2117.      * @since 6.1
  2118.      */
  2119.     public abstract TimeScalarFunction getMeanObliquityFunction() throws OrekitException;

  2121.     /** Get the function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components.
  2122.      * <p>
  2123.      * The returned function computes the two X, Y components of CIP and the S+XY/2 component of the non-rotating CIO.
  2124.      * </p>
  2125.      * @return function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components
  2126.      * @exception OrekitException if table cannot be loaded
  2127.      * @since 6.1
  2128.      */
  2129.     public abstract TimeVectorFunction getXYSpXY2Function()
  2130.         throws OrekitException;

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


  2148.     /** Get the function computing the precession angles.
  2149.      * <p>
  2150.      * The function returned computes the three precession angles
  2151.      * ψ<sub>A</sub> (around Z axis), ω<sub>A</sub> (around X axis)
  2152.      * and χ<sub>A</sub> (around Z axis). The constant angle ε₀
  2153.      * for the fourth rotation (around X axis) can be retrieved by evaluating the
  2154.      * function returned by {@link #getMeanObliquityFunction()} at {@link
  2155.      * #getNutationReferenceEpoch() nutation reference epoch}.
  2156.      * </p>
  2157.      * @return function computing the precession angle
  2158.      * @exception OrekitException if table cannot be loaded
  2159.      * @since 6.1
  2160.      */
  2161.     public abstract TimeVectorFunction getPrecessionFunction() throws OrekitException;

  2163.     /** Get the function computing the nutation angles.
  2164.      * <p>
  2165.      * The function returned computes the two classical angles ΔΨ and Δε,
  2166.      * and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
  2167.      * resolution C7 (the correction is forced to 0 before this date)
  2168.      * </p>
  2169.      * @return function computing the nutation in longitude ΔΨ and Δε
  2170.      * and the correction of equation of equinoxes
  2171.      * @exception OrekitException if table cannot be loaded
  2172.      * @since 6.1
  2173.      */
  2174.     public abstract TimeVectorFunction getNutationFunction()
  2175.         throws OrekitException;

  2177.     /** Get the function computing Greenwich mean sidereal time, in radians.
  2178.      * @param ut1 UT1 time scale
  2179.      * @return function computing Greenwich mean sidereal time
  2180.      * @exception OrekitException if table cannot be loaded
  2181.      * @since 6.1
  2182.      */
  2183.     public abstract TimeScalarFunction getGMSTFunction(TimeScale ut1)
  2184.         throws OrekitException;

  2186.     /** Get the function computing Greenwich mean sidereal time rate, in radians per second.
  2187.      * @param ut1 UT1 time scale
  2188.      * @return function computing Greenwich mean sidereal time rate
  2189.      * @exception OrekitException if table cannot be loaded
  2190.      * @since 9.0
  2191.      */
  2192.     public abstract TimeScalarFunction getGMSTRateFunction(TimeScale ut1)
  2193.         throws OrekitException;

  2195.     /** Get the function computing Greenwich apparent sidereal time, in radians.
  2196.      * @param ut1 UT1 time scale
  2197.      * @param eopHistory EOP history
  2198.      * @return function computing Greenwich apparent sidereal time
  2199.      * @exception OrekitException if table cannot be loaded
  2200.      * @since 6.1
  2201.      */
  2202.     public abstract TimeScalarFunction getGASTFunction(TimeScale ut1, EOPHistory eopHistory)
  2203.         throws OrekitException;

  2205.     /** Get the function computing tidal corrections for Earth Orientation Parameters.
  2206.      * @return function computing tidal corrections for Earth Orientation Parameters,
  2207.      * for xp, yp, ut1 and lod respectively
  2208.      * @exception OrekitException if table cannot be loaded
  2209.      * @since 6.1
  2210.      */
  2211.     public abstract TimeVectorFunction getEOPTidalCorrection()
  2212.         throws OrekitException;

  2214.     /** Get the Love numbers.
  2215.      * @return Love numbers
  2216.      * @exception OrekitException if table cannot be loaded
  2217.      * @since 6.1
  2218.      */
  2219.     public abstract LoveNumbers getLoveNumbers()
  2220.         throws OrekitException;

  2222.     /** Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2223.      * @param ut1 UT1 time scale
  2224.      * @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2225.      * @exception OrekitException if table cannot be loaded
  2226.      * @since 6.1
  2227.      */
  2228.     public abstract TimeVectorFunction getTideFrequencyDependenceFunction(TimeScale ut1)
  2229.         throws OrekitException;

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

  2237.     /** Get the function computing solid pole tide (ΔC₂₁, ΔS₂₁).
  2238.      * @param eopHistory EOP history
  2239.      * @return model for solid pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2240.      * @exception OrekitException if table cannot be loaded
  2241.      * @since 6.1
  2242.      */
  2243.     public abstract TimeVectorFunction getSolidPoleTide(EOPHistory eopHistory)
  2244.         throws OrekitException;

  2246.     /** Get the function computing ocean pole tide (ΔC₂₁, ΔS₂₁).
  2247.      * @param eopHistory EOP history
  2248.      * @return model for ocean pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  2249.      * @exception OrekitException if table cannot be loaded
  2250.      * @since 6.1
  2251.      */
  2252.     public abstract TimeVectorFunction getOceanPoleTide(EOPHistory eopHistory)
  2253.         throws OrekitException;

  2255.     /** Interface for functions converting nutation corrections between
  2256.      * δΔψ/δΔε to δX/δY.
  2257.      * <ul>
  2258.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  2259.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  2260.      * </ul>
  2261.      * @since 6.1
  2262.      */
  2263.     public interface NutationCorrectionConverter {

  2265.         /** Convert nutation corrections.
  2266.          * @param date current date
  2267.          * @param ddPsi δΔψ part of the nutation correction
  2268.          * @param ddEpsilon δΔε part of the nutation correction
  2269.          * @return array containing δX and δY
  2270.          * @exception OrekitException if correction cannot be converted
  2271.          */
  2272.         double[] toNonRotating(AbsoluteDate date, double ddPsi, double ddEpsilon)
  2273.             throws OrekitException;

  2275.         /** Convert nutation corrections.
  2276.          * @param date current date
  2277.          * @param dX δX part of the nutation correction
  2278.          * @param dY δY part of the nutation correction
  2279.          * @return array containing δΔψ and δΔε
  2280.          * @exception OrekitException if correction cannot be converted
  2281.          */
  2282.         double[] toEquinox(AbsoluteDate date, double dX, double dY)
  2283.             throws OrekitException;

  2285.     }

  2287.     /** Create a function converting nutation corrections between
  2288.      * δX/δY and δΔψ/δΔε.
  2289.      * <ul>
  2290.      * <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  2291.      * <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
  2292.      * </ul>
  2293.      * @return a new converter
  2294.      * @exception OrekitException if some convention table cannot be loaded
  2295.      * @since 6.1
  2296.      */
  2297.     public NutationCorrectionConverter getNutationCorrectionConverter()
  2298.         throws OrekitException {

  2300.         // get models parameters
  2301.         final TimeVectorFunction precessionFunction = getPrecessionFunction();
  2302.         final TimeScalarFunction epsilonAFunction = getMeanObliquityFunction();
  2303.         final AbsoluteDate date0 = getNutationReferenceEpoch();
  2304.         final double cosE0 = FastMath.cos(epsilonAFunction.value(date0));

  2306.         return new NutationCorrectionConverter() {

  2308.             /** {@inheritDoc} */
  2309.             @Override
  2310.             public double[] toNonRotating(final AbsoluteDate date,
  2311.                                           final double ddPsi, final double ddEpsilon)
  2312.                 throws OrekitException {
  2313.                 // compute precession angles psiA, omegaA and chiA
  2314.                 final double[] angles = precessionFunction.value(date);

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

  2320.                 // convert nutation corrections (equation 23/IERS-2003 or 5.25/IERS-2010)
  2321.                 return new double[] {
  2322.                     sinEA * ddPsi + c * ddEpsilon,
  2323.                     ddEpsilon - c * sinEA * ddPsi
  2324.                 };

  2326.             }

  2328.             /** {@inheritDoc} */
  2329.             @Override
  2330.             public double[] toEquinox(final AbsoluteDate date,
  2331.                                       final double dX, final double dY)
  2332.                 throws OrekitException {
  2333.                 // compute precession angles psiA, omegaA and chiA
  2334.                 final double[] angles   = precessionFunction.value(date);

  2336.                 // conversion coefficients
  2337.                 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
  2338.                 final double c     = angles[0] * cosE0 - angles[2];
  2339.                 final double opC2  = 1 + c * c;

  2341.                 // convert nutation corrections (inverse of equation 23/IERS-2003 or 5.25/IERS-2010)
  2342.                 return new double[] {
  2343.                     (dX - c * dY) / (sinEA * opC2),
  2344.                     (dY + c * dX) / opC2
  2345.                 };
  2346.             }

  2348.         };

  2350.     }

  2352.     /** Load the Love numbers.
  2353.      * @param nameLove name of the Love number resource
  2354.      * @return Love numbers
  2355.      * @exception OrekitException if the Love numbers embedded in the
  2356.      * library cannot be read
  2357.      */
  2358.     protected LoveNumbers loadLoveNumbers(final String nameLove) throws OrekitException {
  2359.         try {

  2361.             // allocate the three triangular arrays for real, imaginary and time-dependent numbers
  2362.             final double[][] real      = new double[4][];
  2363.             final double[][] imaginary = new double[4][];
  2364.             final double[][] plus      = new double[4][];
  2365.             for (int i = 0; i < real.length; ++i) {
  2366.                 real[i]      = new double[i + 1];
  2367.                 imaginary[i] = new double[i + 1];
  2368.                 plus[i]      = new double[i + 1];
  2369.             }

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

  2373.                 if (stream == null) {
  2374.                     // this should never happen with files embedded within Orekit
  2375.                     throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, nameLove);
  2376.                 }

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

  2381.                     String line = reader.readLine();
  2382.                     int lineNumber = 1;

  2384.                     // look for the Love numbers
  2385.                     while (line != null) {

  2387.                         line = line.trim();
  2388.                         if (!(line.isEmpty() || line.startsWith("#"))) {
  2389.                             final String[] fields = line.split("\\p{Space}+");
  2390.                             if (fields.length != 5) {
  2391.                                 // this should never happen with files embedded within Orekit
  2392.                                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  2393.                                                           lineNumber, nameLove, line);
  2394.                             }
  2395.                             final int n = Integer.parseInt(fields[0]);
  2396.                             final int m = Integer.parseInt(fields[1]);
  2397.                             if ((n < 2) || (n > 3) || (m < 0) || (m > n)) {
  2398.                                 // this should never happen with files embedded within Orekit
  2399.                                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  2400.                                                           lineNumber, nameLove, line);

  2402.                             }
  2403.                             real[n][m]      = Double.parseDouble(fields[2]);
  2404.                             imaginary[n][m] = Double.parseDouble(fields[3]);
  2405.                             plus[n][m]      = Double.parseDouble(fields[4]);
  2406.                         }

  2408.                         // next line
  2409.                         lineNumber++;
  2410.                         line = reader.readLine();

  2412.                     }
  2413.                 }
  2414.             }

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

  2418.         } catch (IOException ioe) {
  2419.             // this should never happen with files embedded within Orekit
  2420.             throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, nameLove);
  2421.         }
  2422.     }

  2424.     /** Get a data stream.
  2425.      * @param name file name of the resource stream
  2426.      * @return stream
  2427.      */
  2428.     private static InputStream getStream(final String name) {
  2429.         return IERSConventions.class.getResourceAsStream(name);
  2430.     }

  2432.     /** Correction to equation of equinoxes.
  2433.      * <p>IAU 1994 resolution C7 added two terms to the equation of equinoxes
  2434.      * taking effect since 1997-02-27 for continuity
  2435.      * </p>
  2436.      */
  2437.     private static class IAU1994ResolutionC7 {

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

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

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

  2451.         /** Evaluate the correction.
  2452.          * @param arguments Delaunay for nutation
  2453.          * @return correction value (0 before 1997-02-27)
  2454.          */
  2455.         public static double value(final DelaunayArguments arguments) {
  2456.             if (arguments.getDate().compareTo(MODEL_START) >= 0) {

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

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

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

  2467.             } else {
  2468.                 return 0.0;
  2469.             }
  2470.         }

  2472.         /** Evaluate the correction.
  2473.          * @param arguments Delaunay for nutation
  2474.          * @param <T> type of the field elements
  2475.          * @return correction value (0 before 1997-02-27)
  2476.          */
  2477.         public static <T extends RealFieldElement<T>> T value(final FieldDelaunayArguments<T> arguments) {
  2478.             if (arguments.getDate().toAbsoluteDate().compareTo(MODEL_START) >= 0) {

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

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

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

  2489.             } else {
  2490.                 return arguments.getDate().getField().getZero();
  2491.             }
  2492.         }

  2494.     };

  2496.     /** Stellar angle model.
  2497.      * <p>
  2498.      * The stellar angle computed here has been defined in the paper "A non-rotating origin on the
  2499.      * instantaneous equator: Definition, properties and use", N. Capitaine, Guinot B., and Souchay J.,
  2500.      * Celestial Mechanics, Volume 39, Issue 3, pp 283-307. It has been proposed as a conventional
  2501.      * conventional relationship between UT1 and Earth rotation in the paper "Definition of the Celestial
  2502.      * Ephemeris origin and of UT1 in the International Celestial Reference Frame", Capitaine, N.,
  2503.      * Guinot, B., and McCarthy, D. D., 2000, Astronomy and Astrophysics, 355(1), pp. 398–405.
  2504.      * </p>
  2505.      * <p>
  2506.      * It is presented simply as stellar angle in IERS conventions 1996 but as since been adopted as
  2507.      * the conventional relationship defining UT1 from ICRF and is called Earth Rotation Angle in
  2508.      * IERS conventions 2003 and 2010.
  2509.      * </p>
  2510.      */
  2511.     private static class StellarAngleCapitaine implements TimeScalarFunction {

  2513.         /** Reference date of Capitaine's Earth Rotation Angle model. */
  2514.         private static final AbsoluteDate REFERENCE_DATE = new AbsoluteDate(DateComponents.J2000_EPOCH,
  2515.                                                                             TimeComponents.H12,
  2516.                                                                             TimeScalesFactory.getTAI());

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

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

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

  2529.         /** UT1 time scale. */
  2530.         private final TimeScale ut1;

  2532.         /** Simple constructor.
  2533.          * @param ut1 UT1 time scale
  2534.          */
  2535.         StellarAngleCapitaine(final TimeScale ut1) {
  2536.             this.ut1 = ut1;
  2537.         }

  2539.         /** Get the rotation rate.
  2540.          * @return rotation rate
  2541.          */
  2542.         public double getRate() {
  2543.             return ERA_1A + ERA_1B;
  2544.         }

  2546.         /** {@inheritDoc} */
  2547.         @Override
  2548.         public double value(final AbsoluteDate date) {

  2550.             // split the date offset as a full number of days plus a smaller part
  2551.             final int secondsInDay = 86400;
  2552.             final double dt  = date.durationFrom(REFERENCE_DATE);
  2553.             final long days  = ((long) dt) / secondsInDay;
  2554.             final double dtA = secondsInDay * days;
  2555.             final double dtB = (dt - dtA) + ut1.offsetFromTAI(date);

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

  2559.         }

  2561.         /** {@inheritDoc} */
  2562.         @Override
  2563.         public <T extends RealFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {

  2565.             // split the date offset as a full number of days plus a smaller part
  2566.             final int secondsInDay = 86400;
  2567.             final T dt  = date.durationFrom(REFERENCE_DATE);
  2568.             final long days  = ((long) dt.getReal()) / secondsInDay;
  2569.             final double dtA = secondsInDay * days;
  2570.             final T dtB = dt.subtract(dtA).add(ut1.offsetFromTAI(date.toAbsoluteDate()));

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

  2574.         }

  2576.     }

  2578.     /** Mean pole. */
  2579.     private static class MeanPole implements TimeStamped, Serializable {

  2581.         /** Serializable UID. */
  2582.         private static final long serialVersionUID = 20131028l;

  2584.         /** Date. */
  2585.         private final AbsoluteDate date;

  2587.         /** X coordinate. */
  2588.         private double x;

  2590.         /** Y coordinate. */
  2591.         private double y;

  2593.         /** Simple constructor.
  2594.          * @param date date
  2595.          * @param x x coordinate
  2596.          * @param y y coordinate
  2597.          */
  2598.         MeanPole(final AbsoluteDate date, final double x, final double y) {
  2599.             this.date = date;
  2600.             this.x    = x;
  2601.             this.y    = y;
  2602.         }

  2604.         /** {@inheritDoc} */
  2605.         @Override
  2606.         public AbsoluteDate getDate() {
  2607.             return date;
  2608.         }

  2610.         /** Get x coordinate.
  2611.          * @return x coordinate
  2612.          */
  2613.         public double getX() {
  2614.             return x;
  2615.         }

  2617.         /** Get y coordinate.
  2618.          * @return y coordinate
  2619.          */
  2620.         public double getY() {
  2621.             return y;
  2622.         }

  2624.     }

  2626.     /** Local class for precession function. */
  2627.     private class PrecessionFunction implements TimeVectorFunction {

  2629.         /** Polynomial nutation for psiA. */
  2630.         private final PolynomialNutation psiA;

  2632.         /** Polynomial nutation for omegaA. */
  2633.         private final PolynomialNutation omegaA;

  2635.         /** Polynomial nutation for chiA. */
  2636.         private final PolynomialNutation chiA;

  2638.         /** Simple constructor.
  2639.          * @param psiA polynomial nutation for psiA
  2640.          * @param omegaA polynomial nutation for omegaA
  2641.          * @param chiA polynomial nutation for chiA
  2642.          */
  2643.         PrecessionFunction(final PolynomialNutation psiA,
  2644.                            final PolynomialNutation omegaA,
  2645.                            final PolynomialNutation chiA) {
  2646.             this.psiA   = psiA;
  2647.             this.omegaA = omegaA;
  2648.             this.chiA   = chiA;
  2649.         }


  2652.         /** {@inheritDoc} */
  2653.         @Override
  2654.         public double[] value(final AbsoluteDate date) {
  2655.             final double tc = evaluateTC(date);
  2656.             return new double[] {
  2657.                 psiA.value(tc), omegaA.value(tc), chiA.value(tc)
  2658.             };
  2659.         }

  2661.         /** {@inheritDoc} */
  2662.         @Override
  2663.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  2664.             final T[] a = MathArrays.buildArray(date.getField(), 3);
  2665.             final T tc = evaluateTC(date);
  2666.             a[0] = psiA.value(tc);
  2667.             a[1] = omegaA.value(tc);
  2668.             a[2] = chiA.value(tc);
  2669.             return a;
  2670.         }

  2672.     }

  2674.     /** Local class for tides frequency function. */
  2675.     private static class TideFrequencyDependenceFunction implements TimeVectorFunction {

  2677.         /** Nutation arguments. */
  2678.         private final FundamentalNutationArguments arguments;

  2680.         /** Correction series. */
  2681.         private final PoissonSeries.CompiledSeries kSeries;

  2683.         /** Simple constructor.
  2684.          * @param arguments nutation arguments
  2685.          * @param c20Series correction series for the C20 term
  2686.          * @param c21Series correction series for the C21 term
  2687.          * @param s21Series correction series for the S21 term
  2688.          * @param c22Series correction series for the C22 term
  2689.          * @param s22Series correction series for the S22 term
  2690.          */
  2691.         TideFrequencyDependenceFunction(final FundamentalNutationArguments arguments,
  2692.                                         final PoissonSeries c20Series,
  2693.                                         final PoissonSeries c21Series,
  2694.                                         final PoissonSeries s21Series,
  2695.                                         final PoissonSeries c22Series,
  2696.                                         final PoissonSeries s22Series) {
  2697.             this.arguments = arguments;
  2698.             this.kSeries   = PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);
  2699.         }


  2702.         /** {@inheritDoc} */
  2703.         @Override
  2704.         public double[] value(final AbsoluteDate date) {
  2705.             return kSeries.value(arguments.evaluateAll(date));
  2706.         }

  2708.         /** {@inheritDoc} */
  2709.         @Override
  2710.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  2711.             return kSeries.value(arguments.evaluateAll(date));
  2712.         }

  2714.     }

  2716.     /** Local class for EOP tidal corrections. */
  2717.     private static class EOPTidalCorrection implements TimeVectorFunction {

  2719.         /** Nutation arguments. */
  2720.         private final FundamentalNutationArguments arguments;

  2722.         /** Correction series. */
  2723.         private final PoissonSeries.CompiledSeries correctionSeries;

  2725.         /** Simple constructor.
  2726.          * @param arguments nutation arguments
  2727.          * @param xSeries correction series for the x coordinate
  2728.          * @param ySeries correction series for the y coordinate
  2729.          * @param ut1Series correction series for the UT1
  2730.          */
  2731.         EOPTidalCorrection(final FundamentalNutationArguments arguments,
  2732.                            final PoissonSeries xSeries,
  2733.                            final PoissonSeries ySeries,
  2734.                            final PoissonSeries ut1Series) {
  2735.             this.arguments        = arguments;
  2736.             this.correctionSeries = PoissonSeries.compile(xSeries, ySeries, ut1Series);
  2737.         }

  2739.         /** {@inheritDoc} */
  2740.         @Override
  2741.         public double[] value(final AbsoluteDate date) {
  2742.             final BodiesElements elements = arguments.evaluateAll(date);
  2743.             final double[] correction    = correctionSeries.value(elements);
  2744.             final double[] correctionDot = correctionSeries.derivative(elements);
  2745.             return new double[] {
  2746.                 correction[0],
  2747.                 correction[1],
  2748.                 correction[2],
  2749.                 correctionDot[2] * (-Constants.JULIAN_DAY)
  2750.             };
  2751.         }

  2753.         /** {@inheritDoc} */
  2754.         @Override
  2755.         public <T extends RealFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {

  2757.             final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
  2758.             final T[] correction    = correctionSeries.value(elements);
  2759.             final T[] correctionDot = correctionSeries.derivative(elements);
  2760.             final T[] a = MathArrays.buildArray(date.getField(), 4);
  2761.             a[0] = correction[0];
  2762.             a[1] = correction[1];
  2763.             a[2] = correction[2];
  2764.             a[3] = correctionDot[2].multiply(-Constants.JULIAN_DAY);
  2765.             return a;
  2766.         }

  2768.     }

  2770. }
Created with JaCoCo 0.7.9.201702052155