IERSConventions.java

  1. /* Copyright 2002-2013 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.apache.commons.math3.analysis.differentiation.DerivativeStructure;
  27. import org.apache.commons.math3.analysis.interpolation.HermiteInterpolator;
  28. import org.apache.commons.math3.util.FastMath;
  29. import org.apache.commons.math3.util.MathUtils;
  30. import org.orekit.data.BodiesElements;
  31. import org.orekit.data.DelaunayArguments;
  32. import org.orekit.data.FieldBodiesElements;
  33. import org.orekit.data.FundamentalNutationArguments;
  34. import org.orekit.data.PoissonSeries;
  35. import org.orekit.data.PoissonSeriesParser;
  36. import org.orekit.data.PolynomialNutation;
  37. import org.orekit.data.PolynomialParser;
  38. import org.orekit.data.PolynomialParser.Unit;
  39. import org.orekit.data.SimpleTimeStampedTableParser;
  40. import org.orekit.errors.OrekitException;
  41. import org.orekit.errors.OrekitMessages;
  42. import org.orekit.errors.TimeStampedCacheException;
  43. import org.orekit.frames.EOPHistory;
  44. import org.orekit.frames.PoleCorrection;
  45. import org.orekit.time.AbsoluteDate;
  46. import org.orekit.time.DateComponents;
  47. import org.orekit.time.TimeComponents;
  48. import org.orekit.time.TimeFunction;
  49. import org.orekit.time.TimeScale;
  50. import org.orekit.time.TimeScalesFactory;
  51. import org.orekit.time.TimeStamped;


  54. /** Supported IERS conventions.
  55.  * @since 6.0
  56.  * @author Luc Maisonobe
  57.  */
  58. public enum IERSConventions {

  60.     /** Constant for IERS 1996 conventions. */
  61.     IERS_1996 {

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

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

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

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

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

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

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

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

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

  90.         /** {@inheritDoc} */
  91.         @Override
  92.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
  93.             throws OrekitException {
  94.             return new FundamentalNutationArguments(this, timeScale,
  95.                                                     getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
  96.         }

  98.         /** {@inheritDoc} */
  99.         @Override
  100.         public TimeFunction<Double> getMeanObliquityFunction() throws OrekitException {

  102.             // value from chapter 5, page 22
  103.             final PolynomialNutation<DerivativeStructure> epsilonA =
  104.                     new PolynomialNutation<DerivativeStructure>(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  105.                                                                   -46.8150   * Constants.ARC_SECONDS_TO_RADIANS,
  106.                                                                    -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  107.                                                                     0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  109.             return new TimeFunction<Double>() {

  111.                 /** {@inheritDoc} */
  112.                 @Override
  113.                 public Double value(final AbsoluteDate date) {
  114.                     return epsilonA.value(evaluateTC(date));
  115.                 }

  117.             };

  119.         }

  121.         /** {@inheritDoc} */
  122.         @Override
  123.         public TimeFunction<double[]> getXYSpXY2Function()
  124.             throws OrekitException {

  126.             // set up nutation arguments
  127.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  129.             // X = 2004.3109″t - 0.42665″t² - 0.198656″t³ + 0.0000140″t⁴
  130.             //     + 0.00006″t² cos Ω + sin ε0 { Σ [(Ai + Ai' t) sin(ARGUMENT) + Ai'' t cos(ARGUMENT)]}
  131.             //     + 0.00204″t² sin Ω + 0.00016″t² sin 2(F - D + Ω),
  132.             final PolynomialNutation<DerivativeStructure> xPolynomial =
  133.                     new PolynomialNutation<DerivativeStructure>(0,
  134.                                                                 2004.3109 * Constants.ARC_SECONDS_TO_RADIANS,
  135.                                                                 -0.42665  * Constants.ARC_SECONDS_TO_RADIANS,
  136.                                                                 -0.198656 * Constants.ARC_SECONDS_TO_RADIANS,
  137.                                                                 0.0000140 * Constants.ARC_SECONDS_TO_RADIANS);

  139.             final double fXCosOm    = 0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
  140.             final double fXSinOm    = 0.00204 * Constants.ARC_SECONDS_TO_RADIANS;
  141.             final double fXSin2FDOm = 0.00016 * Constants.ARC_SECONDS_TO_RADIANS;
  142.             final double sinEps0   = FastMath.sin(getMeanObliquityFunction().value(getNutationReferenceEpoch()));

  144.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  145.             final PoissonSeriesParser<DerivativeStructure> baseParser =
  146.                     new PoissonSeriesParser<DerivativeStructure>(12).withFirstDelaunay(1);

  148.             final PoissonSeriesParser<DerivativeStructure> xParser =
  149.                     baseParser.
  150.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  151.                     withSinCos(1, 8, deciMilliAS,  9, deciMilliAS);
  152.             final PoissonSeries<DerivativeStructure> xSum = xParser.parse(getStream(X_Y_SERIES), X_Y_SERIES);

  154.             // Y = -0.00013″ - 22.40992″t² + 0.001836″t³ + 0.0011130″t⁴
  155.             //     + Σ [(Bi + Bi' t) cos(ARGUMENT) + Bi'' t sin(ARGUMENT)]
  156.             //    - 0.00231″t² cos Ω − 0.00014″t² cos 2(F - D + Ω)
  157.             final PolynomialNutation<DerivativeStructure> yPolynomial =
  158.                     new PolynomialNutation<DerivativeStructure>(-0.00013  * Constants.ARC_SECONDS_TO_RADIANS,
  159.                                                                 0.0,
  160.                                                                 -22.40992 * Constants.ARC_SECONDS_TO_RADIANS,
  161.                                                                 0.001836  * Constants.ARC_SECONDS_TO_RADIANS,
  162.                                                                 0.0011130 * Constants.ARC_SECONDS_TO_RADIANS);

  164.             final double fYCosOm    = -0.00231 * Constants.ARC_SECONDS_TO_RADIANS;
  165.             final double fYCos2FDOm = -0.00014 * Constants.ARC_SECONDS_TO_RADIANS;

  167.             final PoissonSeriesParser<DerivativeStructure> yParser =
  168.                     baseParser.
  169.                     withSinCos(0, -1, deciMilliAS, 10, deciMilliAS).
  170.                     withSinCos(1, 12, deciMilliAS, 11, deciMilliAS);
  171.             final PoissonSeries<DerivativeStructure> ySum = yParser.parse(getStream(X_Y_SERIES), X_Y_SERIES);

  173.             @SuppressWarnings("unchecked")
  174.             final PoissonSeries.CompiledSeries<DerivativeStructure> xySum =
  175.                     PoissonSeries.compile(xSum, ySum);

  177.             // s = -XY/2 + 0.00385″t - 0.07259″t³ - 0.00264″ sin Ω - 0.00006″ sin 2Ω
  178.             //     + 0.00074″t² sin Ω + 0.00006″t² sin 2(F - D + Ω)
  179.             final double fST          =  0.00385 * Constants.ARC_SECONDS_TO_RADIANS;
  180.             final double fST3         = -0.07259 * Constants.ARC_SECONDS_TO_RADIANS;
  181.             final double fSSinOm      = -0.00264 * Constants.ARC_SECONDS_TO_RADIANS;
  182.             final double fSSin2Om     = -0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
  183.             final double fST2SinOm    =  0.00074 * Constants.ARC_SECONDS_TO_RADIANS;
  184.             final double fST2Sin2FDOm =  0.00006 * Constants.ARC_SECONDS_TO_RADIANS;

  186.             return new TimeFunction<double[]>() {

  188.                 /** {@inheritDoc} */
  189.                 @Override
  190.                 public double[] value(final AbsoluteDate date) {

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

  195.                     final double omega     = elements.getOmega();
  196.                     final double f         = elements.getF();
  197.                     final double d         = elements.getD();
  198.                     final double t         = elements.getTC();

  200.                     final double cosOmega  = FastMath.cos(omega);
  201.                     final double sinOmega  = FastMath.sin(omega);
  202.                     final double sin2Omega = FastMath.sin(2 * omega);
  203.                     final double cos2FDOm  = FastMath.cos(2 * (f - d + omega));
  204.                     final double sin2FDOm  = FastMath.sin(2 * (f - d + omega));

  206.                     final double x = xPolynomial.value(t) + sinEps0 * xy[0] +
  207.                             t * t * (fXCosOm * cosOmega + fXSinOm * sinOmega + fXSin2FDOm * cos2FDOm);
  208.                     final double y = yPolynomial.value(t) + xy[1] +
  209.                             t * t * (fYCosOm * cosOmega + fYCos2FDOm * cos2FDOm);
  210.                     final double sPxy2 = fSSinOm * sinOmega + fSSin2Om * sin2Omega +
  211.                             t * (fST + t * (fST2SinOm * sinOmega + fST2Sin2FDOm * sin2FDOm + t * fST3));

  213.                     return new double[] {
  214.                         x, y, sPxy2
  215.                     };

  217.                 }

  219.             };

  221.         }

  223.         /** {@inheritDoc} */
  224.         @Override
  225.         public TimeFunction<double[]> getPrecessionFunction() throws OrekitException {

  227.             // set up the conventional polynomials
  228.             // the following values are from Lieske et al. paper:
  229.             // Expressions for the precession quantities based upon the IAU(1976) system of astronomical constants
  230.             // http://articles.adsabs.harvard.edu/full/1977A%26A....58....1L
  231.             // also available as equation 30 in IERS 2003 conventions
  232.             final PolynomialNutation<DerivativeStructure> psiA =
  233.                     new PolynomialNutation<DerivativeStructure>(   0.0,
  234.                                                                 5038.7784   * Constants.ARC_SECONDS_TO_RADIANS,
  235.                                                                   -1.07259  * Constants.ARC_SECONDS_TO_RADIANS,
  236.                                                                   -0.001147 * Constants.ARC_SECONDS_TO_RADIANS);
  237.             final PolynomialNutation<DerivativeStructure> omegaA =
  238.                     new PolynomialNutation<DerivativeStructure>(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
  239.                                                                  0.0,
  240.                                                                  0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
  241.                                                                 -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
  242.             final PolynomialNutation<DerivativeStructure> chiA =
  243.                     new PolynomialNutation<DerivativeStructure>( 0.0,
  244.                                                                 10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
  245.                                                                 -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
  246.                                                                 -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);

  248.             return new TimeFunction<double[]>() {
  249.                 /** {@inheritDoc} */
  250.                 @Override
  251.                 public double[] value(final AbsoluteDate date) {
  252.                     final double tc = evaluateTC(date);
  253.                     return new double[] {
  254.                         psiA.value(tc), omegaA.value(tc), chiA.value(tc)
  255.                     };
  256.                 }
  257.             };

  259.         }

  261.         /** {@inheritDoc} */
  262.         @Override
  263.         public TimeFunction<double[]> getNutationFunction()
  264.             throws OrekitException {

  266.             // set up nutation arguments
  267.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  269.             // set up Poisson series
  270.             final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
  271.             final PoissonSeriesParser<DerivativeStructure> baseParser =
  272.                     new PoissonSeriesParser<DerivativeStructure>(10).withFirstDelaunay(1);

  274.             final PoissonSeriesParser<DerivativeStructure> psiParser =
  275.                     baseParser.
  276.                     withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
  277.                     withSinCos(1, 8, deciMilliAS, -1, deciMilliAS);
  278.             final PoissonSeries<DerivativeStructure> psiSeries = psiParser.parse(getStream(PSI_EPSILON_SERIES), PSI_EPSILON_SERIES);

  280.             final PoissonSeriesParser<DerivativeStructure> epsilonParser =
  281.                     baseParser.
  282.                     withSinCos(0, -1, deciMilliAS, 9, deciMilliAS).
  283.                     withSinCos(1, -1, deciMilliAS, 10, deciMilliAS);
  284.             final PoissonSeries<DerivativeStructure> epsilonSeries = epsilonParser.parse(getStream(PSI_EPSILON_SERIES), PSI_EPSILON_SERIES);

  286.             @SuppressWarnings("unchecked")
  287.             final PoissonSeries.CompiledSeries<DerivativeStructure> psiEpsilonSeries =
  288.                     PoissonSeries.compile(psiSeries, epsilonSeries);

  290.             return new TimeFunction<double[]>() {
  291.                 /** {@inheritDoc} */
  292.                 @Override
  293.                 public double[] value(final AbsoluteDate date) {
  294.                     final BodiesElements elements = arguments.evaluateAll(date);
  295.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  296.                     return new double[] {
  297.                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements)
  298.                     };
  299.                 }
  300.             };

  302.         }

  304.         /** {@inheritDoc} */
  305.         @Override
  306.         public TimeFunction<DerivativeStructure> getGMSTFunction(final TimeScale ut1)
  307.             throws OrekitException {

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

  312.             // constants from IERS 1996 page 21
  313.             // the underlying model is IAU 1982 GMST-UT1
  314.             final AbsoluteDate gmstReference =
  315.                 new AbsoluteDate(DateComponents.J2000_EPOCH, TimeComponents.H12, TimeScalesFactory.getTAI());
  316.             final double gmst0 = 24110.54841;
  317.             final double gmst1 = 8640184.812866;
  318.             final double gmst2 = 0.093104;
  319.             final double gmst3 = -6.2e-6;

  321.             return new TimeFunction<DerivativeStructure>() {

  323.                 /** {@inheritDoc} */
  324.                 @Override
  325.                 public DerivativeStructure value(final AbsoluteDate date) {

  327.                     // offset in Julian centuries from J2000 epoch (UT1 scale)
  328.                     final double dtai = date.durationFrom(gmstReference);
  329.                     final DerivativeStructure tut1 =
  330.                             new DerivativeStructure(1, 1, dtai + ut1.offsetFromTAI(date), 1.0);
  331.                     final DerivativeStructure tt = tut1.divide(Constants.JULIAN_CENTURY);

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

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

  340.                 }

  342.             };

  344.         }

  346.         /** {@inheritDoc} */
  347.         @Override
  348.         public TimeFunction<DerivativeStructure> getGASTFunction(final TimeScale ut1,
  349.                                                                  final EOPHistory eopHistory)
  350.             throws OrekitException {

  352.             // obliquity
  353.             final TimeFunction<Double> epsilonA = getMeanObliquityFunction();

  355.             // GMST function
  356.             final TimeFunction<DerivativeStructure> gmst = getGMSTFunction(ut1);

  358.             // nutation function
  359.             final TimeFunction<double[]> nutation = getNutationFunction();

  361.             return new TimeFunction<DerivativeStructure>() {

  363.                 /** {@inheritDoc} */
  364.                 @Override
  365.                 public DerivativeStructure value(final AbsoluteDate date) {

  367.                     // compute equation of equinoxes
  368.                     final double[] angles = nutation.value(date);
  369.                     double deltaPsi = angles[0];
  370.                     if (eopHistory != null) {
  371.                         deltaPsi += eopHistory.getEquinoxNutationCorrection(date)[0];
  372.                     }
  373.                     final double eqe = deltaPsi  * FastMath.cos(epsilonA.value(date)) + angles[2];

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

  378.                 }

  380.             };

  382.         }

  384.         /** {@inheritDoc} */
  385.         @Override
  386.         public TimeFunction<double[]> getEOPTidalCorrection()
  387.             throws OrekitException {

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

  396.             // set up Poisson series
  397.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  398.             final PoissonSeriesParser<DerivativeStructure> xyParser = new PoissonSeriesParser<DerivativeStructure>(17).
  399.                     withOptionalColumn(1).
  400.                     withGamma(7).
  401.                     withFirstDelaunay(2);
  402.             final PoissonSeries<DerivativeStructure> xSeries =
  403.                     xyParser.
  404.                     withSinCos(0, 14, milliAS, 15, milliAS).
  405.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
  406.             final PoissonSeries<DerivativeStructure> ySeries =
  407.                     xyParser.
  408.                     withSinCos(0, 16, milliAS, 17, milliAS).
  409.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES),
  410.                                                          TIDAL_CORRECTION_XP_YP_SERIES);

  412.             final double deciMilliS = 1.0e-4;
  413.             final PoissonSeriesParser<DerivativeStructure> ut1Parser = new PoissonSeriesParser<DerivativeStructure>(17).
  414.                     withOptionalColumn(1).
  415.                     withGamma(7).
  416.                     withFirstDelaunay(2).
  417.                     withSinCos(0, 16, deciMilliS, 17, deciMilliS);
  418.             final PoissonSeries<DerivativeStructure> ut1Series =
  419.                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);

  421.             @SuppressWarnings("unchecked")
  422.             final PoissonSeries.CompiledSeries<DerivativeStructure> correctionSeries =
  423.                 PoissonSeries.compile(xSeries, ySeries, ut1Series);

  425.             return new TimeFunction<double[]>() {
  426.                 /** {@inheritDoc} */
  427.                 @Override
  428.                 public double[] value(final AbsoluteDate date) {
  429.                     final FieldBodiesElements<DerivativeStructure> elements =
  430.                             arguments.evaluateDerivative(date);
  431.                     final DerivativeStructure[] correction = correctionSeries.value(elements);
  432.                     return new double[] {
  433.                         correction[0].getValue(),
  434.                         correction[1].getValue(),
  435.                         correction[2].getValue(),
  436.                         -correction[2].getPartialDerivative(1) * Constants.JULIAN_DAY
  437.                     };
  438.                 }
  439.             };

  441.         }

  443.         /** {@inheritDoc} */
  444.         public LoveNumbers getLoveNumbers() throws OrekitException {
  445.             return loadLoveNumbers(LOVE_NUMBERS);
  446.         }

  448.         /** {@inheritDoc} */
  449.         public TimeFunction<double[]> getTideFrequencyDependenceFunction(final TimeScale ut1)
  450.             throws OrekitException {

  452.             // set up nutation arguments
  453.             final FundamentalNutationArguments arguments = getNutationArguments(ut1);

  455.             // set up Poisson series
  456.             final PoissonSeriesParser<DerivativeStructure> k20Parser =
  457.                     new PoissonSeriesParser<DerivativeStructure>(18).
  458.                         withOptionalColumn(1).
  459.                         withDoodson(4, 2).
  460.                         withFirstDelaunay(10);
  461.             final PoissonSeriesParser<DerivativeStructure> k21Parser =
  462.                     new PoissonSeriesParser<DerivativeStructure>(18).
  463.                         withOptionalColumn(1).
  464.                         withDoodson(4, 3).
  465.                         withFirstDelaunay(10);
  466.             final PoissonSeriesParser<DerivativeStructure> k22Parser =
  467.                     new PoissonSeriesParser<DerivativeStructure>(16).
  468.                         withOptionalColumn(1).
  469.                         withDoodson(4, 2).
  470.                         withFirstDelaunay(10);

  472.             final double pico = 1.0e-12;
  473.             final PoissonSeries<DerivativeStructure> c20Series =
  474.                     k20Parser.
  475.                     // TODO: check sign!
  476. //                  withSinCos(0, 18, pico, 16, pico).
  477.                   withSinCos(0, 18, -pico, 16, pico).
  478.                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
  479.             final PoissonSeries<DerivativeStructure> c21Series =
  480.                     k21Parser.
  481.                     withSinCos(0, 17, pico, 18, pico).
  482.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  483.             final PoissonSeries<DerivativeStructure> s21Series =
  484.                     k21Parser.
  485.                     withSinCos(0, 18, -pico, 17, pico).
  486.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  487.             final PoissonSeries<DerivativeStructure> c22Series =
  488.                     k22Parser.
  489.                     withSinCos(0, -1, pico, 16, pico).
  490.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
  491.             final PoissonSeries<DerivativeStructure> s22Series =
  492.                     k22Parser.
  493.                     withSinCos(0, 16, -pico, -1, pico).
  494.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);

  496.             @SuppressWarnings("unchecked")
  497.             final PoissonSeries.CompiledSeries<DerivativeStructure> kSeries =
  498.                 PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);

  500.             return new TimeFunction<double[]>() {
  501.                 /** {@inheritDoc} */
  502.                 @Override
  503.                 public double[] value(final AbsoluteDate date) {
  504.                     return kSeries.value(arguments.evaluateAll(date));
  505.                 }
  506.             };

  508.         }

  510.         /** {@inheritDoc} */
  511.         @Override
  512.         public double getPermanentTide() throws OrekitException {
  513.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  514.         }

  516.         /** {@inheritDoc} */
  517.         @Override
  518.         public TimeFunction<double[]> getSolidPoleTide(final EOPHistory eopHistory) {

  520.             // constants from IERS 1996 page 47
  521.             final double globalFactor = -1.348e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  522.             final double coupling     =  0.00112;

  524.             return new TimeFunction<double[]>() {
  525.                 /** {@inheritDoc} */
  526.                 @Override
  527.                 public double[] value(final AbsoluteDate date) {
  528.                     final PoleCorrection pole = eopHistory.getPoleCorrection(date);
  529.                     return new double[] {
  530.                         globalFactor * (pole.getXp() + coupling * pole.getYp()),
  531.                         globalFactor * (coupling * pole.getXp() - pole.getYp()),
  532.                     };
  533.                 }
  534.             };
  535.         }

  537.         /** {@inheritDoc} */
  538.         @Override
  539.         public TimeFunction<double[]> getOceanPoleTide(final EOPHistory eopHistory)
  540.             throws OrekitException {

  542.             return new TimeFunction<double[]>() {
  543.                 /** {@inheritDoc} */
  544.                 @Override
  545.                 public double[] value(final AbsoluteDate date) {
  546.                     // there are no model for ocean pole tide prior to conventions 2010
  547.                     return new double[] {
  548.                         0.0, 0.0
  549.                     };
  550.                 }
  551.             };
  552.         }

  554.     },

  556.     /** Constant for IERS 2003 conventions. */
  557.     IERS_2003 {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  601.         /** {@inheritDoc} */
  602.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
  603.             throws OrekitException {
  604.             return new FundamentalNutationArguments(this, timeScale,
  605.                                                     getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
  606.         }

  608.         /** {@inheritDoc} */
  609.         @Override
  610.         public TimeFunction<Double> getMeanObliquityFunction() throws OrekitException {

  612.             // epsilon 0 value from chapter 5, page 41, other terms from equation 32 page 45
  613.             final PolynomialNutation<DerivativeStructure> epsilonA =
  614.                     new PolynomialNutation<DerivativeStructure>(84381.448    * Constants.ARC_SECONDS_TO_RADIANS,
  615.                                                                   -46.84024  * Constants.ARC_SECONDS_TO_RADIANS,
  616.                                                                    -0.00059  * Constants.ARC_SECONDS_TO_RADIANS,
  617.                                                                     0.001813 * Constants.ARC_SECONDS_TO_RADIANS);

  619.             return new TimeFunction<Double>() {

  621.                 /** {@inheritDoc} */
  622.                 @Override
  623.                 public Double value(final AbsoluteDate date) {
  624.                     return epsilonA.value(evaluateTC(date));
  625.                 }

  627.             };

  629.         }

  631.         /** {@inheritDoc} */
  632.         @Override
  633.         public TimeFunction<double[]> getXYSpXY2Function()
  634.             throws OrekitException {

  636.             // set up nutation arguments
  637.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  639.             // set up Poisson series
  640.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  641.             final PoissonSeriesParser<DerivativeStructure> parser =
  642.                     new PoissonSeriesParser<DerivativeStructure>(17).
  643.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  644.                         withFirstDelaunay(4).
  645.                         withFirstPlanetary(9).
  646.                         withSinCos(0, 2, microAS, 3, microAS);

  648.             final PoissonSeries<DerivativeStructure> xSeries = parser.parse(getStream(X_SERIES), X_SERIES);
  649.             final PoissonSeries<DerivativeStructure> ySeries = parser.parse(getStream(Y_SERIES), Y_SERIES);
  650.             final PoissonSeries<DerivativeStructure> sSeries = parser.parse(getStream(S_SERIES), S_SERIES);
  651.             @SuppressWarnings("unchecked")
  652.             final PoissonSeries.CompiledSeries<DerivativeStructure> xys = PoissonSeries.compile(xSeries, ySeries, sSeries);

  654.             // create a function evaluating the series
  655.             return new TimeFunction<double[]>() {

  657.                 /** {@inheritDoc} */
  658.                 @Override
  659.                 public double[] value(final AbsoluteDate date) {
  660.                     return xys.value(arguments.evaluateAll(date));
  661.                 }

  663.             };

  665.         }


  668.         /** {@inheritDoc} */
  669.         @Override
  670.         public TimeFunction<double[]> getPrecessionFunction() throws OrekitException {

  672.             // set up the conventional polynomials
  673.             // the following values are from equation 32 in IERS 2003 conventions
  674.             final PolynomialNutation<DerivativeStructure> psiA =
  675.                     new PolynomialNutation<DerivativeStructure>(    0.0,
  676.                                                                  5038.47875   * Constants.ARC_SECONDS_TO_RADIANS,
  677.                                                                    -1.07259   * Constants.ARC_SECONDS_TO_RADIANS,
  678.                                                                    -0.001147  * Constants.ARC_SECONDS_TO_RADIANS);
  679.             final PolynomialNutation<DerivativeStructure> omegaA =
  680.                     new PolynomialNutation<DerivativeStructure>(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
  681.                                                                 -0.02524   * Constants.ARC_SECONDS_TO_RADIANS,
  682.                                                                  0.05127   * Constants.ARC_SECONDS_TO_RADIANS,
  683.                                                                 -0.007726  * Constants.ARC_SECONDS_TO_RADIANS);
  684.             final PolynomialNutation<DerivativeStructure> chiA =
  685.                     new PolynomialNutation<DerivativeStructure>( 0.0,
  686.                                                                 10.5526   * Constants.ARC_SECONDS_TO_RADIANS,
  687.                                                                 -2.38064  * Constants.ARC_SECONDS_TO_RADIANS,
  688.                                                                 -0.001125 * Constants.ARC_SECONDS_TO_RADIANS);

  690.             return new TimeFunction<double[]>() {
  691.                 /** {@inheritDoc} */
  692.                 @Override
  693.                 public double[] value(final AbsoluteDate date) {
  694.                     final double tc = evaluateTC(date);
  695.                     return new double[] {
  696.                         psiA.value(tc), omegaA.value(tc), chiA.value(tc)
  697.                     };
  698.                 }
  699.             };

  701.         }

  703.         /** {@inheritDoc} */
  704.         @Override
  705.         public TimeFunction<double[]> getNutationFunction()
  706.             throws OrekitException {

  708.             // set up nutation arguments
  709.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  711.             // set up Poisson series
  712.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  713.             final PoissonSeriesParser<DerivativeStructure> luniSolarParser =
  714.                     new PoissonSeriesParser<DerivativeStructure>(14).withFirstDelaunay(1);
  715.             final PoissonSeriesParser<DerivativeStructure> luniSolarPsiParser =
  716.                     luniSolarParser.
  717.                     withSinCos(0, 7, milliAS, 11, milliAS).
  718.                     withSinCos(1, 8, milliAS, 12, milliAS);
  719.             final PoissonSeries<DerivativeStructure> psiLuniSolarSeries =
  720.                     luniSolarPsiParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);
  721.             final PoissonSeriesParser<DerivativeStructure> luniSolarEpsilonParser =
  722.                     luniSolarParser.
  723.                     withSinCos(0, 13, milliAS, 9, milliAS).
  724.                     withSinCos(1, 14, milliAS, 10, milliAS);
  725.             final PoissonSeries<DerivativeStructure> epsilonLuniSolarSeries =
  726.                     luniSolarEpsilonParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);

  728.             final PoissonSeriesParser<DerivativeStructure> planetaryParser =
  729.                     new PoissonSeriesParser<DerivativeStructure>(21).
  730.                         withFirstDelaunay(2).
  731.                         withFirstPlanetary(7);
  732.             final PoissonSeriesParser<DerivativeStructure> planetaryPsiParser =
  733.                     planetaryParser.withSinCos(0, 17, milliAS, 18, milliAS);
  734.             final PoissonSeries<DerivativeStructure> psiPlanetarySeries =
  735.                     planetaryPsiParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);
  736.             final PoissonSeriesParser<DerivativeStructure> planetaryEpsilonParser =
  737.                     planetaryParser.withSinCos(0, 19, milliAS, 20, milliAS);
  738.             final PoissonSeries<DerivativeStructure> epsilonPlanetarySeries =
  739.                     planetaryEpsilonParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);

  741.             @SuppressWarnings("unchecked")
  742.             final PoissonSeries.CompiledSeries<DerivativeStructure> luniSolarSeries =
  743.                     PoissonSeries.compile(psiLuniSolarSeries, epsilonLuniSolarSeries);
  744.             @SuppressWarnings("unchecked")
  745.             final PoissonSeries.CompiledSeries<DerivativeStructure> planetarySeries =
  746.                     PoissonSeries.compile(psiPlanetarySeries, epsilonPlanetarySeries);

  748.             return new TimeFunction<double[]>() {
  749.                 /** {@inheritDoc} */
  750.                 @Override
  751.                 public double[] value(final AbsoluteDate date) {
  752.                     final BodiesElements elements = arguments.evaluateAll(date);
  753.                     final double[] luniSolar = luniSolarSeries.value(elements);
  754.                     final double[] planetary = planetarySeries.value(elements);
  755.                     return new double[] {
  756.                         luniSolar[0] + planetary[0], luniSolar[1] + planetary[1],
  757.                         IAU1994ResolutionC7.value(elements)
  758.                     };
  759.                 }
  760.             };

  762.         }

  764.         /** {@inheritDoc} */
  765.         @Override
  766.         public TimeFunction<DerivativeStructure> getGMSTFunction(final TimeScale ut1)
  767.             throws OrekitException {

  769.             // Earth Rotation Angle
  770.             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);

  772.             // Polynomial part of the apparent sidereal time series
  773.             // which is the opposite of Equation of Origins (EO)
  774.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  775.             final PoissonSeriesParser<DerivativeStructure> parser =
  776.                     new PoissonSeriesParser<DerivativeStructure>(17).
  777.                         withFirstDelaunay(4).
  778.                         withFirstPlanetary(9).
  779.                         withSinCos(0, 2, microAS, 3, microAS).
  780.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  781.             final PolynomialNutation<DerivativeStructure> minusEO =
  782.                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();

  784.             // create a function evaluating the series
  785.             return new TimeFunction<DerivativeStructure>() {

  787.                 /** {@inheritDoc} */
  788.                 @Override
  789.                 public DerivativeStructure value(final AbsoluteDate date) {
  790.                     return era.value(date).add(minusEO.value(dsEvaluateTC(date)));
  791.                 }

  793.             };

  795.         }

  797.         /** {@inheritDoc} */
  798.         @Override
  799.         public TimeFunction<DerivativeStructure> getGASTFunction(final TimeScale ut1,
  800.                                                                  final EOPHistory eopHistory)
  801.             throws OrekitException {

  803.             // set up nutation arguments
  804.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  806.             // mean obliquity function
  807.             final TimeFunction<Double> epsilon = getMeanObliquityFunction();

  809.             // set up Poisson series
  810.             final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
  811.             final PoissonSeriesParser<DerivativeStructure> luniSolarPsiParser =
  812.                     new PoissonSeriesParser<DerivativeStructure>(14).
  813.                         withFirstDelaunay(1).
  814.                         withSinCos(0, 7, milliAS, 11, milliAS).
  815.                         withSinCos(1, 8, milliAS, 12, milliAS);
  816.             final PoissonSeries<DerivativeStructure> psiLuniSolarSeries =
  817.                     luniSolarPsiParser.parse(getStream(LUNI_SOLAR_SERIES), LUNI_SOLAR_SERIES);

  819.             final PoissonSeriesParser<DerivativeStructure> planetaryPsiParser =
  820.                     new PoissonSeriesParser<DerivativeStructure>(21).
  821.                         withFirstDelaunay(2).
  822.                         withFirstPlanetary(7).
  823.                         withSinCos(0, 17, milliAS, 18, milliAS);
  824.             final PoissonSeries<DerivativeStructure> psiPlanetarySeries =
  825.                     planetaryPsiParser.parse(getStream(PLANETARY_SERIES), PLANETARY_SERIES);

  827.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  828.             final PoissonSeriesParser<DerivativeStructure> gstParser =
  829.                     new PoissonSeriesParser<DerivativeStructure>(17).
  830.                         withFirstDelaunay(4).
  831.                         withFirstPlanetary(9).
  832.                         withSinCos(0, 2, microAS, 3, microAS).
  833.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  834.             final PoissonSeries<DerivativeStructure> gstSeries = gstParser.parse(getStream(GST_SERIES), GST_SERIES);
  835.             @SuppressWarnings("unchecked")
  836.             final PoissonSeries.CompiledSeries<DerivativeStructure> psiGstSeries =
  837.                     PoissonSeries.compile(psiLuniSolarSeries, psiPlanetarySeries, gstSeries);

  839.             // ERA function
  840.             final TimeFunction<DerivativeStructure> era = getEarthOrientationAngleFunction(ut1);

  842.             return new TimeFunction<DerivativeStructure>() {

  844.                 /** {@inheritDoc} */
  845.                 @Override
  846.                 public DerivativeStructure value(final AbsoluteDate date) {

  848.                     // evaluate equation of origins
  849.                     final BodiesElements elements = arguments.evaluateAll(date);
  850.                     final double[] angles = psiGstSeries.value(elements);
  851.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  852.                     final double deltaPsi = angles[0] + angles[1] + ddPsi;
  853.                     final double epsilonA = epsilon.value(date);

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

  859.                 }

  861.             };

  863.         }

  865.         /** {@inheritDoc} */
  866.         @Override
  867.         public TimeFunction<double[]> getEOPTidalCorrection()
  868.             throws OrekitException {

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

  877.             // set up Poisson series
  878.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  879.             final PoissonSeriesParser<DerivativeStructure> xyParser = new PoissonSeriesParser<DerivativeStructure>(13).
  880.                     withOptionalColumn(1).
  881.                     withGamma(2).
  882.                     withFirstDelaunay(3);
  883.             final PoissonSeries<DerivativeStructure> xSeries =
  884.                     xyParser.
  885.                     withSinCos(0, 10, microAS, 11, microAS).
  886.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
  887.             final PoissonSeries<DerivativeStructure> ySeries =
  888.                     xyParser.
  889.                     withSinCos(0, 12, microAS, 13, microAS).
  890.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);

  892.             final double microS = 1.0e-6;
  893.             final PoissonSeriesParser<DerivativeStructure> ut1Parser = new PoissonSeriesParser<DerivativeStructure>(11).
  894.                     withOptionalColumn(1).
  895.                     withGamma(2).
  896.                     withFirstDelaunay(3).
  897.                     withSinCos(0, 10, microS, 11, microS);
  898.             final PoissonSeries<DerivativeStructure> ut1Series =
  899.                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);

  901.             @SuppressWarnings("unchecked")
  902.             final PoissonSeries.CompiledSeries<DerivativeStructure> correctionSeries =
  903.                 PoissonSeries.compile(xSeries, ySeries, ut1Series);

  905.             return new TimeFunction<double[]>() {
  906.                 /** {@inheritDoc} */
  907.                 @Override
  908.                 public double[] value(final AbsoluteDate date) {
  909.                     final FieldBodiesElements<DerivativeStructure> elements =
  910.                             arguments.evaluateDerivative(date);
  911.                     final DerivativeStructure[] correction = correctionSeries.value(elements);
  912.                     return new double[] {
  913.                         correction[0].getValue(),
  914.                         correction[1].getValue(),
  915.                         correction[2].getValue(),
  916.                         -correction[2].getPartialDerivative(1) * Constants.JULIAN_DAY
  917.                     };
  918.                 }
  919.             };

  921.         }

  923.         /** {@inheritDoc} */
  924.         public LoveNumbers getLoveNumbers() throws OrekitException {
  925.             return loadLoveNumbers(LOVE_NUMBERS);
  926.         }

  928.         /** {@inheritDoc} */
  929.         public TimeFunction<double[]> getTideFrequencyDependenceFunction(final TimeScale ut1)
  930.             throws OrekitException {

  932.             // set up nutation arguments
  933.             final FundamentalNutationArguments arguments = getNutationArguments(ut1);

  935.             // set up Poisson series
  936.             final PoissonSeriesParser<DerivativeStructure> k20Parser =
  937.                     new PoissonSeriesParser<DerivativeStructure>(18).
  938.                         withOptionalColumn(1).
  939.                         withDoodson(4, 2).
  940.                         withFirstDelaunay(10);
  941.             final PoissonSeriesParser<DerivativeStructure> k21Parser =
  942.                     new PoissonSeriesParser<DerivativeStructure>(18).
  943.                         withOptionalColumn(1).
  944.                         withDoodson(4, 3).
  945.                         withFirstDelaunay(10);
  946.             final PoissonSeriesParser<DerivativeStructure> k22Parser =
  947.                     new PoissonSeriesParser<DerivativeStructure>(16).
  948.                         withOptionalColumn(1).
  949.                         withDoodson(4, 2).
  950.                         withFirstDelaunay(10);

  952.             final double pico = 1.0e-12;
  953.             final PoissonSeries<DerivativeStructure> c20Series =
  954.                     k20Parser.
  955.                     // TODO: check sign!
  956. //                  withSinCos(0, 18, pico, 16, pico).
  957.                   withSinCos(0, 18, -pico, 16, pico).
  958.                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
  959.             final PoissonSeries<DerivativeStructure> c21Series =
  960.                     k21Parser.
  961.                     withSinCos(0, 17, pico, 18, pico).
  962.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  963.             final PoissonSeries<DerivativeStructure> s21Series =
  964.                     k21Parser.
  965.                     withSinCos(0, 18, -pico, 17, pico).
  966.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  967.             final PoissonSeries<DerivativeStructure> c22Series =
  968.                     k22Parser.
  969.                     withSinCos(0, -1, pico, 16, pico).
  970.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
  971.             final PoissonSeries<DerivativeStructure> s22Series =
  972.                     k22Parser.
  973.                     withSinCos(0, 16, -pico, -1, pico).
  974.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);

  976.             @SuppressWarnings("unchecked")
  977.             final PoissonSeries.CompiledSeries<DerivativeStructure> kSeries =
  978.                 PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);

  980.             return new TimeFunction<double[]>() {
  981.                 /** {@inheritDoc} */
  982.                 @Override
  983.                 public double[] value(final AbsoluteDate date) {
  984.                     return kSeries.value(arguments.evaluateAll(date));
  985.                 }
  986.             };

  988.         }

  990.         /** {@inheritDoc} */
  991.         @Override
  992.         public double getPermanentTide() throws OrekitException {
  993.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  994.         }

  996.         /** {@inheritDoc} */
  997.         @Override
  998.         public TimeFunction<double[]> getSolidPoleTide(final EOPHistory eopHistory)
  999.             throws OrekitException {

  1001.             // annual pole from ftp://tai.bipm.org/iers/conv2003/chapter7/annual.pole
  1002.             final TimeScale utc = TimeScalesFactory.getUTC();
  1003.             final SimpleTimeStampedTableParser.RowConverter<MeanPole> converter =
  1004.                 new SimpleTimeStampedTableParser.RowConverter<MeanPole>() {
  1005.                     /** {@inheritDoc} */
  1006.                     @Override
  1007.                     public MeanPole convert(final double[] rawFields) throws OrekitException {
  1008.                         return new MeanPole(new AbsoluteDate((int) rawFields[0], 1, 1, utc),
  1009.                                             rawFields[1] * Constants.ARC_SECONDS_TO_RADIANS,
  1010.                                             rawFields[2] * Constants.ARC_SECONDS_TO_RADIANS);
  1011.                     }
  1012.                 };
  1013.             final SimpleTimeStampedTableParser<MeanPole> parser =
  1014.                     new SimpleTimeStampedTableParser<MeanPole>(3, converter);
  1015.             final List<MeanPole> annualPoleList = parser.parse(getStream(ANNUAL_POLE), ANNUAL_POLE);
  1016.             final AbsoluteDate firstAnnualPoleDate = annualPoleList.get(0).getDate();
  1017.             final AbsoluteDate lastAnnualPoleDate  = annualPoleList.get(annualPoleList.size() - 1).getDate();
  1018.             final ImmutableTimeStampedCache<MeanPole> annualCache =
  1019.                     new ImmutableTimeStampedCache<MeanPole>(2, annualPoleList);

  1021.             // polynomial extension from IERS 2003, section 7.1.4, equations 23a and 23b
  1022.             final double xp0    = 0.054   * Constants.ARC_SECONDS_TO_RADIANS;
  1023.             final double xp0Dot = 0.00083 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
  1024.             final double yp0    = 0.357   * Constants.ARC_SECONDS_TO_RADIANS;
  1025.             final double yp0Dot = 0.00395 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;

  1027.             // constants from IERS 2003, section 6.2
  1028.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1029.             final double ratio        =  0.00115;

  1031.             return new TimeFunction<double[]>() {
  1032.                 /** {@inheritDoc} */
  1033.                 @Override
  1034.                 public double[] value(final AbsoluteDate date) {

  1036.                     // we can't compute anything before the range covered by the annual pole file
  1037.                     if (date.compareTo(firstAnnualPoleDate) <= 0) {
  1038.                         return new double[] {
  1039.                             0.0, 0.0
  1040.                         };
  1041.                     }

  1043.                     // evaluate mean pole
  1044.                     double meanPoleX = 0;
  1045.                     double meanPoleY = 0;
  1046.                     if (date.compareTo(lastAnnualPoleDate) <= 0) {
  1047.                         // we are within the range covered by the annual pole file,
  1048.                         // we interpolate within it
  1049.                         try {
  1050.                             final List<MeanPole> neighbors = annualCache.getNeighbors(date);
  1051.                             final HermiteInterpolator interpolator = new HermiteInterpolator();
  1052.                             for (final MeanPole neighbor : neighbors) {
  1053.                                 interpolator.addSamplePoint(neighbor.getDate().durationFrom(date),
  1054.                                                             new double[] {
  1055.                                                                 neighbor.getX(), neighbor.getY()
  1056.                                                             });
  1057.                             }
  1058.                             final double[] interpolated = interpolator.value(0);
  1059.                             meanPoleX = interpolated[0];
  1060.                             meanPoleY = interpolated[1];
  1061.                         } catch (TimeStampedCacheException tsce) {
  1062.                             // this should never happen
  1063.                             throw OrekitException.createInternalError(tsce);
  1064.                         }
  1065.                     } else {

  1067.                         // we are after the range covered by the annual pole file,
  1068.                         // we use the polynomial extension
  1069.                         final double t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
  1070.                         meanPoleX = xp0 + t * xp0Dot;
  1071.                         meanPoleY = yp0 + t * yp0Dot;

  1073.                     }

  1075.                     // evaluate wobble variables
  1076.                     final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1077.                     final double m1 = correction.getXp() - meanPoleX;
  1078.                     final double m2 = meanPoleY - correction.getYp();

  1080.                     return new double[] {
  1081.                         // the following correspond to the equations published in IERS 2003 conventions,
  1082.                         // section 6.2 page 65. In the publication, the equations read:
  1083.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1084.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1085.                         // However, it seems there are sign errors in these equations, which have
  1086.                         // been fixed in IERS 2010 conventions, section 6.4 page 94. In these newer
  1087.                         // publication, the equations read:
  1088.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1089.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1090.                         // the newer equations seem more consistent with the premises as the
  1091.                         // deformation due to the centrifugal potential has the form:
  1092.                         // −Ω²r²/2 sin 2θ Re [k₂(m₁ − im₂) exp(iλ)] where k₂ is the complex
  1093.                         // number 0.3077 + 0.0036i, so the real part in the previous equation is:
  1094.                         // A[Re(k₂) m₁ + Im(k₂) m₂)] cos λ + A[Re(k₂) m₂ - Im(k₂) m₁] sin λ
  1095.                         // identifying this with ∆C₂₁ cos λ + ∆S₂₁ sin λ we get:
  1096.                         // ∆C₂₁ = A Re(k₂) [m₁ + Im(k₂)/Re(k₂) m₂)]
  1097.                         // ∆S₂₁ = A Re(k₂) [m₂ - Im(k₂)/Re(k₂) m₁)]
  1098.                         // and Im(k₂)/Re(k₂) is very close to +0.00115
  1099.                         // As the equation as written in the IERS 2003 conventions are used in
  1100.                         // legacy systems, we have reproduced this alleged error here (and fixed it in
  1101.                         // the IERS 2010 conventions below) for validation purposes. We don't recommend
  1102.                         // using the IERS 2003 conventions for solid pole tide computation other than
  1103.                         // for validation or reproducibility of legacy applications behavior.
  1104.                         // As solid pole tide is small and as the sign change is on the smallest coefficient,
  1105.                         // the effect is quite small. A test case on a propagated orbit showed a position change
  1106.                         // slightly below 0.4m after a 30 days propagation on a Low Earth Orbit
  1107.                         globalFactor * (m1 - ratio * m2),
  1108.                         globalFactor * (m2 + ratio * m1),
  1109.                     };

  1111.                 }
  1112.             };

  1114.         }

  1116.         /** {@inheritDoc} */
  1117.         @Override
  1118.         public TimeFunction<double[]> getOceanPoleTide(final EOPHistory eopHistory)
  1119.             throws OrekitException {

  1121.             return new TimeFunction<double[]>() {
  1122.                 /** {@inheritDoc} */
  1123.                 @Override
  1124.                 public double[] value(final AbsoluteDate date) {
  1125.                     // there are no model for ocean pole tide prior to conventions 2010
  1126.                     return new double[] {
  1127.                         0.0, 0.0
  1128.                     };
  1129.                 }
  1130.             };
  1131.         }

  1133.     },

  1135.     /** Constant for IERS 2010 conventions. */
  1136.     IERS_2010 {

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

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

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

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

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

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

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

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

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

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

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

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

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

  1177.         /** {@inheritDoc} */
  1178.         public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
  1179.             throws OrekitException {
  1180.             return new FundamentalNutationArguments(this, timeScale,
  1181.                                                     getStream(NUTATION_ARGUMENTS), NUTATION_ARGUMENTS);
  1182.         }

  1184.         /** {@inheritDoc} */
  1185.         @Override
  1186.         public TimeFunction<Double> getMeanObliquityFunction() throws OrekitException {

  1188.             // epsilon 0 value from chapter 5, page 56, other terms from equation 5.40 page 65
  1189.             final PolynomialNutation<DerivativeStructure> epsilonA =
  1190.                     new PolynomialNutation<DerivativeStructure>(84381.406        * Constants.ARC_SECONDS_TO_RADIANS,
  1191.                                                                   -46.836769     * Constants.ARC_SECONDS_TO_RADIANS,
  1192.                                                                    -0.0001831    * Constants.ARC_SECONDS_TO_RADIANS,
  1193.                                                                     0.00200340   * Constants.ARC_SECONDS_TO_RADIANS,
  1194.                                                                    -0.000000576  * Constants.ARC_SECONDS_TO_RADIANS,
  1195.                                                                    -0.0000000434 * Constants.ARC_SECONDS_TO_RADIANS);

  1197.             return new TimeFunction<Double>() {

  1199.                 /** {@inheritDoc} */
  1200.                 @Override
  1201.                 public Double value(final AbsoluteDate date) {
  1202.                     return epsilonA.value(evaluateTC(date));
  1203.                 }

  1205.             };

  1207.         }

  1209.         /** {@inheritDoc} */
  1210.         @Override
  1211.         public TimeFunction<double[]> getXYSpXY2Function() throws OrekitException {

  1213.             // set up nutation arguments
  1214.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  1216.             // set up Poisson series
  1217.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1218.             final PoissonSeriesParser<DerivativeStructure> parser =
  1219.                     new PoissonSeriesParser<DerivativeStructure>(17).
  1220.                         withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
  1221.                         withFirstDelaunay(4).
  1222.                         withFirstPlanetary(9).
  1223.                         withSinCos(0, 2, microAS, 3, microAS);
  1224.             final PoissonSeries<DerivativeStructure> xSeries = parser.parse(getStream(X_SERIES), X_SERIES);
  1225.             final PoissonSeries<DerivativeStructure> ySeries = parser.parse(getStream(Y_SERIES), Y_SERIES);
  1226.             final PoissonSeries<DerivativeStructure> sSeries = parser.parse(getStream(S_SERIES), S_SERIES);
  1227.             @SuppressWarnings("unchecked")
  1228.             final PoissonSeries.CompiledSeries<DerivativeStructure> xys = PoissonSeries.compile(xSeries, ySeries, sSeries);

  1230.             // create a function evaluating the series
  1231.             return new TimeFunction<double[]>() {

  1233.                 /** {@inheritDoc} */
  1234.                 @Override
  1235.                 public double[] value(final AbsoluteDate date) {
  1236.                     return xys.value(arguments.evaluateAll(date));
  1237.                 }

  1239.             };

  1241.         }

  1243.         /** {@inheritDoc} */
  1244.         public LoveNumbers getLoveNumbers() throws OrekitException {
  1245.             return loadLoveNumbers(LOVE_NUMBERS);
  1246.         }

  1248.         /** {@inheritDoc} */
  1249.         public TimeFunction<double[]> getTideFrequencyDependenceFunction(final TimeScale ut1)
  1250.             throws OrekitException {

  1252.             // set up nutation arguments
  1253.             final FundamentalNutationArguments arguments = getNutationArguments(ut1);

  1255.             // set up Poisson series
  1256.             final PoissonSeriesParser<DerivativeStructure> k20Parser =
  1257.                     new PoissonSeriesParser<DerivativeStructure>(18).
  1258.                         withOptionalColumn(1).
  1259.                         withDoodson(4, 2).
  1260.                         withFirstDelaunay(10);
  1261.             final PoissonSeriesParser<DerivativeStructure> k21Parser =
  1262.                     new PoissonSeriesParser<DerivativeStructure>(18).
  1263.                         withOptionalColumn(1).
  1264.                         withDoodson(4, 3).
  1265.                         withFirstDelaunay(10);
  1266.             final PoissonSeriesParser<DerivativeStructure> k22Parser =
  1267.                     new PoissonSeriesParser<DerivativeStructure>(16).
  1268.                         withOptionalColumn(1).
  1269.                         withDoodson(4, 2).
  1270.                         withFirstDelaunay(10);

  1272.             final double pico = 1.0e-12;
  1273.             final PoissonSeries<DerivativeStructure> c20Series =
  1274.                     k20Parser.
  1275.                     // TODO: check sign!
  1276. //                    withSinCos(0, 18, pico, 16, pico).
  1277.                     withSinCos(0, 18, -pico, 16, pico).
  1278.                     parse(getStream(K20_FREQUENCY_DEPENDENCE), K20_FREQUENCY_DEPENDENCE);
  1279.             final PoissonSeries<DerivativeStructure> c21Series =
  1280.                     k21Parser.
  1281.                     withSinCos(0, 17, pico, 18, pico).
  1282.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  1283.             final PoissonSeries<DerivativeStructure> s21Series =
  1284.                     k21Parser.
  1285.                     withSinCos(0, 18, -pico, 17, pico).
  1286.                     parse(getStream(K21_FREQUENCY_DEPENDENCE), K21_FREQUENCY_DEPENDENCE);
  1287.             final PoissonSeries<DerivativeStructure> c22Series =
  1288.                     k22Parser.
  1289.                     withSinCos(0, -1, pico, 16, pico).
  1290.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);
  1291.             final PoissonSeries<DerivativeStructure> s22Series =
  1292.                     k22Parser.
  1293.                     withSinCos(0, 16, -pico, -1, pico).
  1294.                     parse(getStream(K22_FREQUENCY_DEPENDENCE), K22_FREQUENCY_DEPENDENCE);

  1296.             @SuppressWarnings("unchecked")
  1297.             final PoissonSeries.CompiledSeries<DerivativeStructure> kSeries =
  1298.                 PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);

  1300.             return new TimeFunction<double[]>() {
  1301.                 /** {@inheritDoc} */
  1302.                 @Override
  1303.                 public double[] value(final AbsoluteDate date) {
  1304.                     return kSeries.value(arguments.evaluateAll(date));
  1305.                 }
  1306.             };

  1308.         }

  1310.         /** {@inheritDoc} */
  1311.         @Override
  1312.         public double getPermanentTide() throws OrekitException {
  1313.             return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
  1314.         }

  1316.         /** Compute pole wobble variables m₁ and m₂.
  1317.          * @param date current date
  1318.          * @param eopHistory EOP history
  1319.          * @return array containing m₁ and m₂
  1320.          */
  1321.         private double[] computePoleWobble(final AbsoluteDate date, final EOPHistory eopHistory) {

  1323.             // polynomial model from IERS 2010, table 7.7
  1324.             final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
  1325.             final double f1 = f0 / Constants.JULIAN_YEAR;
  1326.             final double f2 = f1 / Constants.JULIAN_YEAR;
  1327.             final double f3 = f2 / Constants.JULIAN_YEAR;
  1328.             final AbsoluteDate changeDate = new AbsoluteDate(2010, 1, 1, TimeScalesFactory.getTT());

  1330.             // evaluate mean pole
  1331.             final double[] xPolynomial;
  1332.             final double[] yPolynomial;
  1333.             if (date.compareTo(changeDate) <= 0) {
  1334.                 xPolynomial = new double[] {
  1335.                     55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
  1336.                 };
  1337.                 yPolynomial = new double[] {
  1338.                     346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
  1339.                 };
  1340.             } else {
  1341.                 xPolynomial = new double[] {
  1342.                     23.513 * f0, 7.6141 * f1
  1343.                 };
  1344.                 yPolynomial = new double[] {
  1345.                     358.891 * f0,  -0.6287 * f1
  1346.                 };
  1347.             }
  1348.             double meanPoleX = 0;
  1349.             double meanPoleY = 0;
  1350.             final double t = date.durationFrom(AbsoluteDate.J2000_EPOCH);
  1351.             for (int i = xPolynomial.length - 1; i >= 0; --i) {
  1352.                 meanPoleX = meanPoleX * t + xPolynomial[i];
  1353.             }
  1354.             for (int i = yPolynomial.length - 1; i >= 0; --i) {
  1355.                 meanPoleY = meanPoleY * t + yPolynomial[i];
  1356.             }

  1358.             // evaluate wobble variables
  1359.             final PoleCorrection correction = eopHistory.getPoleCorrection(date);
  1360.             final double m1 = correction.getXp() - meanPoleX;
  1361.             final double m2 = meanPoleY - correction.getYp();

  1363.             return new double[] {
  1364.                 m1, m2
  1365.             };

  1367.         }

  1369.         /** {@inheritDoc} */
  1370.         @Override
  1371.         public TimeFunction<double[]> getSolidPoleTide(final EOPHistory eopHistory)
  1372.             throws OrekitException {

  1374.             // constants from IERS 2010, section 6.4
  1375.             final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
  1376.             final double ratio        =  0.00115;

  1378.             return new TimeFunction<double[]>() {
  1379.                 /** {@inheritDoc} */
  1380.                 @Override
  1381.                 public double[] value(final AbsoluteDate date) {

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

  1386.                     return new double[] {
  1387.                         // the following correspond to the equations published in IERS 2010 conventions,
  1388.                         // section 6.4 page 94. The equations read:
  1389.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
  1390.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
  1391.                         // These equations seem to fix what was probably a sign error in IERS 2003
  1392.                         // conventions section 6.2 page 65. In this older publication, the equations read:
  1393.                         // ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
  1394.                         // ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
  1395.                         globalFactor * (wobbleM[0] + ratio * wobbleM[1]),
  1396.                         globalFactor * (wobbleM[1] - ratio * wobbleM[0])
  1397.                     };

  1399.                 }
  1400.             };

  1402.         }

  1404.         /** {@inheritDoc} */
  1405.         @Override
  1406.         public TimeFunction<double[]> getOceanPoleTide(final EOPHistory eopHistory)
  1407.             throws OrekitException {

  1409.             return new TimeFunction<double[]>() {
  1410.                 /** {@inheritDoc} */
  1411.                 @Override
  1412.                 public double[] value(final AbsoluteDate date) {

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

  1417.                     return new double[] {
  1418.                         // the following correspond to the equations published in IERS 2010 conventions,
  1419.                         // section 6.4 page 94 equation 6.24:
  1420.                         // ∆C₂₁ = −2.1778 × 10⁻¹⁰ (m₁ − 0.01724m₂)
  1421.                         // ∆S₂₁ = −1.7232 × 10⁻¹⁰ (m₂ − 0.03365m₁)
  1422.                         -2.1778e-10 * (wobbleM[0] - 0.01724 * wobbleM[1]) / Constants.ARC_SECONDS_TO_RADIANS,
  1423.                         -1.7232e-10 * (wobbleM[1] - 0.03365 * wobbleM[0]) / Constants.ARC_SECONDS_TO_RADIANS
  1424.                     };

  1426.                 }
  1427.             };

  1429.         }

  1431.         /** {@inheritDoc} */
  1432.         @Override
  1433.         public TimeFunction<double[]> getPrecessionFunction() throws OrekitException {

  1435.             // set up the conventional polynomials
  1436.             // the following values are from equation 5.40 in IERS 2010 conventions
  1437.             final PolynomialNutation<DerivativeStructure> psiA =
  1438.                     new PolynomialNutation<DerivativeStructure>(   0.0,
  1439.                                                                 5038.481507     * Constants.ARC_SECONDS_TO_RADIANS,
  1440.                                                                   -1.0790069    * Constants.ARC_SECONDS_TO_RADIANS,
  1441.                                                                   -0.00114045   * Constants.ARC_SECONDS_TO_RADIANS,
  1442.                                                                    0.000132851  * Constants.ARC_SECONDS_TO_RADIANS,
  1443.                                                                   -0.0000000951 * Constants.ARC_SECONDS_TO_RADIANS);
  1444.             final PolynomialNutation<DerivativeStructure> omegaA =
  1445.                     new PolynomialNutation<DerivativeStructure>(getMeanObliquityFunction().value(getNutationReferenceEpoch()),
  1446.                                                                 -0.025754     * Constants.ARC_SECONDS_TO_RADIANS,
  1447.                                                                  0.0512623    * Constants.ARC_SECONDS_TO_RADIANS,
  1448.                                                                 -0.00772503   * Constants.ARC_SECONDS_TO_RADIANS,
  1449.                                                                 -0.000000467  * Constants.ARC_SECONDS_TO_RADIANS,
  1450.                                                                  0.0000003337 * Constants.ARC_SECONDS_TO_RADIANS);
  1451.             final PolynomialNutation<DerivativeStructure> chiA =
  1452.                     new PolynomialNutation<DerivativeStructure>( 0.0,
  1453.                                                                 10.556403     * Constants.ARC_SECONDS_TO_RADIANS,
  1454.                                                                 -2.3814292    * Constants.ARC_SECONDS_TO_RADIANS,
  1455.                                                                 -0.00121197   * Constants.ARC_SECONDS_TO_RADIANS,
  1456.                                                                  0.000170663  * Constants.ARC_SECONDS_TO_RADIANS,
  1457.                                                                 -0.0000000560 * Constants.ARC_SECONDS_TO_RADIANS);

  1459.             return new TimeFunction<double[]>() {
  1460.                 /** {@inheritDoc} */
  1461.                 @Override
  1462.                 public double[] value(final AbsoluteDate date) {
  1463.                     final double tc = evaluateTC(date);
  1464.                     return new double[] {
  1465.                         psiA.value(tc), omegaA.value(tc), chiA.value(tc)
  1466.                     };
  1467.                 }
  1468.             };

  1470.         }

  1472.          /** {@inheritDoc} */
  1473.         @Override
  1474.         public TimeFunction<double[]> getNutationFunction()
  1475.             throws OrekitException {

  1477.             // set up nutation arguments
  1478.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  1480.             // set up Poisson series
  1481.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1482.             final PoissonSeriesParser<DerivativeStructure> parser =
  1483.                     new PoissonSeriesParser<DerivativeStructure>(17).
  1484.                         withFirstDelaunay(4).
  1485.                         withFirstPlanetary(9).
  1486.                         withSinCos(0, 2, microAS, 3, microAS);
  1487.             final PoissonSeries<DerivativeStructure> psiSeries     = parser.parse(getStream(PSI_SERIES), PSI_SERIES);
  1488.             final PoissonSeries<DerivativeStructure> epsilonSeries = parser.parse(getStream(EPSILON_SERIES), EPSILON_SERIES);
  1489.             @SuppressWarnings("unchecked")
  1490.             final PoissonSeries.CompiledSeries<DerivativeStructure> psiEpsilonSeries =
  1491.                     PoissonSeries.compile(psiSeries, epsilonSeries);

  1493.             return new TimeFunction<double[]>() {
  1494.                 /** {@inheritDoc} */
  1495.                 @Override
  1496.                 public double[] value(final AbsoluteDate date) {
  1497.                     final BodiesElements elements = arguments.evaluateAll(date);
  1498.                     final double[] psiEpsilon = psiEpsilonSeries.value(elements);
  1499.                     return new double[] {
  1500.                         psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements)
  1501.                     };
  1502.                 }
  1503.             };

  1505.         }

  1507.         /** {@inheritDoc} */
  1508.         @Override
  1509.         public TimeFunction<DerivativeStructure> getGMSTFunction(final TimeScale ut1) throws OrekitException {

  1511.             // Earth Rotation Angle
  1512.             final StellarAngleCapitaine era = new StellarAngleCapitaine(ut1);

  1514.             // Polynomial part of the apparent sidereal time series
  1515.             // which is the opposite of Equation of Origins (EO)
  1516.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1517.             final PoissonSeriesParser<DerivativeStructure> parser =
  1518.                     new PoissonSeriesParser<DerivativeStructure>(17).
  1519.                         withFirstDelaunay(4).
  1520.                         withFirstPlanetary(9).
  1521.                         withSinCos(0, 2, microAS, 3, microAS).
  1522.                         withPolynomialPart('t', Unit.ARC_SECONDS);
  1523.             final PolynomialNutation<DerivativeStructure> minusEO =
  1524.                     parser.parse(getStream(GST_SERIES), GST_SERIES).getPolynomial();

  1526.             // create a function evaluating the series
  1527.             return new TimeFunction<DerivativeStructure>() {

  1529.                 /** {@inheritDoc} */
  1530.                 @Override
  1531.                 public DerivativeStructure value(final AbsoluteDate date) {
  1532.                     return era.value(date).add(minusEO.value(dsEvaluateTC(date)));
  1533.                 }

  1535.             };

  1537.         }

  1539.         /** {@inheritDoc} */
  1540.         @Override
  1541.         public TimeFunction<DerivativeStructure> getGASTFunction(final TimeScale ut1,
  1542.                                                                  final EOPHistory eopHistory)
  1543.             throws OrekitException {

  1545.             // set up nutation arguments
  1546.             final FundamentalNutationArguments arguments = getNutationArguments(null);

  1548.             // mean obliquity function
  1549.             final TimeFunction<Double> epsilon = getMeanObliquityFunction();

  1551.             // set up Poisson series
  1552.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1553.             final PoissonSeriesParser<DerivativeStructure> baseParser =
  1554.                     new PoissonSeriesParser<DerivativeStructure>(17).
  1555.                         withFirstDelaunay(4).
  1556.                         withFirstPlanetary(9).
  1557.                         withSinCos(0, 2, microAS, 3, microAS);
  1558.             final PoissonSeriesParser<DerivativeStructure> gstParser  = baseParser.withPolynomialPart('t', Unit.ARC_SECONDS);
  1559.             final PoissonSeries<DerivativeStructure> psiSeries        = baseParser.parse(getStream(PSI_SERIES), PSI_SERIES);
  1560.             final PoissonSeries<DerivativeStructure> gstSeries        = gstParser.parse(getStream(GST_SERIES), GST_SERIES);
  1561.             @SuppressWarnings("unchecked")
  1562.             final PoissonSeries.CompiledSeries<DerivativeStructure> psiGstSeries =
  1563.                     PoissonSeries.compile(psiSeries, gstSeries);

  1565.             // ERA function
  1566.             final TimeFunction<DerivativeStructure> era = getEarthOrientationAngleFunction(ut1);

  1568.             return new TimeFunction<DerivativeStructure>() {

  1570.                 /** {@inheritDoc} */
  1571.                 @Override
  1572.                 public DerivativeStructure value(final AbsoluteDate date) {

  1574.                     // evaluate equation of origins
  1575.                     final BodiesElements elements = arguments.evaluateAll(date);
  1576.                     final double[] angles = psiGstSeries.value(elements);
  1577.                     final double ddPsi    = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
  1578.                     final double deltaPsi = angles[0] + ddPsi;
  1579.                     final double epsilonA = epsilon.value(date);

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

  1585.                 }

  1587.             };

  1589.         }

  1591.         /** {@inheritDoc} */
  1592.         @Override
  1593.         public TimeFunction<double[]> getEOPTidalCorrection()
  1594.             throws OrekitException {

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

  1603.             // set up Poisson series
  1604.             final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
  1605.             final PoissonSeriesParser<DerivativeStructure> xyParser = new PoissonSeriesParser<DerivativeStructure>(13).
  1606.                     withOptionalColumn(1).
  1607.                     withGamma(2).
  1608.                     withFirstDelaunay(3);
  1609.             final PoissonSeries<DerivativeStructure> xSeries =
  1610.                     xyParser.
  1611.                     withSinCos(0, 10, microAS, 11, microAS).
  1612.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);
  1613.             final PoissonSeries<DerivativeStructure> ySeries =
  1614.                     xyParser.
  1615.                     withSinCos(0, 12, microAS, 13, microAS).
  1616.                     parse(getStream(TIDAL_CORRECTION_XP_YP_SERIES), TIDAL_CORRECTION_XP_YP_SERIES);

  1618.             final double microS = 1.0e-6;
  1619.             final PoissonSeriesParser<DerivativeStructure> ut1Parser = new PoissonSeriesParser<DerivativeStructure>(11).
  1620.                     withOptionalColumn(1).
  1621.                     withGamma(2).
  1622.                     withFirstDelaunay(3).
  1623.                     withSinCos(0, 10, microS, 11, microS);
  1624.             final PoissonSeries<DerivativeStructure> ut1Series =
  1625.                     ut1Parser.parse(getStream(TIDAL_CORRECTION_UT1_SERIES), TIDAL_CORRECTION_UT1_SERIES);

  1627.             @SuppressWarnings("unchecked")
  1628.             final PoissonSeries.CompiledSeries<DerivativeStructure> correctionSeries =
  1629.                     PoissonSeries.compile(xSeries, ySeries, ut1Series);

  1631.             return new TimeFunction<double[]>() {
  1632.                 /** {@inheritDoc} */
  1633.                 @Override
  1634.                 public double[] value(final AbsoluteDate date) {
  1635.                     final FieldBodiesElements<DerivativeStructure> elements =
  1636.                             arguments.evaluateDerivative(date);
  1637.                     final DerivativeStructure[] correction = correctionSeries.value(elements);
  1638.                     return new double[] {
  1639.                         correction[0].getValue(),
  1640.                         correction[1].getValue(),
  1641.                         correction[2].getValue(),
  1642.                         -correction[2].getPartialDerivative(1) * Constants.JULIAN_DAY
  1643.                     };
  1644.                 }
  1645.             };

  1647.         }

  1649.     };

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

  1654.     /** Get the reference epoch for fundamental nutation arguments.
  1655.      * @return reference epoch for fundamental nutation arguments
  1656.      * @since 6.1
  1657.      */
  1658.     public AbsoluteDate getNutationReferenceEpoch() {
  1659.         // IERS 1996, IERS 2003 and IERS 2010 use the same J2000.0 reference date
  1660.         return AbsoluteDate.J2000_EPOCH;
  1661.     }

  1663.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  1664.      * @param date current date
  1665.      * @return date offset in Julian centuries
  1666.      * @since 6.1
  1667.      */
  1668.     public double evaluateTC(final AbsoluteDate date) {
  1669.         return date.durationFrom(getNutationReferenceEpoch()) / Constants.JULIAN_CENTURY;
  1670.     }

  1672.     /** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
  1673.      * @param date current date
  1674.      * @return date offset in Julian centuries
  1675.      * @since 6.1
  1676.      */
  1677.     public DerivativeStructure dsEvaluateTC(final AbsoluteDate date) {
  1678.         return new DerivativeStructure(1, 1, evaluateTC(date), 1.0 / Constants.JULIAN_CENTURY);
  1679.     }

  1681.     /** Get the fundamental nutation arguments.
  1682.      * @param timeScale time scale for computing Greenwich Mean Sidereal Time
  1683.      * (typically {@link TimeScalesFactory#getUT1(IERSConventions, boolean) UT1})
  1684.      * @return fundamental nutation arguments
  1685.      * @exception OrekitException if fundamental nutation arguments cannot be loaded
  1686.      * @since 6.1
  1687.      */
  1688.     public abstract FundamentalNutationArguments getNutationArguments(final TimeScale timeScale)
  1689.         throws OrekitException;

  1691.     /** Get the function computing mean obliquity of the ecliptic.
  1692.      * @return function computing mean obliquity of the ecliptic
  1693.      * @exception OrekitException if table cannot be loaded
  1694.      * @since 6.1
  1695.      */
  1696.     public abstract TimeFunction<Double> getMeanObliquityFunction() throws OrekitException;

  1698.     /** Get the function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components.
  1699.      * <p>
  1700.      * The returned function computes the two X, Y components of CIP and the S+XY/2 component of the non-rotating CIO.
  1701.      * </p>
  1702.      * @return function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components
  1703.      * @exception OrekitException if table cannot be loaded
  1704.      * @since 6.1
  1705.      */
  1706.     public abstract TimeFunction<double[]> getXYSpXY2Function()
  1707.         throws OrekitException;

  1709.     /** Get the function computing the raw Earth Orientation Angle.
  1710.      * <p>
  1711.      * The raw angle does not contain any correction. If for example dTU1 correction
  1712.      * due to tidal effect is desired, it must be added afterward by the caller.
  1713.      * The returned value contain the angle as the value and the angular rate as
  1714.      * the first derivative.
  1715.      * </p>
  1716.      * @param ut1 UT1 time scale
  1717.      * @return function computing the rawEarth Orientation Angle, in the non-rotating origin paradigm,
  1718.      * the return value containing both the angle and its first time derivative
  1719.      * @since 6.1
  1720.      */
  1721.     public TimeFunction<DerivativeStructure> getEarthOrientationAngleFunction(final TimeScale ut1) {
  1722.         return new StellarAngleCapitaine(ut1);
  1723.     }


  1726.     /** Get the function computing the precession angles.
  1727.      * <p>
  1728.      * The function returned computes the three precession angles
  1729.      * &psi;<sub>A</sub> (around Z axis), &omega;<sub>A</sub> (around X axis)
  1730.      * and &chi;<sub>A</sub> (around Z axis). The constant angle &epsilon;<sub>0</sub>
  1731.      * for the fourth rotation (around X axis) can be retrieved by evaluating the
  1732.      * function returned by {@link #getMeanObliquityFunction()} at {@link
  1733.      * #getNutationReferenceEpoch() nutation reference epoch}.
  1734.      * </p>
  1735.      * @return function computing the precession angle
  1736.      * @exception OrekitException if table cannot be loaded
  1737.      * @since 6.1
  1738.      */
  1739.     public abstract TimeFunction<double[]> getPrecessionFunction() throws OrekitException;

  1741.     /** Get the function computing the nutation angles.
  1742.      * <p>
  1743.      * The function returned computes the two classical angles &Delta;&Psi; and &Delta;&epsilon;,
  1744.      * and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
  1745.      * resolution C7 (the correction is forced to 0 before this date)
  1746.      * </p>
  1747.      * @return function computing the nutation in longitude &Delta;&Psi; and &Delta;&epsilon;
  1748.      * and the correction of equation of equinoxes
  1749.      * @exception OrekitException if table cannot be loaded
  1750.      * @since 6.1
  1751.      */
  1752.     public abstract TimeFunction<double[]> getNutationFunction()
  1753.         throws OrekitException;

  1755.     /** Get the function computing Greenwich mean sidereal time, in radians.
  1756.      * @param ut1 UT1 time scale
  1757.      * @return function computing Greenwich mean sidereal time,
  1758.      * the return value containing both the angle and its first time derivative
  1759.      * @exception OrekitException if table cannot be loaded
  1760.      * @since 6.1
  1761.      */
  1762.     public abstract TimeFunction<DerivativeStructure> getGMSTFunction(TimeScale ut1)
  1763.         throws OrekitException;

  1765.     /** Get the function computing Greenwich apparent sidereal time, in radians.
  1766.      * @param ut1 UT1 time scale
  1767.      * @param eopHistory EOP history
  1768.      * @return function computing Greenwich apparent sidereal time,
  1769.      * the return value containing both the angle and its first time derivative
  1770.      * @exception OrekitException if table cannot be loaded
  1771.      * @since 6.1
  1772.      */
  1773.     public abstract TimeFunction<DerivativeStructure> getGASTFunction(TimeScale ut1,
  1774.                                                                       EOPHistory eopHistory)
  1775.         throws OrekitException;

  1777.     /** Get the function computing tidal corrections for Earth Orientation Parameters.
  1778.      * @return function computing tidal corrections for Earth Orientation Parameters,
  1779.      * for xp, yp, ut1 and lod respectively
  1780.      * @exception OrekitException if table cannot be loaded
  1781.      * @since 6.1
  1782.      */
  1783.     public abstract TimeFunction<double[]> getEOPTidalCorrection()
  1784.         throws OrekitException;

  1786.     /** Get the Love numbers.
  1787.      * @return Love numbers
  1788.      * @exception OrekitException if table cannot be loaded
  1789.      * @since 6.1
  1790.      */
  1791.     public abstract LoveNumbers getLoveNumbers()
  1792.         throws OrekitException;

  1794.     /** Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  1795.      * @param ut1 UT1 time scale
  1796.      * @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  1797.      * @exception OrekitException if table cannot be loaded
  1798.      * @since 6.1
  1799.      */
  1800.     public abstract TimeFunction<double[]> getTideFrequencyDependenceFunction(TimeScale ut1)
  1801.         throws OrekitException;

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

  1809.     /** Get the function computing solid pole tide (ΔC₂₁, ΔS₂₁).
  1810.      * @param eopHistory EOP history
  1811.      * @return model for solid pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  1812.      * @exception OrekitException if table cannot be loaded
  1813.      * @since 6.1
  1814.      */
  1815.     public abstract TimeFunction<double[]> getSolidPoleTide(EOPHistory eopHistory)
  1816.         throws OrekitException;

  1818.     /** Get the function computing ocean pole tide (ΔC₂₁, ΔS₂₁).
  1819.      * @param eopHistory EOP history
  1820.      * @return model for ocean pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
  1821.      * @exception OrekitException if table cannot be loaded
  1822.      * @since 6.1
  1823.      */
  1824.     public abstract TimeFunction<double[]> getOceanPoleTide(EOPHistory eopHistory)
  1825.         throws OrekitException;

  1827.     /** Interface for functions converting nutation corrections between
  1828.      * &delta;&Delta;&psi;/&delta;&Delta;&epsilon; to &delta;X/&delta;Y.
  1829.      * <ul>
  1830.      * <li>&delta;&Delta;&psi;/&delta;&Delta;&epsilon; nutation corrections are used with the equinox-based paradigm.</li>
  1831.      * <li>&delta;X/&delta;Y nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  1832.      * </ul>
  1833.      * @since 6.1
  1834.      */
  1835.     public interface NutationCorrectionConverter {

  1837.         /** Convert nutation corrections.
  1838.          * @param date current date
  1839.          * @param ddPsi &delta;&Delta;&psi; part of the nutation correction
  1840.          * @param ddEpsilon &delta;&Delta;&epsilon; part of the nutation correction
  1841.          * @return array containing &delta;X and &delta;Y
  1842.          * @exception OrekitException if correction cannot be converted
  1843.          */
  1844.         double[] toNonRotating(AbsoluteDate date, double ddPsi, double ddEpsilon)
  1845.             throws OrekitException;

  1847.         /** Convert nutation corrections.
  1848.          * @param date current date
  1849.          * @param dX &delta;X part of the nutation correction
  1850.          * @param dY &delta;Y part of the nutation correction
  1851.          * @return array containing &delta;&Delta;&psi; and &delta;&Delta;&epsilon;
  1852.          * @exception OrekitException if correction cannot be converted
  1853.          */
  1854.         double[] toEquinox(AbsoluteDate date, double dX, double dY)
  1855.             throws OrekitException;

  1857.     }

  1859.     /** Create a function converting nutation corrections between
  1860.      * &delta;X/&delta;Y and &delta;&Delta;&psi;/&delta;&Delta;&epsilon;.
  1861.      * <ul>
  1862.      * <li>&delta;X/&delta;Y nutation corrections are used with the Non-Rotating Origin paradigm.</li>
  1863.      * <li>&delta;&Delta;&psi;/&delta;&Delta;&epsilon; nutation corrections are used with the equinox-based paradigm.</li>
  1864.      * </ul>
  1865.      * @return a new converter
  1866.      * @exception OrekitException if some convention table cannot be loaded
  1867.      * @since 6.1
  1868.      */
  1869.     public NutationCorrectionConverter getNutationCorrectionConverter()
  1870.         throws OrekitException {

  1872.         // get models parameters
  1873.         final TimeFunction<double[]> precessionFunction = getPrecessionFunction();
  1874.         final TimeFunction<Double> epsilonAFunction = getMeanObliquityFunction();
  1875.         final AbsoluteDate date0 = getNutationReferenceEpoch();
  1876.         final double cosE0 = FastMath.cos(epsilonAFunction.value(date0));

  1878.         return new NutationCorrectionConverter() {

  1880.             /** {@inheritDoc} */
  1881.             @Override
  1882.             public double[] toNonRotating(final AbsoluteDate date,
  1883.                                           final double ddPsi, final double ddEpsilon)
  1884.                 throws OrekitException {
  1885.                 // compute precession angles psiA, omegaA and chiA
  1886.                 final double[] angles = precessionFunction.value(date);

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

  1892.                 // convert nutation corrections (equation 23/IERS-2003 or 5.25/IERS-2010)
  1893.                 return new double[] {
  1894.                     sinEA * ddPsi + c * ddEpsilon,
  1895.                     ddEpsilon - c * sinEA * ddPsi
  1896.                 };

  1898.             }

  1900.             /** {@inheritDoc} */
  1901.             @Override
  1902.             public double[] toEquinox(final AbsoluteDate date,
  1903.                                       final double dX, final double dY)
  1904.                 throws OrekitException {
  1905.                 // compute precession angles psiA, omegaA and chiA
  1906.                 final double[] angles   = precessionFunction.value(date);

  1908.                 // conversion coefficients
  1909.                 final double sinEA = FastMath.sin(epsilonAFunction.value(date));
  1910.                 final double c     = angles[0] * cosE0 - angles[2];
  1911.                 final double opC2  = 1 + c * c;

  1913.                 // convert nutation corrections (inverse of equation 23/IERS-2003 or 5.25/IERS-2010)
  1914.                 return new double[] {
  1915.                     (dX - c * dY) / (sinEA * opC2),
  1916.                     (dY + c * dX) / opC2
  1917.                 };
  1918.             }

  1920.         };

  1922.     }

  1924.     /** Load the Love numbers.
  1925.      * @param nameLove name of the Love number resource
  1926.      * @return Love numbers
  1927.      * @exception OrekitException if the Love numbers embedded in the
  1928.      * library cannot be read
  1929.      */
  1930.     protected LoveNumbers loadLoveNumbers(final String nameLove) throws OrekitException {
  1931.         InputStream stream = null;
  1932.         BufferedReader reader = null;
  1933.         try {

  1935.             // allocate the three triangular arrays for real, imaginary and time-dependent numbers
  1936.             final double[][] real      = new double[4][];
  1937.             final double[][] imaginary = new double[4][];
  1938.             final double[][] plus      = new double[4][];
  1939.             for (int i = 0; i < real.length; ++i) {
  1940.                 real[i]      = new double[i + 1];
  1941.                 imaginary[i] = new double[i + 1];
  1942.                 plus[i]      = new double[i + 1];
  1943.             }

  1945.             stream = IERSConventions.class.getResourceAsStream(nameLove);
  1946.             if (stream == null) {
  1947.                 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, nameLove);
  1948.             }

  1950.             // setup the reader
  1951.             reader = new BufferedReader(new InputStreamReader(stream, "UTF-8"));

  1953.             String line = reader.readLine();
  1954.             int lineNumber = 1;

  1956.             // look for the Love numbers
  1957.             while (line != null) {

  1959.                 line = line.trim();
  1960.                 if (!(line.isEmpty() || line.startsWith("#"))) {
  1961.                     final String[] fields = line.split("\\p{Space}+");
  1962.                     if (fields.length != 5) {
  1963.                         // this should never happen with files embedded within Orekit
  1964.                         throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  1965.                                                   lineNumber, nameLove, line);
  1966.                     }
  1967.                     final int n = Integer.parseInt(fields[0]);
  1968.                     final int m = Integer.parseInt(fields[1]);
  1969.                     if ((n < 2) || (n > 3) || (m < 0) || (m > n)) {
  1970.                         // this should never happen with files embedded within Orekit
  1971.                         throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  1972.                                                   lineNumber, nameLove, line);

  1974.                     }
  1975.                     real[n][m]      = Double.parseDouble(fields[2]);
  1976.                     imaginary[n][m] = Double.parseDouble(fields[3]);
  1977.                     plus[n][m]      = Double.parseDouble(fields[4]);
  1978.                 }

  1980.                 // next line
  1981.                 lineNumber++;
  1982.                 line = reader.readLine();

  1984.             }

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

  1988.         } catch (IOException ioe) {
  1989.             throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, nameLove);
  1990.         } finally {
  1991.             try {
  1992.                 if (stream != null) {
  1993.                     stream.close();
  1994.                 }
  1995.                 if (reader != null) {
  1996.                     reader.close();
  1997.                 }
  1998.             } catch (IOException ioe) {
  1999.                 // ignored here
  2000.             }
  2001.         }
  2002.     }

  2004.     /** Get a data stream.
  2005.      * @param name file name of the resource stream
  2006.      * @return stream
  2007.      */
  2008.     private static InputStream getStream(final String name) {
  2009.         return IERSConventions.class.getResourceAsStream(name);
  2010.     }

  2012.     /** Correction to equation of equinoxes.
  2013.      * <p>IAU 1994 resolution C7 added two terms to the equation of equinoxes
  2014.      * taking effect since 1997-02-27 for continuity
  2015.      * </p>
  2016.      */
  2017.     private static class IAU1994ResolutionC7 {

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

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

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

  2031.         /** Evaluate the correction.
  2032.          * @param arguments Delaunay for nutation
  2033.          * @return correction value (0 before 1997-02-27)
  2034.          */
  2035.         public static double value(final DelaunayArguments arguments) {
  2036.             if (arguments.getDate().compareTo(MODEL_START) >= 0) {

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

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

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

  2047.             } else {
  2048.                 return 0.0;
  2049.             }
  2050.         }

  2052.     };

  2054.     /** Stellar angle model.
  2055.      * <p>
  2056.      * The stellar angle computed here has been defined in the paper "A non-rotating origin on the
  2057.      * instantaneous equator: Definition, properties and use", N. Capitaine, Guinot B., and Souchay J.,
  2058.      * Celestial Mechanics, Volume 39, Issue 3, pp 283-307. It has been proposed as a conventional
  2059.      * conventional relationship between UT1 and Earth rotation in the paper "Definition of the Celestial
  2060.      * Ephemeris origin and of UT1 in the International Celestial Reference Frame", Capitaine, N.,
  2061.      * Guinot, B., and McCarthy, D. D., 2000, “,” Astronomy and Astrophysics, 355(1), pp. 398–405.
  2062.      * </p>
  2063.      * <p>
  2064.      * It is presented simply as stellar angle in IERS conventions 1996 but as since been adopted as
  2065.      * the conventional relationship defining UT1 from ICRF and is called Earth Rotation Angle in
  2066.      * IERS conventions 2003 and 2010.
  2067.      * </p>
  2068.      */
  2069.     private static class StellarAngleCapitaine implements TimeFunction<DerivativeStructure> {

  2071.         /** Reference date of Capitaine's Earth Rotation Angle model. */
  2072.         private static final AbsoluteDate REFERENCE_DATE = new AbsoluteDate(DateComponents.J2000_EPOCH,
  2073.                                                                             TimeComponents.H12,
  2074.                                                                             TimeScalesFactory.getTAI());

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

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

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

  2087.         /** Total rate term of Capitaine's Earth Rotation Angle model. */
  2088.         private static final double ERA_1AB = ERA_1A + ERA_1B;

  2090.         /** UT1 time scale. */
  2091.         private final TimeScale ut1;

  2093.         /** Simple constructor.
  2094.          * @param ut1 UT1 time scale
  2095.          */
  2096.         public StellarAngleCapitaine(final TimeScale ut1) {
  2097.             this.ut1 = ut1;
  2098.         }

  2100.         /** {@inheritDoc} */
  2101.         @Override
  2102.         public DerivativeStructure value(final AbsoluteDate date) {

  2104.             // split the date offset as a full number of days plus a smaller part
  2105.             final int secondsInDay = 86400;
  2106.             final double dt  = date.durationFrom(REFERENCE_DATE);
  2107.             final long days  = ((long) dt) / secondsInDay;
  2108.             final double dtA = secondsInDay * days;
  2109.             final double dtB = (dt - dtA) + ut1.offsetFromTAI(date);

  2111.             return new DerivativeStructure(1, 1,
  2112.                                            ERA_0 + ERA_1A * dtB + ERA_1B * (dtA + dtB),
  2113.                                            ERA_1AB);

  2115.         }

  2117.     }

  2119.     /** Mean pole. */
  2120.     private static class MeanPole implements TimeStamped, Serializable {

  2122.         /** Serializable UID. */
  2123.         private static final long serialVersionUID = 20131028l;

  2125.         /** Date. */
  2126.         private final AbsoluteDate date;

  2128.         /** X coordinate. */
  2129.         private double x;

  2131.         /** Y coordinate. */
  2132.         private double y;

  2134.         /** Simple constructor.
  2135.          * @param date date
  2136.          * @param x x coordinate
  2137.          * @param y y coordinate
  2138.          */
  2139.         public MeanPole(final AbsoluteDate date, final double x, final double y) {
  2140.             this.date = date;
  2141.             this.x    = x;
  2142.             this.y    = y;
  2143.         }

  2145.         /** {@inheritDoc} */
  2146.         @Override
  2147.         public AbsoluteDate getDate() {
  2148.             return date;
  2149.         }

  2151.         /** Get x coordinate.
  2152.          * @return x coordinate
  2153.          */
  2154.         public double getX() {
  2155.             return x;
  2156.         }

  2158.         /** Get y coordinate.
  2159.          * @return y coordinate
  2160.          */
  2161.         public double getY() {
  2162.             return y;
  2163.         }

  2165.     }
  2166. }
Created with JaCoCo 0.6.3.201306030806