EOPHistory.java

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

  19. import java.util.ArrayList;
  20. import java.util.Collection;
  21. import java.util.List;
  22. import java.util.Optional;
  23. import java.util.function.BiFunction;
  24. import java.util.function.Consumer;
  25. import java.util.function.Function;
  26. import java.util.stream.Stream;

  28. import org.hipparchus.CalculusFieldElement;
  29. import org.hipparchus.analysis.interpolation.FieldHermiteInterpolator;
  30. import org.hipparchus.analysis.interpolation.HermiteInterpolator;
  31. import org.hipparchus.util.FastMath;
  32. import org.hipparchus.util.MathArrays;
  33. import org.orekit.annotation.DefaultDataContext;
  34. import org.orekit.data.DataContext;
  35. import org.orekit.errors.OrekitException;
  36. import org.orekit.errors.OrekitInternalError;
  37. import org.orekit.errors.OrekitMessages;
  38. import org.orekit.errors.TimeStampedCacheException;
  39. import org.orekit.time.AbsoluteDate;
  40. import org.orekit.time.FieldAbsoluteDate;
  41. import org.orekit.time.TimeScales;
  42. import org.orekit.time.TimeStamped;
  43. import org.orekit.time.TimeVectorFunction;
  44. import org.orekit.utils.Constants;
  45. import org.orekit.utils.GenericTimeStampedCache;
  46. import org.orekit.utils.IERSConventions;
  47. import org.orekit.utils.ImmutableTimeStampedCache;
  48. import org.orekit.utils.OrekitConfiguration;
  49. import org.orekit.utils.TimeStampedCache;
  50. import org.orekit.utils.TimeStampedGenerator;

  52. /** This class loads any kind of Earth Orientation Parameter data throughout a large time range.
  53.  * @author Pascal Parraud
  54.  * @author Evan Ward
  55.  */
  56. public class EOPHistory {

  58.     /** Default interpolation degree.
  59.      * @since 12.0
  60.      */
  61.     public static final int DEFAULT_INTERPOLATION_DEGREE = 3;

  63.     /** Interpolation degree.
  64.      * @since 12.0
  65.      */
  66.     private final int interpolationDegree;

  68.     /**
  69.      * If this history has any EOP data.
  70.      *
  71.      * @see #hasDataFor(AbsoluteDate)
  72.      */
  73.     private final boolean hasData;

  75.     /** EOP history entries. */
  76.     private final transient ImmutableTimeStampedCache<EOPEntry> cache;

  78.     /** IERS conventions to which EOP refers. */
  79.     private final IERSConventions conventions;

  81.     /** Correction to apply to EOP (may be null). */
  82.     private final transient TimeVectorFunction tidalCorrection;

  84.     /** Time scales to use when computing corrections. */
  85.     private final transient TimeScales timeScales;

  87.     /** Simple constructor.
  88.      *
  89.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  90.      *
  91.      * @param conventions IERS conventions to which EOP refers
  92.      * @param interpolationDegree interpolation degree (must be of the form 4k-1)
  93.      * @param data the EOP data to use
  94.      * @param simpleEOP if true, tidal effects are ignored when interpolating EOP
  95.      * @see #EOPHistory(IERSConventions, int, Collection, boolean, TimeScales)
  96.      */
  97.     @DefaultDataContext
  98.     protected EOPHistory(final IERSConventions conventions,
  99.                          final int interpolationDegree,
  100.                          final Collection<? extends EOPEntry> data,
  101.                          final boolean simpleEOP) {
  102.         this(conventions, interpolationDegree, data, simpleEOP, DataContext.getDefault().getTimeScales());
  103.     }

  105.     /** Simple constructor.
  106.      * @param conventions IERS conventions to which EOP refers
  107.      * @param interpolationDegree interpolation degree (must be of the form 4k-1)
  108.      * @param data the EOP data to use
  109.      * @param simpleEOP if true, tidal effects are ignored when interpolating EOP
  110.      * @param timeScales to use when computing EOP corrections.
  111.      * @since 10.1
  112.      */
  113.     public EOPHistory(final IERSConventions conventions,
  114.                       final int interpolationDegree,
  115.                       final Collection<? extends EOPEntry> data,
  116.                       final boolean simpleEOP,
  117.                       final TimeScales timeScales) {
  118.         this(conventions, interpolationDegree, data,
  119.              simpleEOP ? null : new CachedCorrection(conventions.getEOPTidalCorrection(timeScales)),
  120.              timeScales);
  121.     }

  123.     /** Simple constructor.
  124.      * @param conventions IERS conventions to which EOP refers
  125.      * @param interpolationDegree interpolation degree (must be of the form 4k-1)
  126.      * @param data the EOP data to use
  127.      * @param tidalCorrection correction to apply to EOP
  128.      * @param timeScales to use when computing EOP corrections
  129.      * @since 12
  130.      */
  131.     private EOPHistory(final IERSConventions conventions,
  132.                        final int interpolationDegree,
  133.                        final Collection<? extends EOPEntry> data,
  134.                        final TimeVectorFunction tidalCorrection,
  135.                        final TimeScales timeScales) {

  137.         // check interpolation degree is 4k-1
  138.         final int k = (interpolationDegree + 1) / 4;
  139.         if (interpolationDegree != 4 * k - 1) {
  140.             throw new OrekitException(OrekitMessages.WRONG_EOP_INTERPOLATION_DEGREE, interpolationDegree);
  141.         }

  143.         this.conventions         = conventions;
  144.         this.interpolationDegree = interpolationDegree;
  145.         this.tidalCorrection     = tidalCorrection;
  146.         this.timeScales          = timeScales;
  147.         if (data.size() >= 1) {
  148.             // enough data to interpolate
  149.             if (missSomeDerivatives(data)) {
  150.                 // we need to estimate the missing derivatives
  151.                 final ImmutableTimeStampedCache<EOPEntry> rawCache =
  152.                                 new ImmutableTimeStampedCache<>(FastMath.min(interpolationDegree + 1, data.size()), data);
  153.                 final List<EOPEntry>fixedData = new ArrayList<>();
  154.                 for (final EOPEntry entry : rawCache.getAll()) {
  155.                     fixedData.add(fixDerivatives(entry, rawCache));
  156.                 }
  157.                 cache = new ImmutableTimeStampedCache<>(FastMath.min(interpolationDegree + 1, fixedData.size()), fixedData);
  158.             } else {
  159.                 cache = new ImmutableTimeStampedCache<>(FastMath.min(interpolationDegree + 1, data.size()), data);
  160.             }
  161.             hasData = true;
  162.         } else {
  163.             // not enough data to interpolate -> always use null correction
  164.             cache   = ImmutableTimeStampedCache.emptyCache();
  165.             hasData = false;
  166.         }
  167.     }

  169.     /**
  170.      * Determine if this history uses simplified EOP corrections.
  171.      *
  172.      * @return {@code true} if tidal corrections are ignored, {@code false} otherwise.
  173.      */
  174.     public boolean isSimpleEop() {
  175.         return tidalCorrection == null;
  176.     }

  178.     /** Get interpolation degree.
  179.      * @return interpolation degree
  180.      * @since 12.0
  181.      */
  182.     public int getInterpolationDegree() {
  183.         return interpolationDegree;
  184.     }

  186.     /**
  187.      * Get the time scales used in computing EOP corrections.
  188.      *
  189.      * @return set of time scales.
  190.      * @since 10.1
  191.      */
  192.     public TimeScales getTimeScales() {
  193.         return timeScales;
  194.     }

  196.     /** Get version of the instance that does not cache tidal correction.
  197.      * @return version of the instance that does not cache tidal correction
  198.      * @since 12.0
  199.      */
  200.     public EOPHistory getEOPHistoryWithoutCachedTidalCorrection() {
  201.         return new EOPHistory(conventions, interpolationDegree, getEntries(),
  202.                               conventions.getEOPTidalCorrection(timeScales),
  203.                               timeScales);
  204.     }

  206.     /** Check if the instance caches tidal corrections.
  207.      * @return true if the instance caches tidal corrections
  208.      * @since 12.0
  209.      */
  210.     public boolean cachesTidalCorrection() {
  211.         return tidalCorrection instanceof CachedCorrection;
  212.     }

  214.     /** Get the IERS conventions to which these EOP apply.
  215.      * @return IERS conventions to which these EOP apply
  216.      */
  217.     public IERSConventions getConventions() {
  218.         return conventions;
  219.     }

  221.     /** Get the date of the first available Earth Orientation Parameters.
  222.      * @return the start date of the available data
  223.      */
  224.     public AbsoluteDate getStartDate() {
  225.         return this.cache.getEarliest().getDate();
  226.     }

  228.     /** Get the date of the last available Earth Orientation Parameters.
  229.      * @return the end date of the available data
  230.      */
  231.     public AbsoluteDate getEndDate() {
  232.         return this.cache.getLatest().getDate();
  233.     }

  235.     /** Get the UT1-UTC value.
  236.      * <p>The data provided comes from the IERS files. It is smoothed data.</p>
  237.      * @param date date at which the value is desired
  238.      * @return UT1-UTC in seconds (0 if date is outside covered range)
  239.      */
  240.     public double getUT1MinusUTC(final AbsoluteDate date) {

  242.         //check if there is data for date
  243.         if (!this.hasDataFor(date)) {
  244.             // no EOP data available for this date, we use a default 0.0 offset
  245.             return (tidalCorrection == null) ? 0.0 : tidalCorrection.value(date)[2];
  246.         }

  248.         // we have EOP data -> interpolate offset
  249.         try {
  250.             final DUT1Interpolator interpolator = new DUT1Interpolator(date);
  251.             getNeighbors(date, interpolationDegree + 1).forEach(interpolator);
  252.             double interpolated = interpolator.getInterpolated();
  253.             if (tidalCorrection != null) {
  254.                 interpolated += tidalCorrection.value(date)[2];
  255.             }
  256.             return interpolated;
  257.         } catch (TimeStampedCacheException tce) {
  258.             //this should not happen because of date check above
  259.             throw new OrekitInternalError(tce);
  260.         }

  262.     }

  264.     /** Get the UT1-UTC value.
  265.      * <p>The data provided comes from the IERS files. It is smoothed data.</p>
  266.      * @param date date at which the value is desired
  267.      * @param <T> type of the field elements
  268.      * @return UT1-UTC in seconds (0 if date is outside covered range)
  269.      * @since 9.0
  270.      */
  271.     public <T extends CalculusFieldElement<T>> T getUT1MinusUTC(final FieldAbsoluteDate<T> date) {

  273.         //check if there is data for date
  274.         final AbsoluteDate absDate = date.toAbsoluteDate();
  275.         if (!this.hasDataFor(absDate)) {
  276.             // no EOP data available for this date, we use a default 0.0 offset
  277.             return (tidalCorrection == null) ? date.getField().getZero() : tidalCorrection.value(date)[2];
  278.         }

  280.         // we have EOP data -> interpolate offset
  281.         try {
  282.             final FieldDUT1Interpolator<T> interpolator = new FieldDUT1Interpolator<>(date, absDate);
  283.             getNeighbors(absDate, interpolationDegree + 1).forEach(interpolator);
  284.             T interpolated = interpolator.getInterpolated();
  285.             if (tidalCorrection != null) {
  286.                 interpolated = interpolated.add(tidalCorrection.value(date)[2]);
  287.             }
  288.             return interpolated;
  289.         } catch (TimeStampedCacheException tce) {
  290.             //this should not happen because of date check above
  291.             throw new OrekitInternalError(tce);
  292.         }

  294.     }

  296.     /** Local class for DUT1 interpolation, crossing leaps safely. */
  297.     private static class DUT1Interpolator implements Consumer<EOPEntry> {

  299.         /** DUT at first entry. */
  300.         private double firstDUT;

  302.         /** Indicator for dates just before a leap occurring during the interpolation sample. */
  303.         private boolean beforeLeap;

  305.         /** Interpolator to use. */
  306.         private final HermiteInterpolator interpolator;

  308.         /** Interpolation date. */
  309.         private AbsoluteDate date;

  311.         /** Simple constructor.
  312.          * @param date interpolation date
  313.          */
  314.         DUT1Interpolator(final AbsoluteDate date) {
  315.             this.firstDUT     = Double.NaN;
  316.             this.beforeLeap   = true;
  317.             this.interpolator = new HermiteInterpolator();
  318.             this.date         = date;
  319.         }

  321.         /** {@inheritDoc} */
  322.         @Override
  323.         public void accept(final EOPEntry neighbor) {
  324.             if (Double.isNaN(firstDUT)) {
  325.                 firstDUT = neighbor.getUT1MinusUTC();
  326.             }
  327.             final double dut;
  328.             if (neighbor.getUT1MinusUTC() - firstDUT > 0.9) {
  329.                 // there was a leap second between the entries
  330.                 dut = neighbor.getUT1MinusUTC() - 1.0;
  331.                 // UTCScale considers the discontinuity to occur at the start of the leap
  332.                 // second so this code must use the same convention. EOP entries are time
  333.                 // stamped at midnight UTC so 1 second before is the start of the leap
  334.                 // second.
  335.                 if (neighbor.getDate().shiftedBy(-1).compareTo(date) <= 0) {
  336.                     beforeLeap = false;
  337.                 }
  338.             } else {
  339.                 dut = neighbor.getUT1MinusUTC();
  340.             }
  341.             interpolator.addSamplePoint(neighbor.getDate().durationFrom(date),
  342.                                         new double[] {
  343.                                             dut
  344.                                         });
  345.         }

  347.         /** Get the interpolated value.
  348.          * @return interpolated value
  349.          */
  350.         public double getInterpolated() {
  351.             final double interpolated = interpolator.value(0)[0];
  352.             return beforeLeap ? interpolated : interpolated + 1.0;
  353.         }

  355.     }

  357.     /** Local class for DUT1 interpolation, crossing leaps safely. */
  358.     private static class FieldDUT1Interpolator<T extends CalculusFieldElement<T>> implements Consumer<EOPEntry> {

  360.         /** DUT at first entry. */
  361.         private double firstDUT;

  363.         /** Indicator for dates just before a leap occurring during the interpolation sample. */
  364.         private boolean beforeLeap;

  366.         /** Interpolator to use. */
  367.         private final FieldHermiteInterpolator<T> interpolator;

  369.         /** Interpolation date. */
  370.         private FieldAbsoluteDate<T> date;

  372.         /** Interpolation date. */
  373.         private AbsoluteDate absDate;

  375.         /** Simple constructor.
  376.          * @param date interpolation date
  377.          * @param absDate interpolation date
  378.          */
  379.         FieldDUT1Interpolator(final FieldAbsoluteDate<T> date, final AbsoluteDate absDate) {
  380.             this.firstDUT     = Double.NaN;
  381.             this.beforeLeap   = true;
  382.             this.interpolator = new FieldHermiteInterpolator<>();
  383.             this.date         = date;
  384.             this.absDate      = absDate;
  385.         }

  387.         /** {@inheritDoc} */
  388.         @Override
  389.         public void accept(final EOPEntry neighbor) {
  390.             if (Double.isNaN(firstDUT)) {
  391.                 firstDUT = neighbor.getUT1MinusUTC();
  392.             }
  393.             final double dut;
  394.             if (neighbor.getUT1MinusUTC() - firstDUT > 0.9) {
  395.                 // there was a leap second between the entries
  396.                 dut = neighbor.getUT1MinusUTC() - 1.0;
  397.                 if (neighbor.getDate().compareTo(absDate) <= 0) {
  398.                     beforeLeap = false;
  399.                 }
  400.             } else {
  401.                 dut = neighbor.getUT1MinusUTC();
  402.             }
  403.             final T[] array = MathArrays.buildArray(date.getField(), 1);
  404.             array[0] = date.getField().getZero().newInstance(dut);
  405.             interpolator.addSamplePoint(date.durationFrom(neighbor.getDate()).negate(),
  406.                                         array);
  407.         }

  409.         /** Get the interpolated value.
  410.          * @return interpolated value
  411.          */
  412.         public T getInterpolated() {
  413.             final T interpolated = interpolator.value(date.getField().getZero())[0];
  414.             return beforeLeap ? interpolated : interpolated.add(1.0);
  415.         }

  417.     }

  419.     /**
  420.      * Get the entries surrounding a central date.
  421.      * <p>
  422.      * See {@link #hasDataFor(AbsoluteDate)} to determine if the cache has data
  423.      * for {@code central} without throwing an exception.
  424.      *
  425.      * @param central central date
  426.      * @param n number of neighbors
  427.      * @return array of cached entries surrounding specified date
  428.      * @since 12.0
  429.      */
  430.     protected Stream<EOPEntry> getNeighbors(final AbsoluteDate central, final int n) {
  431.         return cache.getNeighbors(central, n);
  432.     }

  434.     /** Get the LoD (Length of Day) value.
  435.      * <p>The data provided comes from the IERS files. It is smoothed data.</p>
  436.      * @param date date at which the value is desired
  437.      * @return LoD in seconds (0 if date is outside covered range)
  438.      */
  439.     public double getLOD(final AbsoluteDate date) {

  441.         // check if there is data for date
  442.         if (!this.hasDataFor(date)) {
  443.             // no EOP data available for this date, we use a default null correction
  444.             return (tidalCorrection == null) ? 0.0 : tidalCorrection.value(date)[3];
  445.         }

  447.         // we have EOP data for date -> interpolate correction
  448.         double interpolated = interpolate(date, EOPEntry::getLOD);
  449.         if (tidalCorrection != null) {
  450.             interpolated += tidalCorrection.value(date)[3];
  451.         }
  452.         return interpolated;

  454.     }

  456.     /** Get the LoD (Length of Day) value.
  457.      * <p>The data provided comes from the IERS files. It is smoothed data.</p>
  458.      * @param date date at which the value is desired
  459.      * @param <T> type of the filed elements
  460.      * @return LoD in seconds (0 if date is outside covered range)
  461.      * @since 9.0
  462.      */
  463.     public <T extends CalculusFieldElement<T>> T getLOD(final FieldAbsoluteDate<T> date) {

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

  467.         // check if there is data for date
  468.         if (!this.hasDataFor(aDate)) {
  469.             // no EOP data available for this date, we use a default null correction
  470.             return (tidalCorrection == null) ? date.getField().getZero() : tidalCorrection.value(date)[3];
  471.         }

  473.         // we have EOP data for date -> interpolate correction
  474.         T interpolated = interpolate(date, aDate, EOPEntry::getLOD);
  475.         if (tidalCorrection != null) {
  476.             interpolated = interpolated.add(tidalCorrection.value(date)[3]);
  477.         }

  479.         return interpolated;

  481.     }

  483.     /** Get the pole IERS Reference Pole correction.
  484.      * <p>The data provided comes from the IERS files. It is smoothed data.</p>
  485.      * @param date date at which the correction is desired
  486.      * @return pole correction ({@link PoleCorrection#NULL_CORRECTION
  487.      * PoleCorrection.NULL_CORRECTION} if date is outside covered range)
  488.      */
  489.     public PoleCorrection getPoleCorrection(final AbsoluteDate date) {

  491.         // check if there is data for date
  492.         if (!this.hasDataFor(date)) {
  493.             // no EOP data available for this date, we use a default null correction
  494.             if (tidalCorrection == null) {
  495.                 return PoleCorrection.NULL_CORRECTION;
  496.             } else {
  497.                 final double[] correction = tidalCorrection.value(date);
  498.                 return new PoleCorrection(correction[0], correction[1]);
  499.             }
  500.         }

  502.         // we have EOP data for date -> interpolate correction
  503.         final double[] interpolated = interpolate(date,
  504.                 EOPEntry::getX, EOPEntry::getY,
  505.                 EOPEntry::getXRate, EOPEntry::getYRate);
  506.         if (tidalCorrection != null) {
  507.             final double[] correction = tidalCorrection.value(date);
  508.             interpolated[0] += correction[0];
  509.             interpolated[1] += correction[1];
  510.         }
  511.         return new PoleCorrection(interpolated[0], interpolated[1]);

  513.     }

  515.     /** Get the pole IERS Reference Pole correction.
  516.      * <p>The data provided comes from the IERS files. It is smoothed data.</p>
  517.      * @param date date at which the correction is desired
  518.      * @param <T> type of the field elements
  519.      * @return pole correction ({@link PoleCorrection#NULL_CORRECTION
  520.      * PoleCorrection.NULL_CORRECTION} if date is outside covered range)
  521.      */
  522.     public <T extends CalculusFieldElement<T>> FieldPoleCorrection<T> getPoleCorrection(final FieldAbsoluteDate<T> date) {

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

  526.         // check if there is data for date
  527.         if (!this.hasDataFor(aDate)) {
  528.             // no EOP data available for this date, we use a default null correction
  529.             if (tidalCorrection == null) {
  530.                 return new FieldPoleCorrection<>(date.getField().getZero(), date.getField().getZero());
  531.             } else {
  532.                 final T[] correction = tidalCorrection.value(date);
  533.                 return new FieldPoleCorrection<>(correction[0], correction[1]);
  534.             }
  535.         }

  537.         // we have EOP data for date -> interpolate correction
  538.         final T[] interpolated = interpolate(date, aDate, EOPEntry::getX, EOPEntry::getY);
  539.         if (tidalCorrection != null) {
  540.             final T[] correction = tidalCorrection.value(date);
  541.             interpolated[0] = interpolated[0].add(correction[0]);
  542.             interpolated[1] = interpolated[1].add(correction[1]);
  543.         }
  544.         return new FieldPoleCorrection<>(interpolated[0], interpolated[1]);

  546.     }

  548.     /** Get the correction to the nutation parameters for equinox-based paradigm.
  549.      * <p>The data provided comes from the IERS files. It is smoothed data.</p>
  550.      * @param date date at which the correction is desired
  551.      * @return nutation correction in longitude ΔΨ and in obliquity Δε
  552.      * (zero if date is outside covered range)
  553.      */
  554.     public double[] getEquinoxNutationCorrection(final AbsoluteDate date) {

  556.         // check if there is data for date
  557.         if (!this.hasDataFor(date)) {
  558.             // no EOP data available for this date, we use a default null correction
  559.             return new double[2];
  560.         }

  562.         // we have EOP data for date -> interpolate correction
  563.         return interpolate(date, EOPEntry::getDdPsi, EOPEntry::getDdEps);

  565.     }

  567.     /** Get the correction to the nutation parameters for equinox-based paradigm.
  568.      * <p>The data provided comes from the IERS files. It is smoothed data.</p>
  569.      * @param date date at which the correction is desired
  570.      * @param <T> type of the field elements
  571.      * @return nutation correction in longitude ΔΨ and in obliquity Δε
  572.      * (zero if date is outside covered range)
  573.      */
  574.     public <T extends CalculusFieldElement<T>> T[] getEquinoxNutationCorrection(final FieldAbsoluteDate<T> date) {

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

  578.         // check if there is data for date
  579.         if (!this.hasDataFor(aDate)) {
  580.             // no EOP data available for this date, we use a default null correction
  581.             return MathArrays.buildArray(date.getField(), 2);
  582.         }

  584.         // we have EOP data for date -> interpolate correction
  585.         return interpolate(date, aDate, EOPEntry::getDdPsi, EOPEntry::getDdEps);

  587.     }

  589.     /** Get the correction to the nutation parameters for Non-Rotating Origin paradigm.
  590.      * <p>The data provided comes from the IERS files. It is smoothed data.</p>
  591.      * @param date date at which the correction is desired
  592.      * @return nutation correction in Celestial Intermediate Pole coordinates
  593.      * δX and δY (zero if date is outside covered range)
  594.      */
  595.     public double[] getNonRotatinOriginNutationCorrection(final AbsoluteDate date) {

  597.         // check if there is data for date
  598.         if (!this.hasDataFor(date)) {
  599.             // no EOP data available for this date, we use a default null correction
  600.             return new double[2];
  601.         }

  603.         // we have EOP data for date -> interpolate correction
  604.         return interpolate(date, EOPEntry::getDx, EOPEntry::getDy);

  606.     }

  608.     /** Get the correction to the nutation parameters for Non-Rotating Origin paradigm.
  609.      * <p>The data provided comes from the IERS files. It is smoothed data.</p>
  610.      * @param date date at which the correction is desired
  611.      * @param <T> type of the filed elements
  612.      * @return nutation correction in Celestial Intermediate Pole coordinates
  613.      * δX and δY (zero if date is outside covered range)
  614.      */
  615.     public <T extends CalculusFieldElement<T>> T[] getNonRotatinOriginNutationCorrection(final FieldAbsoluteDate<T> date) {

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

  619.         // check if there is data for date
  620.         if (!this.hasDataFor(aDate)) {
  621.             // no EOP data available for this date, we use a default null correction
  622.             return MathArrays.buildArray(date.getField(), 2);
  623.         }

  625.         // we have EOP data for date -> interpolate correction
  626.         return interpolate(date, aDate, EOPEntry::getDx, EOPEntry::getDy);

  628.     }

  630.     /** Get the ITRF version.
  631.      * @param date date at which the value is desired
  632.      * @return ITRF version of the EOP covering the specified date
  633.      * @since 9.2
  634.      */
  635.     public ITRFVersion getITRFVersion(final AbsoluteDate date) {

  637.         // check if there is data for date
  638.         if (!this.hasDataFor(date)) {
  639.             // no EOP data available for this date, we use a default ITRF 2014
  640.             return ITRFVersion.ITRF_2014;
  641.         }

  643.         try {
  644.             // we have EOP data for date
  645.             final Optional<EOPEntry> first = getNeighbors(date, 1).findFirst();
  646.             return first.isPresent() ? first.get().getITRFType() : ITRFVersion.ITRF_2014;

  648.         } catch (TimeStampedCacheException tce) {
  649.             // this should not happen because of date check performed at start
  650.             throw new OrekitInternalError(tce);
  651.         }

  653.     }

  655.     /** Check Earth orientation parameters continuity.
  656.      * @param maxGap maximal allowed gap between entries (in seconds)
  657.      */
  658.     public void checkEOPContinuity(final double maxGap) {
  659.         TimeStamped preceding = null;
  660.         for (final TimeStamped current : this.cache.getAll()) {

  662.             // compare the dates of preceding and current entries
  663.             if (preceding != null && (current.getDate().durationFrom(preceding.getDate())) > maxGap) {
  664.                 throw new OrekitException(OrekitMessages.MISSING_EARTH_ORIENTATION_PARAMETERS_BETWEEN_DATES_GAP,
  665.                                           preceding.getDate(), current.getDate(),
  666.                                           current.getDate().durationFrom(preceding.getDate()));
  667.             }

  669.             // prepare next iteration
  670.             preceding = current;

  672.         }
  673.     }

  675.     /**
  676.      * Check if the cache has data for the given date using
  677.      * {@link #getStartDate()} and {@link #getEndDate()}.
  678.      *
  679.      * @param date the requested date
  680.      * @return true if the {@link #cache} has data for the requested date, false
  681.      *         otherwise.
  682.      */
  683.     protected boolean hasDataFor(final AbsoluteDate date) {
  684.         /*
  685.          * when there is no EOP data, short circuit getStartDate, which will
  686.          * throw an exception
  687.          */
  688.         return this.hasData && this.getStartDate().compareTo(date) <= 0 &&
  689.                date.compareTo(this.getEndDate()) <= 0;
  690.     }

  692.     /** Get a non-modifiable view of the EOP entries.
  693.      * @return non-modifiable view of the EOP entries
  694.      */
  695.     public List<EOPEntry> getEntries() {
  696.         return cache.getAll();
  697.     }

  699.     /** Interpolate a single EOP component.
  700.      * <p>
  701.      * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
  702.      * </p>
  703.      * @param date interpolation date
  704.      * @param selector selector for EOP entry component
  705.      * @return interpolated value
  706.      */
  707.     private double interpolate(final AbsoluteDate date, final Function<EOPEntry, Double> selector) {
  708.         try {
  709.             final HermiteInterpolator interpolator = new HermiteInterpolator();
  710.             getNeighbors(date, interpolationDegree + 1).
  711.                 forEach(entry -> interpolator.addSamplePoint(entry.getDate().durationFrom(date),
  712.                                                              new double[] {
  713.                                                                  selector.apply(entry)
  714.                                                              }));
  715.             return interpolator.value(0)[0];
  716.         } catch (TimeStampedCacheException tce) {
  717.             // this should not happen because of date check performed by caller
  718.             throw new OrekitInternalError(tce);
  719.         }
  720.     }

  722.     /** Interpolate a single EOP component.
  723.      * <p>
  724.      * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
  725.      * </p>
  726.      * @param date interpolation date
  727.      * @param aDate interpolation date, as an {@link AbsoluteDate}
  728.      * @param selector selector for EOP entry component
  729.      * @param <T> type of the field elements
  730.      * @return interpolated value
  731.      */
  732.     private <T extends CalculusFieldElement<T>> T interpolate(final FieldAbsoluteDate<T> date,
  733.                                                               final AbsoluteDate aDate,
  734.                                                               final Function<EOPEntry, Double> selector) {
  735.         try {
  736.             final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
  737.             final T[] y = MathArrays.buildArray(date.getField(), 1);
  738.             final T zero = date.getField().getZero();
  739.             final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
  740.                                                                                        // for example removing derivatives
  741.                                                                                        // if T was DerivativeStructure
  742.             getNeighbors(aDate, interpolationDegree + 1).
  743.                 forEach(entry -> {
  744.                     y[0] = zero.newInstance(selector.apply(entry));
  745.                     interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
  746.                 });
  747.             return interpolator.value(date.durationFrom(central))[0]; // here, we introduce derivatives again (in DerivativeStructure case)
  748.         } catch (TimeStampedCacheException tce) {
  749.             // this should not happen because of date check performed by caller
  750.             throw new OrekitInternalError(tce);
  751.         }
  752.     }

  754.     /** Interpolate two EOP components.
  755.      * <p>
  756.      * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
  757.      * </p>
  758.      * @param date interpolation date
  759.      * @param selector1 selector for first EOP entry component
  760.      * @param selector2 selector for second EOP entry component
  761.      * @return interpolated value
  762.      */
  763.     private double[] interpolate(final AbsoluteDate date,
  764.                                  final Function<EOPEntry, Double> selector1,
  765.                                  final Function<EOPEntry, Double> selector2) {
  766.         try {
  767.             final HermiteInterpolator interpolator = new HermiteInterpolator();
  768.             getNeighbors(date, interpolationDegree + 1).
  769.                 forEach(entry -> interpolator.addSamplePoint(entry.getDate().durationFrom(date),
  770.                                                              new double[] {
  771.                                                                  selector1.apply(entry),
  772.                                                                  selector2.apply(entry)
  773.                                                              }));
  774.             return interpolator.value(0);
  775.         } catch (TimeStampedCacheException tce) {
  776.             // this should not happen because of date check performed by caller
  777.             throw new OrekitInternalError(tce);
  778.         }
  779.     }

  781.     /** Interpolate two EOP components.
  782.      * <p>
  783.      * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
  784.      * </p>
  785.      * @param date interpolation date
  786.      * @param selector1 selector for first EOP entry component
  787.      * @param selector2 selector for second EOP entry component
  788.      * @param selector1Rate selector for first EOP entry component rate
  789.      * @param selector2Rate selector for second EOP entry component rate
  790.      * @return interpolated value
  791.      * @since 12.0
  792.      */
  793.     private double[] interpolate(final AbsoluteDate date,
  794.                                  final Function<EOPEntry, Double> selector1,
  795.                                  final Function<EOPEntry, Double> selector2,
  796.                                  final Function<EOPEntry, Double> selector1Rate,
  797.                                  final Function<EOPEntry, Double> selector2Rate) {
  798.         try {
  799.             final HermiteInterpolator interpolator = new HermiteInterpolator();
  800.             getNeighbors(date, (interpolationDegree + 1) / 2).
  801.                 forEach(entry -> interpolator.addSamplePoint(entry.getDate().durationFrom(date),
  802.                                                              new double[] {
  803.                                                                  selector1.apply(entry),
  804.                                                                  selector2.apply(entry)
  805.                                                              },
  806.                                                              new double[] {
  807.                                                                  selector1Rate.apply(entry),
  808.                                                                  selector2Rate.apply(entry)
  809.                                                              }));
  810.             return interpolator.value(0);
  811.         } catch (TimeStampedCacheException tce) {
  812.             // this should not happen because of date check performed by caller
  813.             throw new OrekitInternalError(tce);
  814.         }
  815.     }

  817.     /** Interpolate two EOP components.
  818.      * <p>
  819.      * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
  820.      * </p>
  821.      * @param date interpolation date
  822.      * @param aDate interpolation date, as an {@link AbsoluteDate}
  823.      * @param selector1 selector for first EOP entry component
  824.      * @param selector2 selector for second EOP entry component
  825.      * @param <T> type of the field elements
  826.      * @return interpolated value
  827.      */
  828.     private <T extends CalculusFieldElement<T>> T[] interpolate(final FieldAbsoluteDate<T> date,
  829.                                                                 final AbsoluteDate aDate,
  830.                                                                 final Function<EOPEntry, Double> selector1,
  831.                                                                 final Function<EOPEntry, Double> selector2) {
  832.         try {
  833.             final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
  834.             final T[] y = MathArrays.buildArray(date.getField(), 2);
  835.             final T zero = date.getField().getZero();
  836.             final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
  837.                                                                                        // for example removing derivatives
  838.                                                                                        // if T was DerivativeStructure
  839.             getNeighbors(aDate, interpolationDegree + 1).
  840.                 forEach(entry -> {
  841.                     y[0] = zero.newInstance(selector1.apply(entry));
  842.                     y[1] = zero.newInstance(selector2.apply(entry));
  843.                     interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
  844.                 });
  845.             return interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
  846.         } catch (TimeStampedCacheException tce) {
  847.             // this should not happen because of date check performed by caller
  848.             throw new OrekitInternalError(tce);
  849.         }
  850.     }

  852.     /** Check if some derivatives are missing.
  853.      * @param raw raw EOP history
  854.      * @return true if some derivatives are missing
  855.      * @since 12.0
  856.      */
  857.     private boolean missSomeDerivatives(final Collection<? extends EOPEntry> raw) {
  858.         for (final EOPEntry entry : raw) {
  859.             if (Double.isNaN(entry.getLOD() + entry.getXRate() + entry.getYRate())) {
  860.                 return true;
  861.             }
  862.         }
  863.         return false;
  864.     }

  866.     /** Fix missing derivatives.
  867.      * @param entry EOP entry to fix
  868.      * @param rawCache raw EOP history cache
  869.      * @return fixed entry
  870.      * @since 12.0
  871.      */
  872.     private EOPEntry fixDerivatives(final EOPEntry entry, final ImmutableTimeStampedCache<EOPEntry> rawCache) {

  874.         // helper function to differentiate some EOP parameters
  875.         final BiFunction<EOPEntry, Function<EOPEntry, Double>, Double> differentiator =
  876.                         (e, selector) -> {
  877.                             final HermiteInterpolator interpolator = new HermiteInterpolator();
  878.                             rawCache.getNeighbors(e.getDate()).
  879.                                 forEach(f -> interpolator.addSamplePoint(f.getDate().durationFrom(e.getDate()),
  880.                                                                          new double[] {
  881.                                                                              selector.apply(f)
  882.                                                                          }));
  883.                             return interpolator.derivatives(0.0, 1)[1][0];
  884.                         };

  886.         if (Double.isNaN(entry.getLOD() + entry.getXRate() + entry.getYRate())) {
  887.             final double lod   = Double.isNaN(entry.getLOD()) ?
  888.                                  -differentiator.apply(entry, EOPEntry::getUT1MinusUTC) :
  889.                                  entry.getLOD();
  890.             final double xRate = Double.isNaN(entry.getXRate()) ?
  891.                                  differentiator.apply(entry, EOPEntry::getX) :
  892.                                  entry.getXRate();
  893.             final double yRate = Double.isNaN(entry.getYRate()) ?
  894.                                  differentiator.apply(entry, EOPEntry::getY) :
  895.                                  entry.getYRate();
  896.             return new EOPEntry(entry.getMjd(),
  897.                                 entry.getUT1MinusUTC(), lod,
  898.                                 entry.getX(), entry.getY(), xRate, yRate,
  899.                                 entry.getDdPsi(), entry.getDdEps(),
  900.                                 entry.getDx(), entry.getDy(),
  901.                                 entry.getITRFType(), entry.getDate());
  902.         } else {
  903.             // the entry already has all derivatives
  904.             return entry;
  905.         }
  906.     }

  908.     /** Internal class for caching tidal correction. */
  909.     private static class TidalCorrectionEntry implements TimeStamped {

  911.         /** Entry date. */
  912.         private final AbsoluteDate date;

  914.         /** Correction. */
  915.         private final double[] correction;

  917.         /** Simple constructor.
  918.          * @param date entry date
  919.          * @param correction correction on the EOP parameters (xp, yp, ut1, lod)
  920.          */
  921.         TidalCorrectionEntry(final AbsoluteDate date, final double[] correction) {
  922.             this.date       = date;
  923.             this.correction = correction;
  924.         }

  926.         /** {@inheritDoc} */
  927.         @Override
  928.         public AbsoluteDate getDate() {
  929.             return date;
  930.         }

  932.     }

  934.     /** Local generator for thread-safe cache. */
  935.     private static class CachedCorrection
  936.         implements TimeVectorFunction, TimeStampedGenerator<TidalCorrectionEntry> {

  938.         /** Correction to apply to EOP (may be null). */
  939.         private final TimeVectorFunction tidalCorrection;

  941.         /** Step between generated entries. */
  942.         private final double step;

  944.         /** Tidal corrections entries cache. */
  945.         private final TimeStampedCache<TidalCorrectionEntry> cache;

  947.         /** Simple constructor.
  948.          * @param tidalCorrection function computing the tidal correction
  949.          */
  950.         CachedCorrection(final TimeVectorFunction tidalCorrection) {
  951.             this.step            = 60 * 60;
  952.             this.tidalCorrection = tidalCorrection;
  953.             this.cache           =
  954.                     new GenericTimeStampedCache<>(8,
  955.                             OrekitConfiguration.getCacheSlotsNumber(),
  956.                             Constants.JULIAN_DAY * 30,
  957.                             Constants.JULIAN_DAY,
  958.                             this);
  959.         }

  961.         /** {@inheritDoc} */
  962.         @Override
  963.         public double[] value(final AbsoluteDate date) {
  964.             try {
  965.                 // set up an interpolator
  966.                 final HermiteInterpolator interpolator = new HermiteInterpolator();
  967.                 cache.getNeighbors(date).forEach(entry -> interpolator.addSamplePoint(entry.date.durationFrom(date), entry.correction));

  969.                 // interpolate to specified date
  970.                 return interpolator.value(0.0);
  971.             } catch (TimeStampedCacheException tsce) {
  972.                 // this should never happen
  973.                 throw new OrekitInternalError(tsce);
  974.             }
  975.         }

  977.         /** {@inheritDoc} */
  978.         @Override
  979.         public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  980.             try {

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

  984.                 final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
  985.                 final T[] y = MathArrays.buildArray(date.getField(), 4);
  986.                 final T zero = date.getField().getZero();
  987.                 final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
  988.                                                                                            // for example removing derivatives
  989.                                                                                            // if T was DerivativeStructure
  990.                 cache.getNeighbors(aDate).forEach(entry -> {
  991.                     for (int i = 0; i < y.length; ++i) {
  992.                         y[i] = zero.newInstance(entry.correction[i]);
  993.                     }
  994.                     interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
  995.                 });

  997.                 // interpolate to specified date
  998.                 return interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)

  1000.             } catch (TimeStampedCacheException tsce) {
  1001.                 // this should never happen
  1002.                 throw new OrekitInternalError(tsce);
  1003.             }
  1004.         }

  1006.         /** {@inheritDoc} */
  1007.         @Override
  1008.         public List<TidalCorrectionEntry> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {

  1010.             final List<TidalCorrectionEntry> generated = new ArrayList<>();

  1012.             if (existingDate == null) {

  1014.                 // no prior existing entries, just generate a first set
  1015.                 for (int i = -cache.getMaxNeighborsSize() / 2; generated.size() < cache.getMaxNeighborsSize(); ++i) {
  1016.                     final AbsoluteDate t = date.shiftedBy(i * step);
  1017.                     generated.add(new TidalCorrectionEntry(t, tidalCorrection.value(t)));
  1018.                 }

  1020.             } else {

  1022.                 // some entries have already been generated
  1023.                 // add the missing ones up to specified date

  1025.                 AbsoluteDate t = existingDate;
  1026.                 if (date.compareTo(t) > 0) {
  1027.                     // forward generation
  1028.                     do {
  1029.                         t = t.shiftedBy(step);
  1030.                         generated.add(new TidalCorrectionEntry(t, tidalCorrection.value(t)));
  1031.                     } while (t.compareTo(date) <= 0);
  1032.                 } else {
  1033.                     // backward generation
  1034.                     do {
  1035.                         t = t.shiftedBy(-step);
  1036.                         generated.add(0, new TidalCorrectionEntry(t, tidalCorrection.value(t)));
  1037.                     } while (t.compareTo(date) >= 0);
  1038.                 }
  1039.             }

  1041.             // return the generated transforms
  1042.             return generated;

  1044.         }
  1045.     }

  1047. }
Created with JaCoCo 0.8.12.202403310830