IERSConventions.java
/* Copyright 2002-2024 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.utils;
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.io.Serializable;
import java.nio.charset.StandardCharsets;
import java.util.List;
import java.util.function.Function;
import java.util.regex.Pattern;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.analysis.interpolation.FieldHermiteInterpolator;
import org.hipparchus.analysis.interpolation.HermiteInterpolator;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.FieldSinCos;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.MathUtils;
import org.hipparchus.util.SinCos;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.data.BodiesElements;
import org.orekit.data.DataContext;
import org.orekit.data.DelaunayArguments;
import org.orekit.data.FieldBodiesElements;
import org.orekit.data.FieldDelaunayArguments;
import org.orekit.data.FundamentalNutationArguments;
import org.orekit.data.PoissonSeries;
import org.orekit.data.PoissonSeries.CompiledSeries;
import org.orekit.data.PoissonSeriesParser;
import org.orekit.data.PolynomialNutation;
import org.orekit.data.PolynomialParser;
import org.orekit.data.PolynomialParser.Unit;
import org.orekit.data.SimpleTimeStampedTableParser;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitInternalError;
import org.orekit.errors.OrekitMessages;
import org.orekit.errors.TimeStampedCacheException;
import org.orekit.frames.EOPHistory;
import org.orekit.frames.FieldPoleCorrection;
import org.orekit.frames.PoleCorrection;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateComponents;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeComponents;
import org.orekit.time.TimeScalarFunction;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScales;
import org.orekit.time.TimeStamped;
import org.orekit.time.TimeVectorFunction;
/** Supported IERS conventions.
* @since 6.0
* @author Luc Maisonobe
*/
public enum IERSConventions {
/** Constant for IERS 1996 conventions. */
IERS_1996 {
/** Nutation arguments resources. */
private static final String NUTATION_ARGUMENTS = IERS_BASE + "1996/nutation-arguments.txt";
/** X series resources. */
private static final String X_Y_SERIES = IERS_BASE + "1996/tab5.4.txt";
/** Psi series resources. */
private static final String PSI_EPSILON_SERIES = IERS_BASE + "1996/tab5.1.txt";
/** Tidal correction for xp, yp series resources. */
private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "1996/tab8.4.txt";
/** Tidal correction for UT1 resources. */
private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "1996/tab8.3.txt";
/** Love numbers resources. */
private static final String LOVE_NUMBERS = IERS_BASE + "1996/tab6.1.txt";
/** Frequency dependence model for k₂₀. */
private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2b.txt";
/** Frequency dependence model for k₂₁. */
private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2a.txt";
/** Frequency dependence model for k₂₂. */
private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "1996/tab6.2c.txt";
/** Tidal displacement frequency correction for diurnal tides. */
private static final String TIDAL_DISPLACEMENT_CORRECTION_DIURNAL = IERS_BASE + "1996/tab7.3a.txt";
/** Tidal displacement frequency correction for zonal tides. */
private static final String TIDAL_DISPLACEMENT_CORRECTION_ZONAL = IERS_BASE + "1996/tab7.3b.txt";
/** {@inheritDoc} */
@Override
public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
final TimeScales timeScales) {
return load(NUTATION_ARGUMENTS, in -> new FundamentalNutationArguments(this, timeScale,
in, NUTATION_ARGUMENTS,
timeScales));
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {
// value from chapter 5, page 22
final PolynomialNutation epsilonA =
new PolynomialNutation(84381.448 * Constants.ARC_SECONDS_TO_RADIANS,
-46.8150 * Constants.ARC_SECONDS_TO_RADIANS,
-0.00059 * Constants.ARC_SECONDS_TO_RADIANS,
0.001813 * Constants.ARC_SECONDS_TO_RADIANS);
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
return epsilonA.value(evaluateTC(date, timeScales));
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
return epsilonA.value(evaluateTC(date, timeScales));
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {
// set up nutation arguments
final FundamentalNutationArguments arguments =
getNutationArguments(timeScales);
// X = 2004.3109″t - 0.42665″t² - 0.198656″t³ + 0.0000140″t⁴
// + 0.00006″t² cos Ω + sin ε0 { Σ [(Ai + Ai' t) sin(ARGUMENT) + Ai'' t cos(ARGUMENT)]}
// + 0.00204″t² sin Ω + 0.00016″t² sin 2(F - D + Ω),
final PolynomialNutation xPolynomial =
new PolynomialNutation(0,
2004.3109 * Constants.ARC_SECONDS_TO_RADIANS,
-0.42665 * Constants.ARC_SECONDS_TO_RADIANS,
-0.198656 * Constants.ARC_SECONDS_TO_RADIANS,
0.0000140 * Constants.ARC_SECONDS_TO_RADIANS);
final double fXCosOm = 0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
final double fXSinOm = 0.00204 * Constants.ARC_SECONDS_TO_RADIANS;
final double fXSin2FDOm = 0.00016 * Constants.ARC_SECONDS_TO_RADIANS;
final double sinEps0 = FastMath.sin(getMeanObliquityFunction(timeScales)
.value(getNutationReferenceEpoch(timeScales)));
final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
final PoissonSeriesParser baseParser =
new PoissonSeriesParser(12).withFirstDelaunay(1);
final PoissonSeriesParser xParser =
baseParser.
withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
withSinCos(1, 8, deciMilliAS, 9, deciMilliAS);
final PoissonSeries xSum = load(X_Y_SERIES, in -> xParser.parse(in, X_Y_SERIES));
// Y = -0.00013″ - 22.40992″t² + 0.001836″t³ + 0.0011130″t⁴
// + Σ [(Bi + Bi' t) cos(ARGUMENT) + Bi'' t sin(ARGUMENT)]
// - 0.00231″t² cos Ω − 0.00014″t² cos 2(F - D + Ω)
final PolynomialNutation yPolynomial =
new PolynomialNutation(-0.00013 * Constants.ARC_SECONDS_TO_RADIANS,
0.0,
-22.40992 * Constants.ARC_SECONDS_TO_RADIANS,
0.001836 * Constants.ARC_SECONDS_TO_RADIANS,
0.0011130 * Constants.ARC_SECONDS_TO_RADIANS);
final double fYCosOm = -0.00231 * Constants.ARC_SECONDS_TO_RADIANS;
final double fYCos2FDOm = -0.00014 * Constants.ARC_SECONDS_TO_RADIANS;
final PoissonSeriesParser yParser =
baseParser.
withSinCos(0, -1, deciMilliAS, 10, deciMilliAS).
withSinCos(1, 12, deciMilliAS, 11, deciMilliAS);
final PoissonSeries ySum = load(X_Y_SERIES, in -> yParser.parse(in, X_Y_SERIES));
final PoissonSeries.CompiledSeries xySum =
PoissonSeries.compile(xSum, ySum);
// s = -XY/2 + 0.00385″t - 0.07259″t³ - 0.00264″ sin Ω - 0.00006″ sin 2Ω
// + 0.00074″t² sin Ω + 0.00006″t² sin 2(F - D + Ω)
final double fST = 0.00385 * Constants.ARC_SECONDS_TO_RADIANS;
final double fST3 = -0.07259 * Constants.ARC_SECONDS_TO_RADIANS;
final double fSSinOm = -0.00264 * Constants.ARC_SECONDS_TO_RADIANS;
final double fSSin2Om = -0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
final double fST2SinOm = 0.00074 * Constants.ARC_SECONDS_TO_RADIANS;
final double fST2Sin2FDOm = 0.00006 * Constants.ARC_SECONDS_TO_RADIANS;
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
final BodiesElements elements = arguments.evaluateAll(date);
final double[] xy = xySum.value(elements);
final double omega = elements.getOmega();
final double f = elements.getF();
final double d = elements.getD();
final double t = elements.getTC();
final SinCos scOmega = FastMath.sinCos(omega);
final SinCos sc2omega = SinCos.sum(scOmega, scOmega);
final SinCos sc2FD0m = FastMath.sinCos(2 * (f - d + omega));
final double cosOmega = scOmega.cos();
final double sinOmega = scOmega.sin();
final double sin2Omega = sc2omega.sin();
final double cos2FDOm = sc2FD0m.cos();
final double sin2FDOm = sc2FD0m.sin();
final double x = xPolynomial.value(t) + sinEps0 * xy[0] +
t * t * (fXCosOm * cosOmega + fXSinOm * sinOmega + fXSin2FDOm * cos2FDOm);
final double y = yPolynomial.value(t) + xy[1] +
t * t * (fYCosOm * cosOmega + fYCos2FDOm * cos2FDOm);
final double sPxy2 = fSSinOm * sinOmega + fSSin2Om * sin2Omega +
t * (fST + t * (fST2SinOm * sinOmega + fST2Sin2FDOm * sin2FDOm + t * fST3));
return new double[] {
x, y, sPxy2
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
final T[] xy = xySum.value(elements);
final T omega = elements.getOmega();
final T f = elements.getF();
final T d = elements.getD();
final T t = elements.getTC();
final T t2 = t.square();
final FieldSinCos<T> scOmega = FastMath.sinCos(omega);
final FieldSinCos<T> sc2omega = FieldSinCos.sum(scOmega, scOmega);
final FieldSinCos<T> sc2FD0m = FastMath.sinCos(f.subtract(d).add(omega).multiply(2));
final T cosOmega = scOmega.cos();
final T sinOmega = scOmega.sin();
final T sin2Omega = sc2omega.sin();
final T cos2FDOm = sc2FD0m.cos();
final T sin2FDOm = sc2FD0m.sin();
final T x = xPolynomial.value(t).
add(xy[0].multiply(sinEps0)).
add(t2.multiply(cosOmega.multiply(fXCosOm).add(sinOmega.multiply(fXSinOm)).add(cos2FDOm.multiply(fXSin2FDOm))));
final T y = yPolynomial.value(t).
add(xy[1]).
add(t2.multiply(cosOmega.multiply(fYCosOm).add(cos2FDOm.multiply(fYCos2FDOm))));
final T sPxy2 = sinOmega.multiply(fSSinOm).
add(sin2Omega.multiply(fSSin2Om)).
add(t.multiply(fST3).add(sinOmega.multiply(fST2SinOm)).add(sin2FDOm.multiply(fST2Sin2FDOm)).multiply(t).add(fST).multiply(t));
final T[] a = MathArrays.buildArray(date.getField(), 3);
a[0] = x;
a[1] = y;
a[2] = sPxy2;
return a;
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {
// set up the conventional polynomials
// the following values are from Lieske et al. paper:
// Expressions for the precession quantities based upon the IAU(1976) system of astronomical constants
// http://articles.adsabs.harvard.edu/full/1977A%26A....58....1L
// also available as equation 30 in IERS 2003 conventions
final PolynomialNutation psiA =
new PolynomialNutation( 0.0,
5038.7784 * Constants.ARC_SECONDS_TO_RADIANS,
-1.07259 * Constants.ARC_SECONDS_TO_RADIANS,
-0.001147 * Constants.ARC_SECONDS_TO_RADIANS);
final PolynomialNutation omegaA =
new PolynomialNutation(getMeanObliquityFunction(timeScales)
.value(getNutationReferenceEpoch(timeScales)),
0.0,
0.05127 * Constants.ARC_SECONDS_TO_RADIANS,
-0.007726 * Constants.ARC_SECONDS_TO_RADIANS);
final PolynomialNutation chiA =
new PolynomialNutation( 0.0,
10.5526 * Constants.ARC_SECONDS_TO_RADIANS,
-2.38064 * Constants.ARC_SECONDS_TO_RADIANS,
-0.001125 * Constants.ARC_SECONDS_TO_RADIANS);
return new PrecessionFunction(psiA, omegaA, chiA, timeScales);
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {
// set up nutation arguments
final FundamentalNutationArguments arguments =
getNutationArguments(timeScales);
// set up Poisson series
final double deciMilliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-4;
final PoissonSeriesParser baseParser =
new PoissonSeriesParser(10).withFirstDelaunay(1);
final PoissonSeriesParser psiParser =
baseParser.
withSinCos(0, 7, deciMilliAS, -1, deciMilliAS).
withSinCos(1, 8, deciMilliAS, -1, deciMilliAS);
final PoissonSeries psiSeries = load(PSI_EPSILON_SERIES, in -> psiParser.parse(in, PSI_EPSILON_SERIES));
final PoissonSeriesParser epsilonParser =
baseParser.
withSinCos(0, -1, deciMilliAS, 9, deciMilliAS).
withSinCos(1, -1, deciMilliAS, 10, deciMilliAS);
final PoissonSeries epsilonSeries = load(PSI_EPSILON_SERIES, in -> epsilonParser.parse(in, PSI_EPSILON_SERIES));
final PoissonSeries.CompiledSeries psiEpsilonSeries =
PoissonSeries.compile(psiSeries, epsilonSeries);
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
final BodiesElements elements = arguments.evaluateAll(date);
final double[] psiEpsilon = psiEpsilonSeries.value(elements);
return new double[] {
psiEpsilon[0], psiEpsilon[1], IAU1994ResolutionC7.value(elements, timeScales.getTAI())
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
final T[] psiEpsilon = psiEpsilonSeries.value(elements);
final T[] result = MathArrays.buildArray(date.getField(), 3);
result[0] = psiEpsilon[0];
result[1] = psiEpsilon[1];
result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
return result;
}
};
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
final TimeScales timeScales) {
// Radians per second of time
final double radiansPerSecond = MathUtils.TWO_PI / Constants.JULIAN_DAY;
// constants from IERS 1996 page 21
// the underlying model is IAU 1982 GMST-UT1
final AbsoluteDate gmstReference = new AbsoluteDate(
DateComponents.J2000_EPOCH, TimeComponents.H12, timeScales.getTAI());
final double gmst0 = 24110.54841;
final double gmst1 = 8640184.812866;
final double gmst2 = 0.093104;
final double gmst3 = -6.2e-6;
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
// offset in Julian centuries from J2000 epoch (UT1 scale)
final double dtai = date.durationFrom(gmstReference);
final double tut1 = dtai + ut1.offsetFromTAI(date);
final double tt = tut1 / Constants.JULIAN_CENTURY;
// Seconds in the day, adjusted by 12 hours because the
// UT1 is supplied as a Julian date beginning at noon.
final double sd = FastMath.IEEEremainder(tut1 + Constants.JULIAN_DAY / 2, Constants.JULIAN_DAY);
// compute Greenwich mean sidereal time, in radians
return ((((((tt * gmst3 + gmst2) * tt) + gmst1) * tt) + gmst0) + sd) * radiansPerSecond;
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
// offset in Julian centuries from J2000 epoch (UT1 scale)
final T dtai = date.durationFrom(gmstReference);
final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
final T tt = tut1.divide(Constants.JULIAN_CENTURY);
// Seconds in the day, adjusted by 12 hours because the
// UT1 is supplied as a Julian date beginning at noon.
final T sd = tut1.add(Constants.JULIAN_DAY / 2).remainder(Constants.JULIAN_DAY);
// compute Greenwich mean sidereal time, in radians
return tt.multiply(gmst3).add(gmst2).multiply(tt).add(gmst1).multiply(tt).add(gmst0).add(sd).multiply(radiansPerSecond);
}
};
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
final TimeScales timeScales) {
// Radians per second of time
final double radiansPerSecond = MathUtils.TWO_PI / Constants.JULIAN_DAY;
// constants from IERS 1996 page 21
// the underlying model is IAU 1982 GMST-UT1
final AbsoluteDate gmstReference = new AbsoluteDate(
DateComponents.J2000_EPOCH, TimeComponents.H12, timeScales.getTAI());
final double gmst1 = 8640184.812866;
final double gmst2 = 0.093104;
final double gmst3 = -6.2e-6;
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
// offset in Julian centuries from J2000 epoch (UT1 scale)
final double dtai = date.durationFrom(gmstReference);
final double tut1 = dtai + ut1.offsetFromTAI(date);
final double tt = tut1 / Constants.JULIAN_CENTURY;
// compute Greenwich mean sidereal time rate, in radians per second
return ((((tt * 3 * gmst3 + 2 * gmst2) * tt) + gmst1) / Constants.JULIAN_CENTURY + 1) * radiansPerSecond;
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
// offset in Julian centuries from J2000 epoch (UT1 scale)
final T dtai = date.durationFrom(gmstReference);
final T tut1 = dtai.add(ut1.offsetFromTAI(date.toAbsoluteDate()));
final T tt = tut1.divide(Constants.JULIAN_CENTURY);
// compute Greenwich mean sidereal time, in radians
return tt.multiply(3 * gmst3).add(2 * gmst2).multiply(tt).add(gmst1).divide(Constants.JULIAN_CENTURY).add(1).multiply(radiansPerSecond);
}
};
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getGASTFunction(final TimeScale ut1,
final EOPHistory eopHistory,
final TimeScales timeScales) {
// obliquity
final TimeScalarFunction epsilonA = getMeanObliquityFunction(timeScales);
// GMST function
final TimeScalarFunction gmst = getGMSTFunction(ut1, timeScales);
// nutation function
final TimeVectorFunction nutation = getNutationFunction(timeScales);
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
// compute equation of equinoxes
final double[] angles = nutation.value(date);
double deltaPsi = angles[0];
if (eopHistory != null) {
deltaPsi += eopHistory.getEquinoxNutationCorrection(date)[0];
}
final double eqe = deltaPsi * FastMath.cos(epsilonA.value(date)) + angles[2];
// add mean sidereal time and equation of equinoxes
return gmst.value(date) + eqe;
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
// compute equation of equinoxes
final T[] angles = nutation.value(date);
T deltaPsi = angles[0];
if (eopHistory != null) {
deltaPsi = deltaPsi.add(eopHistory.getEquinoxNutationCorrection(date)[0]);
}
final T eqe = deltaPsi.multiply(epsilonA.value(date).cos()).add(angles[2]);
// add mean sidereal time and equation of equinoxes
return gmst.value(date).add(eqe);
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {
// set up nutation arguments
// BEWARE! Using TT as the time scale here and not UT1 is intentional!
// as this correction is used to compute UT1 itself, it is not surprising we cannot use UT1 yet,
// however, using the close UTC as would seem logical make the comparison with interp.f from IERS fail
// looking in the interp.f code, the same TT scale is used for both Delaunay and gamma argument
final FundamentalNutationArguments arguments =
getNutationArguments(timeScales.getTT(), timeScales);
// set up Poisson series
final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
final PoissonSeriesParser xyParser = new PoissonSeriesParser(17).
withOptionalColumn(1).
withGamma(7).
withFirstDelaunay(2);
final PoissonSeries xSeries =
load(TIDAL_CORRECTION_XP_YP_SERIES, in -> xyParser.
withSinCos(0, 14, milliAS, 15, milliAS).
parse(in, TIDAL_CORRECTION_XP_YP_SERIES));
final PoissonSeries ySeries =
load(TIDAL_CORRECTION_XP_YP_SERIES, in -> xyParser.
withSinCos(0, 16, milliAS, 17, milliAS).
parse(in, TIDAL_CORRECTION_XP_YP_SERIES));
final double deciMilliS = 1.0e-4;
final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(17).
withOptionalColumn(1).
withGamma(7).
withFirstDelaunay(2);
final PoissonSeries ut1Series =
load(TIDAL_CORRECTION_UT1_SERIES, in -> ut1Parser.
withSinCos(0, 16, deciMilliS, 17, deciMilliS).
parse(in, TIDAL_CORRECTION_UT1_SERIES));
return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
}
/** {@inheritDoc} */
public LoveNumbers getLoveNumbers() {
return loadLoveNumbers(LOVE_NUMBERS);
}
/** {@inheritDoc} */
public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
final TimeScales timeScales) {
// set up nutation arguments
final FundamentalNutationArguments arguments = getNutationArguments(ut1, timeScales);
// set up Poisson series
final PoissonSeriesParser k20Parser =
new PoissonSeriesParser(18).
withOptionalColumn(1).
withDoodson(4, 2).
withFirstDelaunay(10);
final PoissonSeriesParser k21Parser =
new PoissonSeriesParser(18).
withOptionalColumn(1).
withDoodson(4, 3).
withFirstDelaunay(10);
final PoissonSeriesParser k22Parser =
new PoissonSeriesParser(16).
withOptionalColumn(1).
withDoodson(4, 2).
withFirstDelaunay(10);
final double pico = 1.0e-12;
final PoissonSeries c20Series =
load(K20_FREQUENCY_DEPENDENCE, in -> k20Parser.
withSinCos(0, 18, -pico, 16, pico).
parse(in, K20_FREQUENCY_DEPENDENCE));
final PoissonSeries c21Series =
load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
withSinCos(0, 17, pico, 18, pico).
parse(in, K21_FREQUENCY_DEPENDENCE));
final PoissonSeries s21Series =
load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
withSinCos(0, 18, -pico, 17, pico).
parse(in, K21_FREQUENCY_DEPENDENCE));
final PoissonSeries c22Series =
load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
withSinCos(0, -1, pico, 16, pico).
parse(in, K22_FREQUENCY_DEPENDENCE));
final PoissonSeries s22Series =
load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
withSinCos(0, 16, -pico, -1, pico).
parse(in, K22_FREQUENCY_DEPENDENCE));
return new TideFrequencyDependenceFunction(arguments,
c20Series,
c21Series, s21Series,
c22Series, s22Series);
}
/** {@inheritDoc} */
@Override
public double getPermanentTide() {
return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {
// constants from IERS 1996 page 47
final double globalFactor = -1.348e-9 / Constants.ARC_SECONDS_TO_RADIANS;
final double coupling = 0.0112;
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
final PoleCorrection pole = eopHistory.getPoleCorrection(date);
return new double[] {
globalFactor * (pole.getXp() + coupling * pole.getYp()),
globalFactor * (coupling * pole.getXp() - pole.getYp()),
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
final FieldPoleCorrection<T> pole = eopHistory.getPoleCorrection(date);
final T[] a = MathArrays.buildArray(date.getField(), 2);
a[0] = pole.getXp().add(pole.getYp().multiply(coupling)).multiply(globalFactor);
a[1] = pole.getXp().multiply(coupling).subtract(pole.getYp()).multiply(globalFactor);
return a;
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
// there are no model for ocean pole tide prior to conventions 2010
return new double[] {
0.0, 0.0
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
// there are no model for ocean pole tide prior to conventions 2010
return MathArrays.buildArray(date.getField(), 2);
}
};
}
/** {@inheritDoc} */
@Override
public double[] getNominalTidalDisplacement() {
// // elastic Earth values
// return new double[] {
// // h⁽⁰⁾, h⁽²⁾, h₃, hI diurnal, hI semi-diurnal,
// 0.6026, -0.0006, 0.292, -0.0025, -0.0022,
// // l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾, l₃, lI diurnal, lI semi-diurnal
// 0.0831, 0.0012, 0.0024, 0.0002, 0.015, -0.0007, -0.0007,
// // H₀
// -0.31460
// };
// anelastic Earth values
return new double[] {
// h⁽⁰⁾, h⁽²⁾, h₃, hI diurnal, hI semi-diurnal,
0.6078, -0.0006, 0.292, -0.0025, -0.0022,
// l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾, l₃, lI diurnal, lI semi-diurnal
0.0847, 0.0012, 0.0024, 0.0002, 0.015, -0.0007, -0.0007,
// H₀
-0.31460
};
}
/** {@inheritDoc} */
@Override
public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
18, 17, -1, 18, -1);
}
/** {@inheritDoc} */
@Override
public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
20, 17, 19, 18, 20);
}
},
/** Constant for IERS 2003 conventions. */
IERS_2003 {
/** Nutation arguments resources. */
private static final String NUTATION_ARGUMENTS = IERS_BASE + "2003/nutation-arguments.txt";
/** X series resources. */
private static final String X_SERIES = IERS_BASE + "2003/tab5.2a.txt";
/** Y series resources. */
private static final String Y_SERIES = IERS_BASE + "2003/tab5.2b.txt";
/** S series resources. */
private static final String S_SERIES = IERS_BASE + "2003/tab5.2c.txt";
/** Luni-solar series resources. */
private static final String LUNI_SOLAR_SERIES = IERS_BASE + "2003/tab5.3a-first-table.txt";
/** Planetary series resources. */
private static final String PLANETARY_SERIES = IERS_BASE + "2003/tab5.3b.txt";
/** Greenwhich sidereal time series resources. */
private static final String GST_SERIES = IERS_BASE + "2003/tab5.4.txt";
/** Tidal correction for xp, yp series resources. */
private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "2003/tab8.2ab.txt";
/** Tidal correction for UT1 resources. */
private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "2003/tab8.3ab.txt";
/** Love numbers resources. */
private static final String LOVE_NUMBERS = IERS_BASE + "2003/tab6.1.txt";
/** Frequency dependence model for k₂₀. */
private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3b.txt";
/** Frequency dependence model for k₂₁. */
private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3a.txt";
/** Frequency dependence model for k₂₂. */
private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "2003/tab6.3c.txt";
/** Annual pole table. */
private static final String ANNUAL_POLE = IERS_BASE + "2003/annual.pole";
/** Tidal displacement frequency correction for diurnal tides. */
private static final String TIDAL_DISPLACEMENT_CORRECTION_DIURNAL = IERS_BASE + "2003/tab7.5a.txt";
/** Tidal displacement frequency correction for zonal tides. */
private static final String TIDAL_DISPLACEMENT_CORRECTION_ZONAL = IERS_BASE + "2003/tab7.5b.txt";
/** {@inheritDoc} */
public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
final TimeScales timeScales) {
return load(NUTATION_ARGUMENTS, in -> new FundamentalNutationArguments(this, timeScale,
in, NUTATION_ARGUMENTS,
timeScales));
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {
// epsilon 0 value from chapter 5, page 41, other terms from equation 32 page 45
final PolynomialNutation epsilonA =
new PolynomialNutation(84381.448 * Constants.ARC_SECONDS_TO_RADIANS,
-46.84024 * Constants.ARC_SECONDS_TO_RADIANS,
-0.00059 * Constants.ARC_SECONDS_TO_RADIANS,
0.001813 * Constants.ARC_SECONDS_TO_RADIANS);
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
return epsilonA.value(evaluateTC(date, timeScales));
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
return epsilonA.value(evaluateTC(date, timeScales));
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {
// set up nutation arguments
final FundamentalNutationArguments arguments =
getNutationArguments(timeScales);
// set up Poisson series
final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
final PoissonSeriesParser parser =
new PoissonSeriesParser(17).
withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
withFirstDelaunay(4).
withFirstPlanetary(9).
withSinCos(0, 2, microAS, 3, microAS);
final PoissonSeries xSeries = load(X_SERIES, in -> parser.parse(in, X_SERIES));
final PoissonSeries ySeries = load(Y_SERIES, in -> parser.parse(in, Y_SERIES));
final PoissonSeries sSeries = load(S_SERIES, in -> parser.parse(in, S_SERIES));
final PoissonSeries.CompiledSeries xys = PoissonSeries.compile(xSeries, ySeries, sSeries);
// create a function evaluating the series
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
return xys.value(arguments.evaluateAll(date));
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
return xys.value(arguments.evaluateAll(date));
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {
// set up the conventional polynomials
// the following values are from equation 32 in IERS 2003 conventions
final PolynomialNutation psiA =
new PolynomialNutation( 0.0,
5038.47875 * Constants.ARC_SECONDS_TO_RADIANS,
-1.07259 * Constants.ARC_SECONDS_TO_RADIANS,
-0.001147 * Constants.ARC_SECONDS_TO_RADIANS);
final PolynomialNutation omegaA =
new PolynomialNutation(getMeanObliquityFunction(timeScales)
.value(getNutationReferenceEpoch(timeScales)),
-0.02524 * Constants.ARC_SECONDS_TO_RADIANS,
0.05127 * Constants.ARC_SECONDS_TO_RADIANS,
-0.007726 * Constants.ARC_SECONDS_TO_RADIANS);
final PolynomialNutation chiA =
new PolynomialNutation( 0.0,
10.5526 * Constants.ARC_SECONDS_TO_RADIANS,
-2.38064 * Constants.ARC_SECONDS_TO_RADIANS,
-0.001125 * Constants.ARC_SECONDS_TO_RADIANS);
return new PrecessionFunction(psiA, omegaA, chiA, timeScales);
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {
// set up nutation arguments
final FundamentalNutationArguments arguments =
getNutationArguments(timeScales);
// set up Poisson series
final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
final PoissonSeriesParser luniSolarParser =
new PoissonSeriesParser(14).withFirstDelaunay(1);
final PoissonSeriesParser luniSolarPsiParser =
luniSolarParser.
withSinCos(0, 7, milliAS, 11, milliAS).
withSinCos(1, 8, milliAS, 12, milliAS);
final PoissonSeries psiLuniSolarSeries =
load(LUNI_SOLAR_SERIES, in -> luniSolarPsiParser.parse(in, LUNI_SOLAR_SERIES));
final PoissonSeriesParser luniSolarEpsilonParser = luniSolarParser.
withSinCos(0, 13, milliAS, 9, milliAS).
withSinCos(1, 14, milliAS, 10, milliAS);
final PoissonSeries epsilonLuniSolarSeries =
load(LUNI_SOLAR_SERIES, in -> luniSolarEpsilonParser.parse(in, LUNI_SOLAR_SERIES));
final PoissonSeriesParser planetaryParser =
new PoissonSeriesParser(21).
withFirstDelaunay(2).
withFirstPlanetary(7);
final PoissonSeriesParser planetaryPsiParser =
planetaryParser.withSinCos(0, 17, milliAS, 18, milliAS);
final PoissonSeries psiPlanetarySeries =
load(PLANETARY_SERIES, in -> planetaryPsiParser.parse(in, PLANETARY_SERIES));
final PoissonSeriesParser planetaryEpsilonParser =
planetaryParser.withSinCos(0, 19, milliAS, 20, milliAS);
final PoissonSeries epsilonPlanetarySeries =
load(PLANETARY_SERIES, in -> planetaryEpsilonParser.parse(in, PLANETARY_SERIES));
final PoissonSeries.CompiledSeries luniSolarSeries =
PoissonSeries.compile(psiLuniSolarSeries, epsilonLuniSolarSeries);
final PoissonSeries.CompiledSeries planetarySeries =
PoissonSeries.compile(psiPlanetarySeries, epsilonPlanetarySeries);
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
final BodiesElements elements = arguments.evaluateAll(date);
final double[] luniSolar = luniSolarSeries.value(elements);
final double[] planetary = planetarySeries.value(elements);
return new double[] {
luniSolar[0] + planetary[0], luniSolar[1] + planetary[1],
IAU1994ResolutionC7.value(elements, timeScales.getTAI())
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
final T[] luniSolar = luniSolarSeries.value(elements);
final T[] planetary = planetarySeries.value(elements);
final T[] result = MathArrays.buildArray(date.getField(), 3);
result[0] = luniSolar[0].add(planetary[0]);
result[1] = luniSolar[1].add(planetary[1]);
result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
return result;
}
};
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
final TimeScales timeScales) {
// Earth Rotation Angle
final StellarAngleCapitaine era =
new StellarAngleCapitaine(ut1, timeScales.getTAI());
// Polynomial part of the apparent sidereal time series
// which is the opposite of Equation of Origins (EO)
final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
final PoissonSeriesParser parser =
new PoissonSeriesParser(17).
withFirstDelaunay(4).
withFirstPlanetary(9).
withSinCos(0, 2, microAS, 3, microAS).
withPolynomialPart('t', Unit.ARC_SECONDS);
final PolynomialNutation minusEO = load(GST_SERIES, in -> parser.parse(in, GST_SERIES).getPolynomial());
// create a function evaluating the series
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
return era.value(date) + minusEO.value(evaluateTC(date, timeScales));
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
return era.value(date).add(minusEO.value(evaluateTC(date, timeScales)));
}
};
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
final TimeScales timeScales) {
// Earth Rotation Angle
final StellarAngleCapitaine era =
new StellarAngleCapitaine(ut1, timeScales.getTAI());
// Polynomial part of the apparent sidereal time series
// which is the opposite of Equation of Origins (EO)
final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
final PoissonSeriesParser parser =
new PoissonSeriesParser(17).
withFirstDelaunay(4).
withFirstPlanetary(9).
withSinCos(0, 2, microAS, 3, microAS).
withPolynomialPart('t', Unit.ARC_SECONDS);
final PolynomialNutation minusEO = load(GST_SERIES, in -> parser.parse(in, GST_SERIES).getPolynomial());
// create a function evaluating the series
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
return era.getRate() + minusEO.derivative(evaluateTC(date, timeScales));
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
return minusEO.derivative(evaluateTC(date, timeScales)).add(era.getRate());
}
};
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getGASTFunction(final TimeScale ut1,
final EOPHistory eopHistory,
final TimeScales timeScales) {
// set up nutation arguments
final FundamentalNutationArguments arguments =
getNutationArguments(timeScales);
// mean obliquity function
final TimeScalarFunction epsilon = getMeanObliquityFunction(timeScales);
// set up Poisson series
final double milliAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-3;
final PoissonSeriesParser luniSolarPsiParser =
new PoissonSeriesParser(14).
withFirstDelaunay(1).
withSinCos(0, 7, milliAS, 11, milliAS).
withSinCos(1, 8, milliAS, 12, milliAS);
final PoissonSeries psiLuniSolarSeries =
load(LUNI_SOLAR_SERIES, in -> luniSolarPsiParser.parse(in, LUNI_SOLAR_SERIES));
final PoissonSeriesParser planetaryPsiParser =
new PoissonSeriesParser(21).
withFirstDelaunay(2).
withFirstPlanetary(7).
withSinCos(0, 17, milliAS, 18, milliAS);
final PoissonSeries psiPlanetarySeries =
load(PLANETARY_SERIES, in -> planetaryPsiParser.parse(in, PLANETARY_SERIES));
final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
final PoissonSeriesParser gstParser =
new PoissonSeriesParser(17).
withFirstDelaunay(4).
withFirstPlanetary(9).
withSinCos(0, 2, microAS, 3, microAS).
withPolynomialPart('t', Unit.ARC_SECONDS);
final PoissonSeries gstSeries = load(GST_SERIES, in -> gstParser.parse(in, GST_SERIES));
final PoissonSeries.CompiledSeries psiGstSeries =
PoissonSeries.compile(psiLuniSolarSeries, psiPlanetarySeries, gstSeries);
// ERA function
final TimeScalarFunction era =
getEarthOrientationAngleFunction(ut1, timeScales.getTAI());
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
// evaluate equation of origins
final BodiesElements elements = arguments.evaluateAll(date);
final double[] angles = psiGstSeries.value(elements);
final double ddPsi = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
final double deltaPsi = angles[0] + angles[1] + ddPsi;
final double epsilonA = epsilon.value(date);
// subtract equation of origin from EA
// (hence add the series above which have the sign included)
return era.value(date) + deltaPsi * FastMath.cos(epsilonA) + angles[2];
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
// evaluate equation of origins
final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
final T[] angles = psiGstSeries.value(elements);
final T ddPsi = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
final T deltaPsi = angles[0].add(angles[1]).add(ddPsi);
final T epsilonA = epsilon.value(date);
// subtract equation of origin from EA
// (hence add the series above which have the sign included)
return era.value(date).add(deltaPsi.multiply(epsilonA.cos())).add(angles[2]);
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {
// set up nutation arguments
// BEWARE! Using TT as the time scale here and not UT1 is intentional!
// as this correction is used to compute UT1 itself, it is not surprising we cannot use UT1 yet,
// however, using the close UTC as would seem logical make the comparison with interp.f from IERS fail
// looking in the interp.f code, the same TT scale is used for both Delaunay and gamma argument
final FundamentalNutationArguments arguments =
getNutationArguments(timeScales.getTT(), timeScales);
// set up Poisson series
final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
withOptionalColumn(1).
withGamma(2).
withFirstDelaunay(3);
final double microS = 1.0e-6;
final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
withOptionalColumn(1).
withGamma(2).
withFirstDelaunay(3);
final PoissonSeries xSeries =
load(TIDAL_CORRECTION_XP_YP_SERIES, in -> xyParser.
withSinCos(0, 10, microAS, 11, microAS).
parse(in, TIDAL_CORRECTION_XP_YP_SERIES));
final PoissonSeries ySeries =
load(TIDAL_CORRECTION_XP_YP_SERIES, in -> xyParser.
withSinCos(0, 12, microAS, 13, microAS).
parse(in, TIDAL_CORRECTION_XP_YP_SERIES));
final PoissonSeries ut1Series =
load(TIDAL_CORRECTION_UT1_SERIES, in -> ut1Parser.
withSinCos(0, 10, microS, 11, microS).
parse(in, TIDAL_CORRECTION_UT1_SERIES));
return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
}
/** {@inheritDoc} */
public LoveNumbers getLoveNumbers() {
return loadLoveNumbers(LOVE_NUMBERS);
}
/** {@inheritDoc} */
public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
final TimeScales timeScales) {
// set up nutation arguments
final FundamentalNutationArguments arguments = getNutationArguments(ut1, timeScales);
// set up Poisson series
final PoissonSeriesParser k20Parser =
new PoissonSeriesParser(18).
withOptionalColumn(1).
withDoodson(4, 2).
withFirstDelaunay(10);
final PoissonSeriesParser k21Parser =
new PoissonSeriesParser(18).
withOptionalColumn(1).
withDoodson(4, 3).
withFirstDelaunay(10);
final PoissonSeriesParser k22Parser =
new PoissonSeriesParser(16).
withOptionalColumn(1).
withDoodson(4, 2).
withFirstDelaunay(10);
final double pico = 1.0e-12;
final PoissonSeries c20Series =
load(K20_FREQUENCY_DEPENDENCE, in -> k20Parser.
withSinCos(0, 18, -pico, 16, pico).
parse(in, K20_FREQUENCY_DEPENDENCE));
final PoissonSeries c21Series =
load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
withSinCos(0, 17, pico, 18, pico).
parse(in, K21_FREQUENCY_DEPENDENCE));
final PoissonSeries s21Series =
load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
withSinCos(0, 18, -pico, 17, pico).
parse(in, K21_FREQUENCY_DEPENDENCE));
final PoissonSeries c22Series =
load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
withSinCos(0, -1, pico, 16, pico).
parse(in, K22_FREQUENCY_DEPENDENCE));
final PoissonSeries s22Series =
load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
withSinCos(0, 16, -pico, -1, pico).
parse(in, K22_FREQUENCY_DEPENDENCE));
return new TideFrequencyDependenceFunction(arguments,
c20Series,
c21Series, s21Series,
c22Series, s22Series);
}
/** {@inheritDoc} */
@Override
public double getPermanentTide() {
return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {
// annual pole from ftp://tai.bipm.org/iers/conv2003/chapter7/annual.pole
final TimeScale utc = eopHistory.getTimeScales().getUTC();
final SimpleTimeStampedTableParser.RowConverter<MeanPole> converter =
new SimpleTimeStampedTableParser.RowConverter<MeanPole>() {
/** {@inheritDoc} */
@Override
public MeanPole convert(final double[] rawFields) {
return new MeanPole(new AbsoluteDate((int) rawFields[0], 1, 1, utc),
rawFields[1] * Constants.ARC_SECONDS_TO_RADIANS,
rawFields[2] * Constants.ARC_SECONDS_TO_RADIANS);
}
};
final SimpleTimeStampedTableParser<MeanPole> parser =
new SimpleTimeStampedTableParser<MeanPole>(3, converter);
final List<MeanPole> annualPoleList = load(ANNUAL_POLE, in -> parser.parse(in, ANNUAL_POLE));
final AbsoluteDate firstAnnualPoleDate = annualPoleList.get(0).getDate();
final AbsoluteDate lastAnnualPoleDate = annualPoleList.get(annualPoleList.size() - 1).getDate();
final ImmutableTimeStampedCache<MeanPole> annualCache =
new ImmutableTimeStampedCache<MeanPole>(2, annualPoleList);
// polynomial extension from IERS 2003, section 7.1.4, equations 23a and 23b
final double xp0 = 0.054 * Constants.ARC_SECONDS_TO_RADIANS;
final double xp0Dot = 0.00083 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
final double yp0 = 0.357 * Constants.ARC_SECONDS_TO_RADIANS;
final double yp0Dot = 0.00395 * Constants.ARC_SECONDS_TO_RADIANS / Constants.JULIAN_YEAR;
// constants from IERS 2003, section 6.2
final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
final double ratio = 0.0115;
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
// we can't compute anything before the range covered by the annual pole file
if (date.compareTo(firstAnnualPoleDate) <= 0) {
return new double[] {
0.0, 0.0
};
}
// evaluate mean pole
double meanPoleX = 0;
double meanPoleY = 0;
if (date.compareTo(lastAnnualPoleDate) <= 0) {
// we are within the range covered by the annual pole file,
// we interpolate within it
try {
final HermiteInterpolator interpolator = new HermiteInterpolator();
annualCache.getNeighbors(date).forEach(neighbor ->
interpolator.addSamplePoint(neighbor.getDate().durationFrom(date),
new double[] {
neighbor.getX(), neighbor.getY()
}));
final double[] interpolated = interpolator.value(0);
meanPoleX = interpolated[0];
meanPoleY = interpolated[1];
} catch (TimeStampedCacheException tsce) {
// this should never happen
throw new OrekitInternalError(tsce);
}
} else {
// we are after the range covered by the annual pole file,
// we use the polynomial extension
final double t = date.durationFrom(
eopHistory.getTimeScales().getJ2000Epoch());
meanPoleX = xp0 + t * xp0Dot;
meanPoleY = yp0 + t * yp0Dot;
}
// evaluate wobble variables
final PoleCorrection correction = eopHistory.getPoleCorrection(date);
final double m1 = correction.getXp() - meanPoleX;
final double m2 = meanPoleY - correction.getYp();
return new double[] {
// the following correspond to the equations published in IERS 2003 conventions,
// section 6.2 page 65. In the publication, the equations read:
// ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
// ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
// However, it seems there are sign errors in these equations, which have
// been fixed in IERS 2010 conventions, section 6.4 page 94. In these newer
// publication, the equations read:
// ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
// ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
// the newer equations seem more consistent with the premises as the
// deformation due to the centrifugal potential has the form:
// −Ω²r²/2 sin 2θ Re [k₂(m₁ − im₂) exp(iλ)] where k₂ is the complex
// number 0.3077 + 0.0036i, so the real part in the previous equation is:
// A[Re(k₂) m₁ + Im(k₂) m₂)] cos λ + A[Re(k₂) m₂ - Im(k₂) m₁] sin λ
// identifying this with ∆C₂₁ cos λ + ∆S₂₁ sin λ we get:
// ∆C₂₁ = A Re(k₂) [m₁ + Im(k₂)/Re(k₂) m₂)]
// ∆S₂₁ = A Re(k₂) [m₂ - Im(k₂)/Re(k₂) m₁)]
// and Im(k₂)/Re(k₂) is very close to +0.0115
// As the equation as written in the IERS 2003 conventions are used in
// legacy systems, we have reproduced this alleged error here (and fixed it in
// the IERS 2010 conventions below) for validation purposes. We don't recommend
// using the IERS 2003 conventions for solid pole tide computation other than
// for validation or reproducibility of legacy applications behavior.
// As solid pole tide is small and as the sign change is on the smallest coefficient,
// the effect is quite small. A test case on a propagated orbit showed a position change
// slightly below 0.4m after a 30 days propagation on a Low Earth Orbit
globalFactor * (m1 - ratio * m2),
globalFactor * (m2 + ratio * m1),
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
final AbsoluteDate aDate = date.toAbsoluteDate();
// we can't compute anything before the range covered by the annual pole file
if (aDate.compareTo(firstAnnualPoleDate) <= 0) {
return MathArrays.buildArray(date.getField(), 2);
}
// evaluate mean pole
T meanPoleX = date.getField().getZero();
T meanPoleY = date.getField().getZero();
if (aDate.compareTo(lastAnnualPoleDate) <= 0) {
// we are within the range covered by the annual pole file,
// we interpolate within it
try {
final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
final T[] y = MathArrays.buildArray(date.getField(), 2);
final T zero = date.getField().getZero();
final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
// for example removing derivatives
// if T was DerivativeStructure
annualCache.getNeighbors(aDate).forEach(neighbor -> {
y[0] = zero.newInstance(neighbor.getX());
y[1] = zero.newInstance(neighbor.getY());
interpolator.addSamplePoint(central.durationFrom(neighbor.getDate()).negate(), y);
});
final T[] interpolated = interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
meanPoleX = interpolated[0];
meanPoleY = interpolated[1];
} catch (TimeStampedCacheException tsce) {
// this should never happen
throw new OrekitInternalError(tsce);
}
} else {
// we are after the range covered by the annual pole file,
// we use the polynomial extension
final T t = date.durationFrom(
eopHistory.getTimeScales().getJ2000Epoch());
meanPoleX = t.multiply(xp0Dot).add(xp0);
meanPoleY = t.multiply(yp0Dot).add(yp0);
}
// evaluate wobble variables
final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
final T m1 = correction.getXp().subtract(meanPoleX);
final T m2 = meanPoleY.subtract(correction.getYp());
final T[] a = MathArrays.buildArray(date.getField(), 2);
// the following correspond to the equations published in IERS 2003 conventions,
// section 6.2 page 65. In the publication, the equations read:
// ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
// ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
// However, it seems there are sign errors in these equations, which have
// been fixed in IERS 2010 conventions, section 6.4 page 94. In these newer
// publication, the equations read:
// ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
// ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
// the newer equations seem more consistent with the premises as the
// deformation due to the centrifugal potential has the form:
// −Ω²r²/2 sin 2θ Re [k₂(m₁ − im₂) exp(iλ)] where k₂ is the complex
// number 0.3077 + 0.0036i, so the real part in the previous equation is:
// A[Re(k₂) m₁ + Im(k₂) m₂)] cos λ + A[Re(k₂) m₂ - Im(k₂) m₁] sin λ
// identifying this with ∆C₂₁ cos λ + ∆S₂₁ sin λ we get:
// ∆C₂₁ = A Re(k₂) [m₁ + Im(k₂)/Re(k₂) m₂)]
// ∆S₂₁ = A Re(k₂) [m₂ - Im(k₂)/Re(k₂) m₁)]
// and Im(k₂)/Re(k₂) is very close to +0.0115
// As the equation as written in the IERS 2003 conventions are used in
// legacy systems, we have reproduced this alleged error here (and fixed it in
// the IERS 2010 conventions below) for validation purposes. We don't recommend
// using the IERS 2003 conventions for solid pole tide computation other than
// for validation or reproducibility of legacy applications behavior.
// As solid pole tide is small and as the sign change is on the smallest coefficient,
// the effect is quite small. A test case on a propagated orbit showed a position change
// slightly below 0.4m after a 30 days propagation on a Low Earth Orbit
a[0] = m1.add(m2.multiply(-ratio)).multiply(globalFactor);
a[1] = m2.add(m1.multiply( ratio)).multiply(globalFactor);
return a;
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
// there are no model for ocean pole tide prior to conventions 2010
return new double[] {
0.0, 0.0
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
// there are no model for ocean pole tide prior to conventions 2010
return MathArrays.buildArray(date.getField(), 2);
}
};
}
/** {@inheritDoc} */
@Override
public double[] getNominalTidalDisplacement() {
return new double[] {
// h⁽⁰⁾, h⁽²⁾, h₃, hI diurnal, hI semi-diurnal,
0.6078, -0.0006, 0.292, -0.0025, -0.0022,
// l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾, l₃, lI diurnal, lI semi-diurnal
0.0847, 0.0012, 0.0024, 0.0002, 0.015, -0.0007, -0.0007,
// H₀
-0.31460
};
}
/** {@inheritDoc} */
@Override
public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
18, 15, 16, 17, 18);
}
/** {@inheritDoc} */
@Override
public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
18, 15, 16, 17, 18);
}
},
/** Constant for IERS 2010 conventions. */
IERS_2010 {
/** Nutation arguments resources. */
private static final String NUTATION_ARGUMENTS = IERS_BASE + "2010/nutation-arguments.txt";
/** X series resources. */
private static final String X_SERIES = IERS_BASE + "2010/tab5.2a.txt";
/** Y series resources. */
private static final String Y_SERIES = IERS_BASE + "2010/tab5.2b.txt";
/** S series resources. */
private static final String S_SERIES = IERS_BASE + "2010/tab5.2d.txt";
/** Psi series resources. */
private static final String PSI_SERIES = IERS_BASE + "2010/tab5.3a.txt";
/** Epsilon series resources. */
private static final String EPSILON_SERIES = IERS_BASE + "2010/tab5.3b.txt";
/** Greenwhich sidereal time series resources. */
private static final String GST_SERIES = IERS_BASE + "2010/tab5.2e.txt";
/** Tidal correction for xp, yp series resources. */
private static final String TIDAL_CORRECTION_XP_YP_SERIES = IERS_BASE + "2010/tab8.2ab.txt";
/** Tidal correction for UT1 resources. */
private static final String TIDAL_CORRECTION_UT1_SERIES = IERS_BASE + "2010/tab8.3ab.txt";
/** Love numbers resources. */
private static final String LOVE_NUMBERS = IERS_BASE + "2010/tab6.3.txt";
/** Frequency dependence model for k₂₀. */
private static final String K20_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5b.txt";
/** Frequency dependence model for k₂₁. */
private static final String K21_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5a.txt";
/** Frequency dependence model for k₂₂. */
private static final String K22_FREQUENCY_DEPENDENCE = IERS_BASE + "2010/tab6.5c.txt";
/** Tidal displacement frequency correction for diurnal tides. */
private static final String TIDAL_DISPLACEMENT_CORRECTION_DIURNAL = IERS_BASE + "2010/tab7.3a.txt";
/** Tidal displacement frequency correction for zonal tides. */
private static final String TIDAL_DISPLACEMENT_CORRECTION_ZONAL = IERS_BASE + "2010/tab7.3b.txt";
/** {@inheritDoc} */
public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale,
final TimeScales timeScales) {
return load(NUTATION_ARGUMENTS, in -> new FundamentalNutationArguments(this, timeScale,
in, NUTATION_ARGUMENTS,
timeScales));
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getMeanObliquityFunction(final TimeScales timeScales) {
// epsilon 0 value from chapter 5, page 56, other terms from equation 5.40 page 65
final PolynomialNutation epsilonA =
new PolynomialNutation(84381.406 * Constants.ARC_SECONDS_TO_RADIANS,
-46.836769 * Constants.ARC_SECONDS_TO_RADIANS,
-0.0001831 * Constants.ARC_SECONDS_TO_RADIANS,
0.00200340 * Constants.ARC_SECONDS_TO_RADIANS,
-0.000000576 * Constants.ARC_SECONDS_TO_RADIANS,
-0.0000000434 * Constants.ARC_SECONDS_TO_RADIANS);
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
return epsilonA.value(evaluateTC(date, timeScales));
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
return epsilonA.value(evaluateTC(date, timeScales));
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getXYSpXY2Function(final TimeScales timeScales) {
// set up nutation arguments
final FundamentalNutationArguments arguments =
getNutationArguments(timeScales);
// set up Poisson series
final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
final PoissonSeriesParser parser =
new PoissonSeriesParser(17).
withPolynomialPart('t', PolynomialParser.Unit.MICRO_ARC_SECONDS).
withFirstDelaunay(4).
withFirstPlanetary(9).
withSinCos(0, 2, microAS, 3, microAS);
final PoissonSeries xSeries = load(X_SERIES, in -> parser.parse(in, X_SERIES));
final PoissonSeries ySeries = load(Y_SERIES, in -> parser.parse(in, Y_SERIES));
final PoissonSeries sSeries = load(S_SERIES, in -> parser.parse(in, S_SERIES));
final PoissonSeries.CompiledSeries xys = PoissonSeries.compile(xSeries, ySeries, sSeries);
// create a function evaluating the series
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
return xys.value(arguments.evaluateAll(date));
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
return xys.value(arguments.evaluateAll(date));
}
};
}
/** {@inheritDoc} */
public LoveNumbers getLoveNumbers() {
return loadLoveNumbers(LOVE_NUMBERS);
}
/** {@inheritDoc} */
public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1,
final TimeScales timeScales) {
// set up nutation arguments
final FundamentalNutationArguments arguments = getNutationArguments(ut1, timeScales);
// set up Poisson series
final PoissonSeriesParser k20Parser =
new PoissonSeriesParser(18).
withOptionalColumn(1).
withDoodson(4, 2).
withFirstDelaunay(10);
final PoissonSeriesParser k21Parser =
new PoissonSeriesParser(18).
withOptionalColumn(1).
withDoodson(4, 3).
withFirstDelaunay(10);
final PoissonSeriesParser k22Parser =
new PoissonSeriesParser(16).
withOptionalColumn(1).
withDoodson(4, 2).
withFirstDelaunay(10);
final double pico = 1.0e-12;
final PoissonSeries c20Series = load(K20_FREQUENCY_DEPENDENCE, in -> k20Parser.
withSinCos(0, 18, -pico, 16, pico).
parse(in, K20_FREQUENCY_DEPENDENCE));
final PoissonSeries c21Series = load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
withSinCos(0, 17, pico, 18, pico).
parse(in, K21_FREQUENCY_DEPENDENCE));
final PoissonSeries s21Series = load(K21_FREQUENCY_DEPENDENCE, in -> k21Parser.
withSinCos(0, 18, -pico, 17, pico).
parse(in, K21_FREQUENCY_DEPENDENCE));
final PoissonSeries c22Series = load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
withSinCos(0, -1, pico, 16, pico).
parse(in, K22_FREQUENCY_DEPENDENCE));
final PoissonSeries s22Series = load(K22_FREQUENCY_DEPENDENCE, in -> k22Parser.
withSinCos(0, 16, -pico, -1, pico).
parse(in, K22_FREQUENCY_DEPENDENCE));
return new TideFrequencyDependenceFunction(arguments,
c20Series,
c21Series, s21Series,
c22Series, s22Series);
}
/** {@inheritDoc} */
@Override
public double getPermanentTide() {
return 4.4228e-8 * -0.31460 * getLoveNumbers().getReal(2, 0);
}
/** Compute pole wobble variables m₁ and m₂.
* @param date current date
* @param eopHistory EOP history
* @return array containing m₁ and m₂
*/
private double[] computePoleWobble(final AbsoluteDate date, final EOPHistory eopHistory) {
// polynomial model from IERS 2010, table 7.7
final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
final double f1 = f0 / Constants.JULIAN_YEAR;
final double f2 = f1 / Constants.JULIAN_YEAR;
final double f3 = f2 / Constants.JULIAN_YEAR;
final AbsoluteDate changeDate =
new AbsoluteDate(2010, 1, 1, eopHistory.getTimeScales().getTT());
// evaluate mean pole
final double[] xPolynomial;
final double[] yPolynomial;
if (date.compareTo(changeDate) <= 0) {
xPolynomial = new double[] {
55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
};
yPolynomial = new double[] {
346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
};
} else {
xPolynomial = new double[] {
23.513 * f0, 7.6141 * f1
};
yPolynomial = new double[] {
358.891 * f0, -0.6287 * f1
};
}
double meanPoleX = 0;
double meanPoleY = 0;
final double t = date.durationFrom(
eopHistory.getTimeScales().getJ2000Epoch());
for (int i = xPolynomial.length - 1; i >= 0; --i) {
meanPoleX = meanPoleX * t + xPolynomial[i];
}
for (int i = yPolynomial.length - 1; i >= 0; --i) {
meanPoleY = meanPoleY * t + yPolynomial[i];
}
// evaluate wobble variables
final PoleCorrection correction = eopHistory.getPoleCorrection(date);
final double m1 = correction.getXp() - meanPoleX;
final double m2 = meanPoleY - correction.getYp();
return new double[] {
m1, m2
};
}
/** Compute pole wobble variables m₁ and m₂.
* @param date current date
* @param <T> type of the field elements
* @param eopHistory EOP history
* @return array containing m₁ and m₂
*/
private <T extends CalculusFieldElement<T>> T[] computePoleWobble(final FieldAbsoluteDate<T> date, final EOPHistory eopHistory) {
// polynomial model from IERS 2010, table 7.7
final double f0 = Constants.ARC_SECONDS_TO_RADIANS / 1000.0;
final double f1 = f0 / Constants.JULIAN_YEAR;
final double f2 = f1 / Constants.JULIAN_YEAR;
final double f3 = f2 / Constants.JULIAN_YEAR;
final AbsoluteDate changeDate =
new AbsoluteDate(2010, 1, 1, eopHistory.getTimeScales().getTT());
// evaluate mean pole
final double[] xPolynomial;
final double[] yPolynomial;
if (date.toAbsoluteDate().compareTo(changeDate) <= 0) {
xPolynomial = new double[] {
55.974 * f0, 1.8243 * f1, 0.18413 * f2, 0.007024 * f3
};
yPolynomial = new double[] {
346.346 * f0, 1.7896 * f1, -0.10729 * f2, -0.000908 * f3
};
} else {
xPolynomial = new double[] {
23.513 * f0, 7.6141 * f1
};
yPolynomial = new double[] {
358.891 * f0, -0.6287 * f1
};
}
T meanPoleX = date.getField().getZero();
T meanPoleY = date.getField().getZero();
final T t = date.durationFrom(
eopHistory.getTimeScales().getJ2000Epoch());
for (int i = xPolynomial.length - 1; i >= 0; --i) {
meanPoleX = meanPoleX.multiply(t).add(xPolynomial[i]);
}
for (int i = yPolynomial.length - 1; i >= 0; --i) {
meanPoleY = meanPoleY.multiply(t).add(yPolynomial[i]);
}
// evaluate wobble variables
final FieldPoleCorrection<T> correction = eopHistory.getPoleCorrection(date);
final T[] m = MathArrays.buildArray(date.getField(), 2);
m[0] = correction.getXp().subtract(meanPoleX);
m[1] = meanPoleY.subtract(correction.getYp());
return m;
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getSolidPoleTide(final EOPHistory eopHistory) {
// constants from IERS 2010, section 6.4
final double globalFactor = -1.333e-9 / Constants.ARC_SECONDS_TO_RADIANS;
final double ratio = 0.0115;
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
// evaluate wobble variables
final double[] wobbleM = computePoleWobble(date, eopHistory);
return new double[] {
// the following correspond to the equations published in IERS 2010 conventions,
// section 6.4 page 94. The equations read:
// ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
// ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
// These equations seem to fix what was probably a sign error in IERS 2003
// conventions section 6.2 page 65. In this older publication, the equations read:
// ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
// ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
globalFactor * (wobbleM[0] + ratio * wobbleM[1]),
globalFactor * (wobbleM[1] - ratio * wobbleM[0])
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
// evaluate wobble variables
final T[] wobbleM = computePoleWobble(date, eopHistory);
final T[] a = MathArrays.buildArray(date.getField(), 2);
// the following correspond to the equations published in IERS 2010 conventions,
// section 6.4 page 94. The equations read:
// ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ + 0.0115m₂)
// ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ − 0.0115m₁)
// These equations seem to fix what was probably a sign error in IERS 2003
// conventions section 6.2 page 65. In this older publication, the equations read:
// ∆C₂₁ = −1.333 × 10⁻⁹ (m₁ − 0.0115m₂)
// ∆S₂₁ = −1.333 × 10⁻⁹ (m₂ + 0.0115m₁)
a[0] = wobbleM[0].add(wobbleM[1].multiply( ratio)).multiply(globalFactor);
a[1] = wobbleM[1].add(wobbleM[0].multiply(-ratio)).multiply(globalFactor);
return a;
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getOceanPoleTide(final EOPHistory eopHistory) {
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
// evaluate wobble variables
final double[] wobbleM = computePoleWobble(date, eopHistory);
return new double[] {
// the following correspond to the equations published in IERS 2010 conventions,
// section 6.4 page 94 equation 6.24:
// ∆C₂₁ = −2.1778 × 10⁻¹⁰ (m₁ − 0.01724m₂)
// ∆S₂₁ = −1.7232 × 10⁻¹⁰ (m₂ − 0.03365m₁)
-2.1778e-10 * (wobbleM[0] - 0.01724 * wobbleM[1]) / Constants.ARC_SECONDS_TO_RADIANS,
-1.7232e-10 * (wobbleM[1] - 0.03365 * wobbleM[0]) / Constants.ARC_SECONDS_TO_RADIANS
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
final T[] v = MathArrays.buildArray(date.getField(), 2);
// evaluate wobble variables
final T[] wobbleM = computePoleWobble(date, eopHistory);
// the following correspond to the equations published in IERS 2010 conventions,
// section 6.4 page 94 equation 6.24:
// ∆C₂₁ = −2.1778 × 10⁻¹⁰ (m₁ − 0.01724m₂)
// ∆S₂₁ = −1.7232 × 10⁻¹⁰ (m₂ − 0.03365m₁)
v[0] = wobbleM[0].subtract(wobbleM[1].multiply(0.01724)).multiply(-2.1778e-10 / Constants.ARC_SECONDS_TO_RADIANS);
v[1] = wobbleM[1].subtract(wobbleM[0].multiply(0.03365)).multiply(-1.7232e-10 / Constants.ARC_SECONDS_TO_RADIANS);
return v;
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getPrecessionFunction(final TimeScales timeScales) {
// set up the conventional polynomials
// the following values are from equation 5.40 in IERS 2010 conventions
final PolynomialNutation psiA =
new PolynomialNutation( 0.0,
5038.481507 * Constants.ARC_SECONDS_TO_RADIANS,
-1.0790069 * Constants.ARC_SECONDS_TO_RADIANS,
-0.00114045 * Constants.ARC_SECONDS_TO_RADIANS,
0.000132851 * Constants.ARC_SECONDS_TO_RADIANS,
-0.0000000951 * Constants.ARC_SECONDS_TO_RADIANS);
final PolynomialNutation omegaA =
new PolynomialNutation(getMeanObliquityFunction(timeScales)
.value(getNutationReferenceEpoch(timeScales)),
-0.025754 * Constants.ARC_SECONDS_TO_RADIANS,
0.0512623 * Constants.ARC_SECONDS_TO_RADIANS,
-0.00772503 * Constants.ARC_SECONDS_TO_RADIANS,
-0.000000467 * Constants.ARC_SECONDS_TO_RADIANS,
0.0000003337 * Constants.ARC_SECONDS_TO_RADIANS);
final PolynomialNutation chiA =
new PolynomialNutation( 0.0,
10.556403 * Constants.ARC_SECONDS_TO_RADIANS,
-2.3814292 * Constants.ARC_SECONDS_TO_RADIANS,
-0.00121197 * Constants.ARC_SECONDS_TO_RADIANS,
0.000170663 * Constants.ARC_SECONDS_TO_RADIANS,
-0.0000000560 * Constants.ARC_SECONDS_TO_RADIANS);
return new PrecessionFunction(psiA, omegaA, chiA, timeScales);
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getNutationFunction(final TimeScales timeScales) {
// set up nutation arguments
final FundamentalNutationArguments arguments =
getNutationArguments(timeScales);
// set up Poisson series
final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
final PoissonSeriesParser parser =
new PoissonSeriesParser(17).
withFirstDelaunay(4).
withFirstPlanetary(9).
withSinCos(0, 2, microAS, 3, microAS);
final PoissonSeries psiSeries = load(PSI_SERIES, in -> parser.parse(in, PSI_SERIES));
final PoissonSeries epsilonSeries = load(EPSILON_SERIES, in -> parser.parse(in, EPSILON_SERIES));
final PoissonSeries.CompiledSeries psiEpsilonSeries = PoissonSeries.compile(psiSeries, epsilonSeries);
return new TimeVectorFunction() {
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
final BodiesElements elements = arguments.evaluateAll(date);
final double[] psiEpsilon = psiEpsilonSeries.value(elements);
return new double[] {
psiEpsilon[0], psiEpsilon[1],
IAU1994ResolutionC7.value(elements, timeScales.getTAI())
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
final T[] psiEpsilon = psiEpsilonSeries.value(elements);
final T[] result = MathArrays.buildArray(date.getField(), 3);
result[0] = psiEpsilon[0];
result[1] = psiEpsilon[1];
result[2] = IAU1994ResolutionC7.value(elements, timeScales.getTAI());
return result;
}
};
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getGMSTFunction(final TimeScale ut1,
final TimeScales timeScales) {
// Earth Rotation Angle
final StellarAngleCapitaine era =
new StellarAngleCapitaine(ut1, timeScales.getTAI());
// Polynomial part of the apparent sidereal time series
// which is the opposite of Equation of Origins (EO)
final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
final PoissonSeriesParser parser =
new PoissonSeriesParser(17).
withFirstDelaunay(4).
withFirstPlanetary(9).
withSinCos(0, 2, microAS, 3, microAS).
withPolynomialPart('t', Unit.ARC_SECONDS);
final PolynomialNutation minusEO = load(GST_SERIES, in -> parser.parse(in, GST_SERIES).getPolynomial());
// create a function evaluating the series
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
return era.value(date) + minusEO.value(evaluateTC(date, timeScales));
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
return era.value(date).add(minusEO.value(evaluateTC(date, timeScales)));
}
};
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1,
final TimeScales timeScales) {
// Earth Rotation Angle
final StellarAngleCapitaine era =
new StellarAngleCapitaine(ut1, timeScales.getTAI());
// Polynomial part of the apparent sidereal time series
// which is the opposite of Equation of Origins (EO)
final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
final PoissonSeriesParser parser =
new PoissonSeriesParser(17).
withFirstDelaunay(4).
withFirstPlanetary(9).
withSinCos(0, 2, microAS, 3, microAS).
withPolynomialPart('t', Unit.ARC_SECONDS);
final PolynomialNutation minusEO = load(GST_SERIES, in -> parser.parse(in, GST_SERIES).getPolynomial());
// create a function evaluating the series
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
return era.getRate() + minusEO.derivative(evaluateTC(date, timeScales));
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
return minusEO.derivative(evaluateTC(date, timeScales)).add(era.getRate());
}
};
}
/** {@inheritDoc} */
@Override
public TimeScalarFunction getGASTFunction(final TimeScale ut1,
final EOPHistory eopHistory,
final TimeScales timeScales) {
// set up nutation arguments
final FundamentalNutationArguments arguments =
getNutationArguments(timeScales);
// mean obliquity function
final TimeScalarFunction epsilon = getMeanObliquityFunction(timeScales);
// set up Poisson series
final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
final PoissonSeriesParser baseParser =
new PoissonSeriesParser(17).
withFirstDelaunay(4).
withFirstPlanetary(9).
withSinCos(0, 2, microAS, 3, microAS);
final PoissonSeriesParser gstParser = baseParser.withPolynomialPart('t', Unit.ARC_SECONDS);
final PoissonSeries psiSeries = load(PSI_SERIES, in -> baseParser.parse(in, PSI_SERIES));
final PoissonSeries gstSeries = load(GST_SERIES, in -> gstParser.parse(in, GST_SERIES));
final PoissonSeries.CompiledSeries psiGstSeries = PoissonSeries.compile(psiSeries, gstSeries);
// ERA function
final TimeScalarFunction era =
getEarthOrientationAngleFunction(ut1, timeScales.getTAI());
return new TimeScalarFunction() {
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
// evaluate equation of origins
final BodiesElements elements = arguments.evaluateAll(date);
final double[] angles = psiGstSeries.value(elements);
final double ddPsi = (eopHistory == null) ? 0 : eopHistory.getEquinoxNutationCorrection(date)[0];
final double deltaPsi = angles[0] + ddPsi;
final double epsilonA = epsilon.value(date);
// subtract equation of origin from EA
// (hence add the series above which have the sign included)
return era.value(date) + deltaPsi * FastMath.cos(epsilonA) + angles[1];
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
// evaluate equation of origins
final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
final T[] angles = psiGstSeries.value(elements);
final T ddPsi = (eopHistory == null) ? date.getField().getZero() : eopHistory.getEquinoxNutationCorrection(date)[0];
final T deltaPsi = angles[0].add(ddPsi);
final T epsilonA = epsilon.value(date);
// subtract equation of origin from EA
// (hence add the series above which have the sign included)
return era.value(date).add(deltaPsi.multiply(epsilonA.cos())).add(angles[1]);
}
};
}
/** {@inheritDoc} */
@Override
public TimeVectorFunction getEOPTidalCorrection(final TimeScales timeScales) {
// set up nutation arguments
// BEWARE! Using TT as the time scale here and not UT1 is intentional!
// as this correction is used to compute UT1 itself, it is not surprising we cannot use UT1 yet,
// however, using the close UTC as would seem logical make the comparison with interp.f from IERS fail
// looking in the interp.f code, the same TT scale is used for both Delaunay and gamma argument
final FundamentalNutationArguments arguments =
getNutationArguments(timeScales.getTT(), timeScales);
// set up Poisson series
final double microAS = Constants.ARC_SECONDS_TO_RADIANS * 1.0e-6;
final PoissonSeriesParser xyParser = new PoissonSeriesParser(13).
withOptionalColumn(1).
withGamma(2).
withFirstDelaunay(3);
final double microS = 1.0e-6;
final PoissonSeriesParser ut1Parser = new PoissonSeriesParser(11).
withOptionalColumn(1).
withGamma(2).
withFirstDelaunay(3);
final PoissonSeries xSeries =
load(TIDAL_CORRECTION_XP_YP_SERIES, xpIn -> xyParser.
withSinCos(0, 10, microAS, 11, microAS).
parse(xpIn, TIDAL_CORRECTION_XP_YP_SERIES));
final PoissonSeries ySeries =
load(TIDAL_CORRECTION_XP_YP_SERIES, ypIn -> xyParser.
withSinCos(0, 12, microAS, 13, microAS).
parse(ypIn, TIDAL_CORRECTION_XP_YP_SERIES));
final PoissonSeries ut1Series =
load(TIDAL_CORRECTION_UT1_SERIES, ut1In -> ut1Parser.
withSinCos(0, 10, microS, 11, microS).
parse(ut1In, TIDAL_CORRECTION_UT1_SERIES));
return new EOPTidalCorrection(arguments, xSeries, ySeries, ut1Series);
}
/** {@inheritDoc} */
@Override
public double[] getNominalTidalDisplacement() {
return new double[] {
// h⁽⁰⁾, h⁽²⁾, h₃, hI diurnal, hI semi-diurnal,
0.6078, -0.0006, 0.292, -0.0025, -0.0022,
// l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾, l₃, lI diurnal, lI semi-diurnal
0.0847, 0.0012, 0.0024, 0.0002, 0.015, -0.0007, -0.0007,
// H₀
-0.31460
};
}
/** {@inheritDoc} */
@Override
public CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal() {
return getTidalDisplacementFrequencyCorrectionDiurnal(TIDAL_DISPLACEMENT_CORRECTION_DIURNAL,
18, 15, 16, 17, 18);
}
/** {@inheritDoc} */
@Override
public CompiledSeries getTidalDisplacementFrequencyCorrectionZonal() {
return getTidalDisplacementFrequencyCorrectionZonal(TIDAL_DISPLACEMENT_CORRECTION_ZONAL,
18, 15, 16, 17, 18);
}
};
/** Pattern for delimiting regular expressions. */
private static final Pattern SEPARATOR = Pattern.compile("\\p{Space}+");
/** IERS conventions resources base directory. */
private static final String IERS_BASE = "/assets/org/orekit/IERS-conventions/";
/** Get the reference epoch for fundamental nutation arguments.
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @return reference epoch for fundamental nutation arguments
* @since 6.1
* @see #getNutationReferenceEpoch(TimeScales)
*/
@DefaultDataContext
public AbsoluteDate getNutationReferenceEpoch() {
return getNutationReferenceEpoch(DataContext.getDefault().getTimeScales());
}
/**
* Get the reference epoch for fundamental nutation arguments.
*
* @param timeScales to use for the reference epoch.
* @return reference epoch for fundamental nutation arguments
* @since 10.1
*/
public AbsoluteDate getNutationReferenceEpoch(final TimeScales timeScales) {
// IERS 1996, IERS 2003 and IERS 2010 use the same J2000.0 reference date
return timeScales.getJ2000Epoch();
}
/** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @param date current date
* @return date offset in Julian centuries
* @since 6.1
* @see #evaluateTC(AbsoluteDate, TimeScales)
*/
@DefaultDataContext
public double evaluateTC(final AbsoluteDate date) {
return evaluateTC(date, DataContext.getDefault().getTimeScales());
}
/**
* Evaluate the date offset between the current date and the {@link
* #getNutationReferenceEpoch() reference date}.
*
* @param date current date
* @param timeScales used in the evaluation.
* @return date offset in Julian centuries
* @since 10.1
*/
public double evaluateTC(final AbsoluteDate date, final TimeScales timeScales) {
return date.durationFrom(getNutationReferenceEpoch(timeScales)) /
Constants.JULIAN_CENTURY;
}
/** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @param date current date
* @param <T> type of the field elements
* @return date offset in Julian centuries
* @since 9.0
* @see #evaluateTC(FieldAbsoluteDate, TimeScales)
*/
@DefaultDataContext
public <T extends CalculusFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date) {
return evaluateTC(date, DataContext.getDefault().getTimeScales());
}
/** Evaluate the date offset between the current date and the {@link #getNutationReferenceEpoch() reference date}.
* @param <T> type of the field elements
* @param date current date
* @param timeScales used in the evaluation.
* @return date offset in Julian centuries
* @since 10.1
*/
public <T extends CalculusFieldElement<T>> T evaluateTC(final FieldAbsoluteDate<T> date,
final TimeScales timeScales) {
return date.durationFrom(getNutationReferenceEpoch(timeScales))
.divide(Constants.JULIAN_CENTURY);
}
/**
* Get the fundamental nutation arguments. Does not compute GMST based values: gamma,
* gammaDot.
*
* @param timeScales other time scales used in the computation including TAI and TT.
* @return fundamental nutation arguments
* @see #getNutationArguments(TimeScale, TimeScales)
* @since 10.1
*/
protected FundamentalNutationArguments getNutationArguments(
final TimeScales timeScales) {
return getNutationArguments(null, timeScales);
}
/** Get the fundamental nutation arguments.
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @param timeScale time scale for computing Greenwich Mean Sidereal Time
* (typically {@link TimeScales#getUT1(IERSConventions, boolean) UT1})
* @return fundamental nutation arguments
* @since 6.1
* @see #getNutationArguments(TimeScale, TimeScales)
* @see #getNutationArguments(TimeScales)
*/
@DefaultDataContext
public FundamentalNutationArguments getNutationArguments(final TimeScale timeScale) {
return getNutationArguments(timeScale, DataContext.getDefault().getTimeScales());
}
/**
* Get the fundamental nutation arguments.
*
* @param timeScale time scale for computing Greenwich Mean Sidereal Time (typically
* {@link TimeScales#getUT1(IERSConventions, boolean) UT1})
* @param timeScales other time scales used in the computation including TAI and TT.
* @return fundamental nutation arguments
* @since 10.1
*/
public abstract FundamentalNutationArguments getNutationArguments(
TimeScale timeScale,
TimeScales timeScales);
/** Get the function computing mean obliquity of the ecliptic.
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @return function computing mean obliquity of the ecliptic
* @since 6.1
* @see #getMeanObliquityFunction(TimeScales)
*/
@DefaultDataContext
public TimeScalarFunction getMeanObliquityFunction() {
return getMeanObliquityFunction(DataContext.getDefault().getTimeScales());
}
/**
* Get the function computing mean obliquity of the ecliptic.
*
* @param timeScales used in computing the function.
* @return function computing mean obliquity of the ecliptic
* @since 10.1
*/
public abstract TimeScalarFunction getMeanObliquityFunction(TimeScales timeScales);
/** Get the function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components.
* <p>
* The returned function computes the two X, Y components of CIP and the S+XY/2 component of the non-rotating CIO.
* </p>
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @return function computing the Celestial Intermediate Pole and Celestial Intermediate Origin components
* @since 6.1
* @see #getXYSpXY2Function(TimeScales)
*/
@DefaultDataContext
public TimeVectorFunction getXYSpXY2Function() {
return getXYSpXY2Function(DataContext.getDefault().getTimeScales());
}
/**
* Get the function computing the Celestial Intermediate Pole and Celestial
* Intermediate Origin components.
* <p>
* The returned function computes the two X, Y components of CIP and the S+XY/2
* component of the non-rotating CIO.
* </p>
*
* @param timeScales used to define the function.
* @return function computing the Celestial Intermediate Pole and Celestial
* Intermediate Origin components
* @since 10.1
*/
public abstract TimeVectorFunction getXYSpXY2Function(TimeScales timeScales);
/** Get the function computing the raw Earth Orientation Angle.
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* <p>
* The raw angle does not contain any correction. If for example dTU1 correction
* due to tidal effect is desired, it must be added afterward by the caller.
* The returned value contain the angle as the value and the angular rate as
* the first derivative.
* </p>
* @param ut1 UT1 time scale
* @return function computing the rawEarth Orientation Angle, in the non-rotating origin paradigm
* @since 6.1
* @see #getEarthOrientationAngleFunction(TimeScale, TimeScale)
*/
@DefaultDataContext
public TimeScalarFunction getEarthOrientationAngleFunction(final TimeScale ut1) {
return getEarthOrientationAngleFunction(ut1,
DataContext.getDefault().getTimeScales().getTAI());
}
/** Get the function computing the raw Earth Orientation Angle.
* <p>
* The raw angle does not contain any correction. If for example dTU1 correction
* due to tidal effect is desired, it must be added afterward by the caller.
* The returned value contain the angle as the value and the angular rate as
* the first derivative.
* </p>
* @param ut1 UT1 time scale
* @param tai TAI time scale
* @return function computing the rawEarth Orientation Angle, in the non-rotating origin paradigm
* @since 10.1
*/
public TimeScalarFunction getEarthOrientationAngleFunction(final TimeScale ut1,
final TimeScale tai) {
return new StellarAngleCapitaine(ut1, tai);
}
/** Get the function computing the precession angles.
* <p>
* The function returned computes the three precession angles
* ψ<sub>A</sub> (around Z axis), ω<sub>A</sub> (around X axis)
* and χ<sub>A</sub> (around Z axis). The constant angle ε₀
* for the fourth rotation (around X axis) can be retrieved by evaluating the
* function returned by {@link #getMeanObliquityFunction()} at {@link
* #getNutationReferenceEpoch() nutation reference epoch}.
* </p>
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @return function computing the precession angle
* @since 6.1
* @see #getPrecessionFunction(TimeScales)
*/
@DefaultDataContext
public TimeVectorFunction getPrecessionFunction()
{
return getPrecessionFunction(DataContext.getDefault().getTimeScales());
}
/** Get the function computing the precession angles.
* <p>
* The function returned computes the three precession angles
* ψ<sub>A</sub> (around Z axis), ω<sub>A</sub> (around X axis)
* and χ<sub>A</sub> (around Z axis). The constant angle ε₀
* for the fourth rotation (around X axis) can be retrieved by evaluating the
* function returned by {@link #getMeanObliquityFunction()} at {@link
* #getNutationReferenceEpoch() nutation reference epoch}.
* </p>
* @return function computing the precession angle
* @since 10.1
* @param timeScales used to define the function.
*/
public abstract TimeVectorFunction getPrecessionFunction(TimeScales timeScales);
/** Get the function computing the nutation angles.
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* <p>
* The function returned computes the two classical angles ΔΨ and Δε,
* and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
* resolution C7 (the correction is forced to 0 before this date)
* </p>
* @return function computing the nutation in longitude ΔΨ and Δε
* and the correction of equation of equinoxes
* @since 6.1
*/
@DefaultDataContext
public TimeVectorFunction getNutationFunction() {
return getNutationFunction(DataContext.getDefault().getTimeScales());
}
/** Get the function computing the nutation angles.
* <p>
* The function returned computes the two classical angles ΔΨ and Δε,
* and the correction to the equation of equinoxes introduced since 1997-02-27 by IAU 1994
* resolution C7 (the correction is forced to 0 before this date)
* </p>
* @return function computing the nutation in longitude ΔΨ and Δε
* and the correction of equation of equinoxes
* @param timeScales used in the computation including TAI and TT.
* @since 10.1
*/
public abstract TimeVectorFunction getNutationFunction(TimeScales timeScales);
/** Get the function computing Greenwich mean sidereal time, in radians.
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @param ut1 UT1 time scale
* @return function computing Greenwich mean sidereal time
* @since 6.1
* @see #getGMSTFunction(TimeScale, TimeScales)
*/
@DefaultDataContext
public TimeScalarFunction getGMSTFunction(final TimeScale ut1) {
return getGMSTFunction(ut1, DataContext.getDefault().getTimeScales());
}
/**
* Get the function computing Greenwich mean sidereal time, in radians.
*
* @param ut1 UT1 time scale
* @param timeScales other time scales used in the computation including TAI and TT.
* @return function computing Greenwich mean sidereal time
* @since 10.1
*/
public abstract TimeScalarFunction getGMSTFunction(TimeScale ut1,
TimeScales timeScales);
/** Get the function computing Greenwich mean sidereal time rate, in radians per second.
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @param ut1 UT1 time scale
* @return function computing Greenwich mean sidereal time rate
* @since 9.0
* @see #getGMSTRateFunction(TimeScale, TimeScales)
*/
@DefaultDataContext
public TimeScalarFunction getGMSTRateFunction(final TimeScale ut1) {
return getGMSTRateFunction(ut1,
DataContext.getDefault().getTimeScales());
}
/**
* Get the function computing Greenwich mean sidereal time rate, in radians per
* second.
*
* @param ut1 UT1 time scale
* @param timeScales other time scales used in the computation including TAI and TT.
* @return function computing Greenwich mean sidereal time rate
* @since 10.1
*/
public abstract TimeScalarFunction getGMSTRateFunction(TimeScale ut1,
TimeScales timeScales);
/**
* Get the function computing Greenwich apparent sidereal time, in radians.
*
* <p>This method uses the {@link DataContext#getDefault() default data context} if
* {@code eopHistory == null}.
*
* @param ut1 UT1 time scale
* @param eopHistory EOP history. If {@code null} then no nutation correction is
* applied for EOP.
* @return function computing Greenwich apparent sidereal time
* @since 6.1
* @see #getGASTFunction(TimeScale, EOPHistory, TimeScales)
*/
@DefaultDataContext
public TimeScalarFunction getGASTFunction(final TimeScale ut1,
final EOPHistory eopHistory) {
final TimeScales timeScales = eopHistory != null ?
eopHistory.getTimeScales() :
DataContext.getDefault().getTimeScales();
return getGASTFunction(ut1, eopHistory, timeScales);
}
/**
* Get the function computing Greenwich apparent sidereal time, in radians.
*
* @param ut1 UT1 time scale
* @param eopHistory EOP history. If {@code null} then no nutation correction is
* applied for EOP.
* @param timeScales TAI time scale.
* @return function computing Greenwich apparent sidereal time
* @since 10.1
*/
public abstract TimeScalarFunction getGASTFunction(TimeScale ut1,
EOPHistory eopHistory,
TimeScales timeScales);
/** Get the function computing tidal corrections for Earth Orientation Parameters.
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @return function computing tidal corrections for Earth Orientation Parameters,
* for xp, yp, ut1 and lod respectively
* @since 6.1
* @see #getEOPTidalCorrection(TimeScales)
*/
@DefaultDataContext
public TimeVectorFunction getEOPTidalCorrection() {
return getEOPTidalCorrection(DataContext.getDefault().getTimeScales());
}
/**
* Get the function computing tidal corrections for Earth Orientation Parameters.
*
* @param timeScales used in the computation. The TT and TAI scales are used.
* @return function computing tidal corrections for Earth Orientation Parameters, for
* xp, yp, ut1 and lod respectively
* @since 10.1
*/
public abstract TimeVectorFunction getEOPTidalCorrection(TimeScales timeScales);
/** Get the Love numbers.
* @return Love numbers
* @since 6.1
*/
public abstract LoveNumbers getLoveNumbers();
/** Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @param ut1 UT1 time scale
* @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
* @since 6.1
* @see #getTideFrequencyDependenceFunction(TimeScale, TimeScales)
*/
@DefaultDataContext
public TimeVectorFunction getTideFrequencyDependenceFunction(final TimeScale ut1) {
return getTideFrequencyDependenceFunction(ut1,
DataContext.getDefault().getTimeScales());
}
/**
* Get the function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂,
* ΔS₂₂).
*
* @param ut1 UT1 time scale
* @param timeScales other time scales used in the computation including TAI and TT.
* @return frequency dependence model for tides computation (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂,
* ΔS₂₂).
* @since 10.1
*/
public abstract TimeVectorFunction getTideFrequencyDependenceFunction(
TimeScale ut1,
TimeScales timeScales);
/** Get the permanent tide to be <em>removed</em> from ΔC₂₀ when zero-tide potentials are used.
* @return permanent tide to remove
*/
public abstract double getPermanentTide();
/** Get the function computing solid pole tide (ΔC₂₁, ΔS₂₁).
* @param eopHistory EOP history
* @return model for solid pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
* @since 6.1
*/
public abstract TimeVectorFunction getSolidPoleTide(EOPHistory eopHistory);
/** Get the function computing ocean pole tide (ΔC₂₁, ΔS₂₁).
* @param eopHistory EOP history
* @return model for ocean pole tide (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂).
* @since 6.1
*/
public abstract TimeVectorFunction getOceanPoleTide(EOPHistory eopHistory);
/** Get the nominal values of the displacement numbers.
* @return an array containing h⁽⁰⁾, h⁽²⁾, h₃, hI diurnal, hI semi-diurnal,
* l⁽⁰⁾, l⁽¹⁾ diurnal, l⁽¹⁾ semi-diurnal, l⁽²⁾, l₃, lI diurnal, lI semi-diurnal,
* H₀ permanent deformation amplitude
* @since 9.1
*/
public abstract double[] getNominalTidalDisplacement();
/** Get the correction function for tidal displacement for diurnal tides.
* <ul>
* <li>f[0]: radial correction, longitude cosine part</li>
* <li>f[1]: radial correction, longitude sine part</li>
* <li>f[2]: North correction, longitude cosine part</li>
* <li>f[3]: North correction, longitude sine part</li>
* <li>f[4]: East correction, longitude cosine part</li>
* <li>f[5]: East correction, longitude sine part</li>
* </ul>
* @return correction function for tidal displacement
* @since 9.1
*/
public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal();
/** Get the correction function for tidal displacement for diurnal tides.
* <ul>
* <li>f[0]: radial correction, longitude cosine part</li>
* <li>f[1]: radial correction, longitude sine part</li>
* <li>f[2]: North correction, longitude cosine part</li>
* <li>f[3]: North correction, longitude sine part</li>
* <li>f[4]: East correction, longitude cosine part</li>
* <li>f[5]: East correction, longitude sine part</li>
* </ul>
* @param tableName name for the diurnal tides table
* @param cols total number of columns of the diurnal tides table
* @param rIp column holding ∆Rf(ip) in the diurnal tides table, counting from 1
* @param rOp column holding ∆Rf(op) in the diurnal tides table, counting from 1
* @param tIp column holding ∆Tf(ip) in the diurnal tides table, counting from 1
* @param tOp column holding ∆Tf(op) in the diurnal tides table, counting from 1
* @return correction function for tidal displacement for diurnal tides
* @since 9.1
*/
protected static CompiledSeries getTidalDisplacementFrequencyCorrectionDiurnal(final String tableName, final int cols,
final int rIp, final int rOp,
final int tIp, final int tOp) {
// radial component, missing the sin 2φ factor; this corresponds to:
// - equation 15a in IERS conventions 1996, chapter 7
// - equation 16a in IERS conventions 2003, chapter 7
// - equation 7.12a in IERS conventions 2010, chapter 7
final PoissonSeries drCos = load(tableName, in -> new PoissonSeriesParser(cols).
withOptionalColumn(1).
withDoodson(4, 3).
withFirstDelaunay(10).
withSinCos(0, rIp, +1.0e-3, rOp, +1.0e-3).
parse(in, tableName));
final PoissonSeries drSin = load(tableName, in -> new PoissonSeriesParser(cols).
withOptionalColumn(1).
withDoodson(4, 3).
withFirstDelaunay(10).
withSinCos(0, rOp, -1.0e-3, rIp, +1.0e-3).
parse(in, tableName));
// North component, missing the cos 2φ factor; this corresponds to:
// - equation 15b in IERS conventions 1996, chapter 7
// - equation 16b in IERS conventions 2003, chapter 7
// - equation 7.12b in IERS conventions 2010, chapter 7
final PoissonSeries dnCos = load(tableName, in -> new PoissonSeriesParser(cols).
withOptionalColumn(1).
withDoodson(4, 3).
withFirstDelaunay(10).
withSinCos(0, tIp, +1.0e-3, tOp, +1.0e-3).
parse(in, tableName));
final PoissonSeries dnSin = load(tableName, in -> new PoissonSeriesParser(cols).
withOptionalColumn(1).
withDoodson(4, 3).
withFirstDelaunay(10).
withSinCos(0, tOp, -1.0e-3, tIp, +1.0e-3).
parse(in, tableName));
// East component, missing the sin φ factor; this corresponds to:
// - equation 15b in IERS conventions 1996, chapter 7
// - equation 16b in IERS conventions 2003, chapter 7
// - equation 7.12b in IERS conventions 2010, chapter 7
final PoissonSeries deCos = load(tableName, in -> new PoissonSeriesParser(cols).
withOptionalColumn(1).
withDoodson(4, 3).
withFirstDelaunay(10).
withSinCos(0, tOp, -1.0e-3, tIp, +1.0e-3).
parse(in, tableName));
final PoissonSeries deSin = load(tableName, in -> new PoissonSeriesParser(cols).
withOptionalColumn(1).
withDoodson(4, 3).
withFirstDelaunay(10).
withSinCos(0, tIp, -1.0e-3, tOp, -1.0e-3).
parse(in, tableName));
return PoissonSeries.compile(drCos, drSin, dnCos, dnSin, deCos, deSin);
}
/** Get the correction function for tidal displacement for zonal tides.
* <ul>
* <li>f[0]: radial correction</li>
* <li>f[1]: North correction</li>
* </ul>
* @return correction function for tidal displacement
* @since 9.1
*/
public abstract CompiledSeries getTidalDisplacementFrequencyCorrectionZonal();
/** Get the correction function for tidal displacement for zonal tides.
* <ul>
* <li>f[0]: radial correction</li>
* <li>f[1]: North correction</li>
* </ul>
* @param tableName name for the zonal tides table
* @param cols total number of columns of the table
* @param rIp column holding ∆Rf(ip) in the table, counting from 1
* @param rOp column holding ∆Rf(op) in the table, counting from 1
* @param tIp column holding ∆Tf(ip) in the table, counting from 1
* @param tOp column holding ∆Tf(op) in the table, counting from 1
* @return correction function for tidal displacement for zonal tides
* @since 9.1
*/
protected static CompiledSeries getTidalDisplacementFrequencyCorrectionZonal(final String tableName, final int cols,
final int rIp, final int rOp,
final int tIp, final int tOp) {
// radial component, missing the 3⁄2 sin² φ - 1⁄2 factor; this corresponds to:
// - equation 16a in IERS conventions 1996, chapter 7
// - equation 17a in IERS conventions 2003, chapter 7
// - equation 7.13a in IERS conventions 2010, chapter 7
final PoissonSeries dr = load(tableName, in -> new PoissonSeriesParser(cols).
withOptionalColumn(1).
withDoodson(4, 3).
withFirstDelaunay(10).
withSinCos(0, rOp, +1.0e-3, rIp, +1.0e-3).
parse(in, tableName));
// North component, missing the sin 2φ factor; this corresponds to:
// - equation 16b in IERS conventions 1996, chapter 7
// - equation 17b in IERS conventions 2003, chapter 7
// - equation 7.13b in IERS conventions 2010, chapter 7
final PoissonSeries dn = load(tableName, in -> new PoissonSeriesParser(cols).
withOptionalColumn(1).
withDoodson(4, 3).
withFirstDelaunay(10).
withSinCos(0, tOp, +1.0e-3, tIp, +1.0e-3).
parse(in, tableName));
return PoissonSeries.compile(dr, dn);
}
/** Interface for functions converting nutation corrections between
* δΔψ/δΔε to δX/δY.
* <ul>
* <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
* <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
* </ul>
* @since 6.1
*/
public interface NutationCorrectionConverter {
/** Convert nutation corrections.
* @param date current date
* @param ddPsi δΔψ part of the nutation correction
* @param ddEpsilon δΔε part of the nutation correction
* @return array containing δX and δY
*/
double[] toNonRotating(AbsoluteDate date, double ddPsi, double ddEpsilon);
/** Convert nutation corrections.
* @param date current date
* @param dX δX part of the nutation correction
* @param dY δY part of the nutation correction
* @return array containing δΔψ and δΔε
*/
double[] toEquinox(AbsoluteDate date, double dX, double dY);
}
/** Create a function converting nutation corrections between
* δX/δY and δΔψ/δΔε.
* <ul>
* <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
* <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
* </ul>
*
* <p>This method uses the {@link DataContext#getDefault() default data context}.
*
* @return a new converter
* @since 6.1
* @see #getNutationCorrectionConverter(TimeScales)
*/
@DefaultDataContext
public NutationCorrectionConverter getNutationCorrectionConverter() {
return getNutationCorrectionConverter(DataContext.getDefault().getTimeScales());
}
/** Create a function converting nutation corrections between
* δX/δY and δΔψ/δΔε.
* <ul>
* <li>δX/δY nutation corrections are used with the Non-Rotating Origin paradigm.</li>
* <li>δΔψ/δΔε nutation corrections are used with the equinox-based paradigm.</li>
* </ul>
* @return a new converter
* @since 10.1
* @param timeScales used to define the conversion.
*/
public NutationCorrectionConverter getNutationCorrectionConverter(
final TimeScales timeScales) {
// get models parameters
final TimeVectorFunction precessionFunction = getPrecessionFunction(timeScales);
final TimeScalarFunction epsilonAFunction = getMeanObliquityFunction(timeScales);
final AbsoluteDate date0 = getNutationReferenceEpoch(timeScales);
final double cosE0 = FastMath.cos(epsilonAFunction.value(date0));
return new NutationCorrectionConverter() {
/** {@inheritDoc} */
@Override
public double[] toNonRotating(final AbsoluteDate date,
final double ddPsi, final double ddEpsilon) {
// compute precession angles psiA, omegaA and chiA
final double[] angles = precessionFunction.value(date);
// conversion coefficients
final double sinEA = FastMath.sin(epsilonAFunction.value(date));
final double c = angles[0] * cosE0 - angles[2];
// convert nutation corrections (equation 23/IERS-2003 or 5.25/IERS-2010)
return new double[] {
sinEA * ddPsi + c * ddEpsilon,
ddEpsilon - c * sinEA * ddPsi
};
}
/** {@inheritDoc} */
@Override
public double[] toEquinox(final AbsoluteDate date,
final double dX, final double dY) {
// compute precession angles psiA, omegaA and chiA
final double[] angles = precessionFunction.value(date);
// conversion coefficients
final double sinEA = FastMath.sin(epsilonAFunction.value(date));
final double c = angles[0] * cosE0 - angles[2];
final double opC2 = 1 + c * c;
// convert nutation corrections (inverse of equation 23/IERS-2003 or 5.25/IERS-2010)
return new double[] {
(dX - c * dY) / (sinEA * opC2),
(dY + c * dX) / opC2
};
}
};
}
/** Load the Love numbers.
* @param nameLove name of the Love number resource
* @return Love numbers
*/
protected LoveNumbers loadLoveNumbers(final String nameLove) {
try {
// allocate the three triangular arrays for real, imaginary and time-dependent numbers
final double[][] real = new double[4][];
final double[][] imaginary = new double[4][];
final double[][] plus = new double[4][];
for (int i = 0; i < real.length; ++i) {
real[i] = new double[i + 1];
imaginary[i] = new double[i + 1];
plus[i] = new double[i + 1];
}
try (InputStream stream = IERSConventions.class.getResourceAsStream(nameLove)) {
if (stream == null) {
// this should never happen with files embedded within Orekit
throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, nameLove);
}
int lineNumber = 1;
String line = null;
// setup the reader
try (BufferedReader reader = new BufferedReader(new InputStreamReader(stream, StandardCharsets.UTF_8))) {
line = reader.readLine();
// look for the Love numbers
while (line != null) {
line = line.trim();
if (!(line.isEmpty() || line.startsWith("#"))) {
final String[] fields = SEPARATOR.split(line);
if (fields.length != 5) {
// this should never happen with files embedded within Orekit
throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
lineNumber, nameLove, line);
}
final int n = Integer.parseInt(fields[0]);
final int m = Integer.parseInt(fields[1]);
if (n < 2 || n > 3 || m < 0 || m > n) {
// this should never happen with files embedded within Orekit
throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
lineNumber, nameLove, line);
}
real[n][m] = Double.parseDouble(fields[2]);
imaginary[n][m] = Double.parseDouble(fields[3]);
plus[n][m] = Double.parseDouble(fields[4]);
}
// next line
lineNumber++;
line = reader.readLine();
}
} catch (NumberFormatException nfe) {
// this should never happen with files embedded within Orekit
throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
lineNumber, nameLove, line);
}
}
return new LoveNumbers(real, imaginary, plus);
} catch (IOException ioe) {
// this should never happen with files embedded within Orekit
throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, nameLove);
}
}
/** Load resources.
* @param name name of the resource
* @param loader loaader for the resource
* @param <T> type of the processed data
* @return processed data
*/
private static <T> T load(final String name, final Function<InputStream, T> loader) {
try (InputStream is = IERSConventions.class.getResourceAsStream(name)) {
return loader.apply(is);
} catch (IOException ioe) {
// this should never happen with internal streams
throw new OrekitException(OrekitMessages.INTERNAL_ERROR, ioe);
}
}
/** Correction to equation of equinoxes.
* <p>IAU 1994 resolution C7 added two terms to the equation of equinoxes
* taking effect since 1997-02-27 for continuity
* </p>
*/
private static class IAU1994ResolutionC7 {
/** First Moon correction term for the Equation of the Equinoxes. */
private static final double EQE1 = 0.00264 * Constants.ARC_SECONDS_TO_RADIANS;
/** Second Moon correction term for the Equation of the Equinoxes. */
private static final double EQE2 = 0.000063 * Constants.ARC_SECONDS_TO_RADIANS;
/** Evaluate the correction.
* @param arguments Delaunay for nutation
* @param tai TAI time scale.
* @return correction value (0 before 1997-02-27)
*/
public static double value(final DelaunayArguments arguments,
final TimeScale tai) {
/* Start date for applying Moon corrections to the equation of the equinoxes.
* This date corresponds to 1997-02-27T00:00:00 UTC, hence the 30s offset from TAI.
*/
final AbsoluteDate modelStart = new AbsoluteDate(1997, 2, 27, 0, 0, 30, tai);
if (arguments.getDate().compareTo(modelStart) >= 0) {
// IAU 1994 resolution C7 added two terms to the equation of equinoxes
// taking effect since 1997-02-27 for continuity
// Mean longitude of the ascending node of the Moon
final double om = arguments.getOmega();
// add the two correction terms
return EQE1 * FastMath.sin(om) + EQE2 * FastMath.sin(om + om);
} else {
return 0.0;
}
}
/** Evaluate the correction.
* @param arguments Delaunay for nutation
* @param tai TAI time scale.
* @param <T> type of the field elements
* @return correction value (0 before 1997-02-27)
*/
public static <T extends CalculusFieldElement<T>> T value(
final FieldDelaunayArguments<T> arguments,
final TimeScale tai) {
/* Start date for applying Moon corrections to the equation of the equinoxes.
* This date corresponds to 1997-02-27T00:00:00 UTC, hence the 30s offset from TAI.
*/
final AbsoluteDate modelStart = new AbsoluteDate(1997, 2, 27, 0, 0, 30, tai);
if (arguments.getDate().toAbsoluteDate().compareTo(modelStart) >= 0) {
// IAU 1994 resolution C7 added two terms to the equation of equinoxes
// taking effect since 1997-02-27 for continuity
// Mean longitude of the ascending node of the Moon
final T om = arguments.getOmega();
// add the two correction terms
return om.sin().multiply(EQE1).add(om.add(om).sin().multiply(EQE2));
} else {
return arguments.getDate().getField().getZero();
}
}
};
/** Stellar angle model.
* <p>
* The stellar angle computed here has been defined in the paper "A non-rotating origin on the
* instantaneous equator: Definition, properties and use", N. Capitaine, Guinot B., and Souchay J.,
* Celestial Mechanics, Volume 39, Issue 3, pp 283-307. It has been proposed as a conventional
* conventional relationship between UT1 and Earth rotation in the paper "Definition of the Celestial
* Ephemeris origin and of UT1 in the International Celestial Reference Frame", Capitaine, N.,
* Guinot, B., and McCarthy, D. D., 2000, Astronomy and Astrophysics, 355(1), pp. 398–405.
* </p>
* <p>
* It is presented simply as stellar angle in IERS conventions 1996 but as since been adopted as
* the conventional relationship defining UT1 from ICRF and is called Earth Rotation Angle in
* IERS conventions 2003 and 2010.
* </p>
*/
private static class StellarAngleCapitaine implements TimeScalarFunction {
/** Constant term of Capitaine's Earth Rotation Angle model. */
private static final double ERA_0 = MathUtils.TWO_PI * 0.7790572732640;
/** Rate term of Capitaine's Earth Rotation Angle model.
* (radians per day, main part) */
private static final double ERA_1A = MathUtils.TWO_PI / Constants.JULIAN_DAY;
/** Rate term of Capitaine's Earth Rotation Angle model.
* (radians per day, fractional part) */
private static final double ERA_1B = ERA_1A * 0.00273781191135448;
/** UT1 time scale. */
private final TimeScale ut1;
/** Reference date of Capitaine's Earth Rotation Angle model. */
private final AbsoluteDate referenceDate;
/** Simple constructor.
* @param ut1 UT1 time scale
* @param tai TAI time scale
*/
StellarAngleCapitaine(final TimeScale ut1, final TimeScale tai) {
this.ut1 = ut1;
referenceDate = new AbsoluteDate(
DateComponents.J2000_EPOCH,
TimeComponents.H12,
tai);
}
/** Get the rotation rate.
* @return rotation rate
*/
public double getRate() {
return ERA_1A + ERA_1B;
}
/** {@inheritDoc} */
@Override
public double value(final AbsoluteDate date) {
// split the date offset as a full number of days plus a smaller part
final int secondsInDay = 86400;
final double dt = date.durationFrom(referenceDate);
final long days = ((long) dt) / secondsInDay;
final double dtA = secondsInDay * days;
final double dtB = (dt - dtA) + ut1.offsetFromTAI(date);
return ERA_0 + ERA_1A * dtB + ERA_1B * (dtA + dtB);
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T value(final FieldAbsoluteDate<T> date) {
// split the date offset as a full number of days plus a smaller part
final int secondsInDay = 86400;
final T dt = date.durationFrom(referenceDate);
final long days = ((long) dt.getReal()) / secondsInDay;
final double dtA = secondsInDay * days;
final T dtB = dt.subtract(dtA).add(ut1.offsetFromTAI(date.toAbsoluteDate()));
return dtB.add(dtA).multiply(ERA_1B).add(dtB.multiply(ERA_1A)).add(ERA_0);
}
}
/** Mean pole. */
private static class MeanPole implements TimeStamped, Serializable {
/** Serializable UID. */
private static final long serialVersionUID = 20131028l;
/** Date. */
private final AbsoluteDate date;
/** X coordinate. */
private double x;
/** Y coordinate. */
private double y;
/** Simple constructor.
* @param date date
* @param x x coordinate
* @param y y coordinate
*/
MeanPole(final AbsoluteDate date, final double x, final double y) {
this.date = date;
this.x = x;
this.y = y;
}
/** {@inheritDoc} */
@Override
public AbsoluteDate getDate() {
return date;
}
/** Get x coordinate.
* @return x coordinate
*/
public double getX() {
return x;
}
/** Get y coordinate.
* @return y coordinate
*/
public double getY() {
return y;
}
}
/** Local class for precession function. */
private class PrecessionFunction implements TimeVectorFunction {
/** Polynomial nutation for psiA. */
private final PolynomialNutation psiA;
/** Polynomial nutation for omegaA. */
private final PolynomialNutation omegaA;
/** Polynomial nutation for chiA. */
private final PolynomialNutation chiA;
/** Time scales to use. */
private final TimeScales timeScales;
/** Simple constructor.
* @param psiA polynomial nutation for psiA
* @param omegaA polynomial nutation for omegaA
* @param chiA polynomial nutation for chiA
* @param timeScales used in the computation.
*/
PrecessionFunction(final PolynomialNutation psiA,
final PolynomialNutation omegaA,
final PolynomialNutation chiA,
final TimeScales timeScales) {
this.psiA = psiA;
this.omegaA = omegaA;
this.chiA = chiA;
this.timeScales = timeScales;
}
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
final double tc = evaluateTC(date, timeScales);
return new double[] {
psiA.value(tc), omegaA.value(tc), chiA.value(tc)
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
final T[] a = MathArrays.buildArray(date.getField(), 3);
final T tc = evaluateTC(date, timeScales);
a[0] = psiA.value(tc);
a[1] = omegaA.value(tc);
a[2] = chiA.value(tc);
return a;
}
}
/** Local class for tides frequency function. */
private static class TideFrequencyDependenceFunction implements TimeVectorFunction {
/** Nutation arguments. */
private final FundamentalNutationArguments arguments;
/** Correction series. */
private final PoissonSeries.CompiledSeries kSeries;
/** Simple constructor.
* @param arguments nutation arguments
* @param c20Series correction series for the C20 term
* @param c21Series correction series for the C21 term
* @param s21Series correction series for the S21 term
* @param c22Series correction series for the C22 term
* @param s22Series correction series for the S22 term
*/
TideFrequencyDependenceFunction(final FundamentalNutationArguments arguments,
final PoissonSeries c20Series,
final PoissonSeries c21Series,
final PoissonSeries s21Series,
final PoissonSeries c22Series,
final PoissonSeries s22Series) {
this.arguments = arguments;
this.kSeries = PoissonSeries.compile(c20Series, c21Series, s21Series, c22Series, s22Series);
}
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
return kSeries.value(arguments.evaluateAll(date));
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
return kSeries.value(arguments.evaluateAll(date));
}
}
/** Local class for EOP tidal corrections. */
private static class EOPTidalCorrection implements TimeVectorFunction {
/** Nutation arguments. */
private final FundamentalNutationArguments arguments;
/** Correction series. */
private final PoissonSeries.CompiledSeries correctionSeries;
/** Simple constructor.
* @param arguments nutation arguments
* @param xSeries correction series for the x coordinate
* @param ySeries correction series for the y coordinate
* @param ut1Series correction series for the UT1
*/
EOPTidalCorrection(final FundamentalNutationArguments arguments,
final PoissonSeries xSeries,
final PoissonSeries ySeries,
final PoissonSeries ut1Series) {
this.arguments = arguments;
this.correctionSeries = PoissonSeries.compile(xSeries, ySeries, ut1Series);
}
/** {@inheritDoc} */
@Override
public double[] value(final AbsoluteDate date) {
final BodiesElements elements = arguments.evaluateAll(date);
final double[] correction = correctionSeries.value(elements);
final double[] correctionDot = correctionSeries.derivative(elements);
return new double[] {
correction[0],
correction[1],
correction[2],
correctionDot[2] * (-Constants.JULIAN_DAY)
};
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
final FieldBodiesElements<T> elements = arguments.evaluateAll(date);
final T[] correction = correctionSeries.value(elements);
final T[] correctionDot = correctionSeries.derivative(elements);
final T[] a = MathArrays.buildArray(date.getField(), 4);
a[0] = correction[0];
a[1] = correction[1];
a[2] = correction[2];
a[3] = correctionDot[2].multiply(-Constants.JULIAN_DAY);
return a;
}
}
}