EOPHistory.java

  1. /* Copyright 2002-2024 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.io.Serializable;
  20. import java.util.ArrayList;
  21. import java.util.Collection;
  22. import java.util.List;
  23. import java.util.Optional;
  24. import java.util.function.BiFunction;
  25. import java.util.function.Consumer;
  26. import java.util.function.Function;
  27. import java.util.stream.Stream;

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

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

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

  64.     /** Serializable UID. */
  65.     private static final long serialVersionUID = 20231022L;

  67.     /** Interpolation degree.
  68.      * @since 12.0
  69.      */
  70.     private final int interpolationDegree;

  72.     /**
  73.      * If this history has any EOP data.
  74.      *
  75.      * @see #hasDataFor(AbsoluteDate)
  76.      */
  77.     private final boolean hasData;

  79.     /** EOP history entries. */
  80.     private final transient ImmutableTimeStampedCache<EOPEntry> cache;

  82.     /** IERS conventions to which EOP refers. */
  83.     private final IERSConventions conventions;

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

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

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

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

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

  141.         // check interpolation degree is 4k-1
  142.         final int k = (interpolationDegree + 1) / 4;
  143.         if (interpolationDegree != 4 * k - 1) {
  144.             throw new OrekitException(OrekitMessages.WRONG_EOP_INTERPOLATION_DEGREE, interpolationDegree);
  145.         }

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

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

  182.     /** Get interpolation degree.
  183.      * @return interpolation degree
  184.      * @since 12.0
  185.      */
  186.     public int getInterpolationDegree() {
  187.         return interpolationDegree;
  188.     }

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

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

  210.     /** Check if the instance caches tidal corrections.
  211.      * @return true if the instance caches tidal corrections
  212.      * @since 12.0
  213.      */
  214.     public boolean cachesTidalCorrection() {
  215.         return tidalCorrection instanceof CachedCorrection;
  216.     }

  218.     /** Get the IERS conventions to which these EOP apply.
  219.      * @return IERS conventions to which these EOP apply
  220.      */
  221.     public IERSConventions getConventions() {
  222.         return conventions;
  223.     }

  225.     /** Get the date of the first available Earth Orientation Parameters.
  226.      * @return the start date of the available data
  227.      */
  228.     public AbsoluteDate getStartDate() {
  229.         return this.cache.getEarliest().getDate();
  230.     }

  232.     /** Get the date of the last available Earth Orientation Parameters.
  233.      * @return the end date of the available data
  234.      */
  235.     public AbsoluteDate getEndDate() {
  236.         return this.cache.getLatest().getDate();
  237.     }

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

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

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

  266.     }

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

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

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

  298.     }

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

  303.         /** DUT at first entry. */
  304.         private double firstDUT;

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

  309.         /** Interpolator to use. */
  310.         private final HermiteInterpolator interpolator;

  312.         /** Interpolation date. */
  313.         private AbsoluteDate date;

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

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

  351.         /** Get the interpolated value.
  352.          * @return interpolated value
  353.          */
  354.         public double getInterpolated() {
  355.             final double interpolated = interpolator.value(0)[0];
  356.             return beforeLeap ? interpolated : interpolated + 1.0;
  357.         }

  359.     }

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

  364.         /** DUT at first entry. */
  365.         private double firstDUT;

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

  370.         /** Interpolator to use. */
  371.         private final FieldHermiteInterpolator<T> interpolator;

  373.         /** Interpolation date. */
  374.         private FieldAbsoluteDate<T> date;

  376.         /** Interpolation date. */
  377.         private AbsoluteDate absDate;

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

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

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

  421.     }

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

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

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

  451.         // we have EOP data for date -> interpolate correction
  452.         double interpolated = interpolate(date, entry -> entry.getLOD());
  453.         if (tidalCorrection != null) {
  454.             interpolated += tidalCorrection.value(date)[3];
  455.         }
  456.         return interpolated;

  458.     }

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

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

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

  477.         // we have EOP data for date -> interpolate correction
  478.         T interpolated = interpolate(date, aDate, entry -> entry.getLOD());
  479.         if (tidalCorrection != null) {
  480.             interpolated = interpolated.add(tidalCorrection.value(date)[3]);
  481.         }

  483.         return interpolated;

  485.     }

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

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

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

  517.     }

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

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

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

  541.         // we have EOP data for date -> interpolate correction
  542.         final T[] interpolated = interpolate(date, aDate, entry -> entry.getX(), entry -> entry.getY());
  543.         if (tidalCorrection != null) {
  544.             final T[] correction = tidalCorrection.value(date);
  545.             interpolated[0] = interpolated[0].add(correction[0]);
  546.             interpolated[1] = interpolated[1].add(correction[1]);
  547.         }
  548.         return new FieldPoleCorrection<>(interpolated[0], interpolated[1]);

  550.     }

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

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

  566.         // we have EOP data for date -> interpolate correction
  567.         return interpolate(date, entry -> entry.getDdPsi(), entry -> entry.getDdEps());

  569.     }

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

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

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

  588.         // we have EOP data for date -> interpolate correction
  589.         return interpolate(date, aDate, entry -> entry.getDdPsi(), entry -> entry.getDdEps());

  591.     }

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

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

  607.         // we have EOP data for date -> interpolate correction
  608.         return interpolate(date, entry -> entry.getDx(), entry -> entry.getDy());

  610.     }

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

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

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

  629.         // we have EOP data for date -> interpolate correction
  630.         return interpolate(date, aDate, entry -> entry.getDx(), entry -> entry.getDy());

  632.     }

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

  641.         // check if there is data for date
  642.         if (!this.hasDataFor(date)) {
  643.             // no EOP data available for this date, we use a default ITRF 2014
  644.             return ITRFVersion.ITRF_2014;
  645.         }

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

  652.         } catch (TimeStampedCacheException tce) {
  653.             // this should not happen because of date check performed at start
  654.             throw new OrekitInternalError(tce);
  655.         }

  657.     }

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

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

  673.             // prepare next iteration
  674.             preceding = current;

  676.         }
  677.     }

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

  696.     /** Get a non-modifiable view of the EOP entries.
  697.      * @return non-modifiable view of the EOP entries
  698.      */
  699.     public List<EOPEntry> getEntries() {
  700.         return cache.getAll();
  701.     }

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

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

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

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

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

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

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

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

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

  912.     /** Replace the instance with a data transfer object for serialization.
  913.      * @return data transfer object that will be serialized
  914.      */
  915.     @DefaultDataContext
  916.     private Object writeReplace() {
  917.         return new DataTransferObject(conventions, interpolationDegree, getEntries(), tidalCorrection == null);
  918.     }

  920.     /** Internal class used only for serialization. */
  921.     @DefaultDataContext
  922.     private static class DataTransferObject implements Serializable {

  924.         /** Serializable UID. */
  925.         private static final long serialVersionUID = 20231122L;

  927.         /** IERS conventions. */
  928.         private final IERSConventions conventions;

  930.         /** Interpolation degree.
  931.          * @since 12.0
  932.          */
  933.         private final int interpolationDegree;

  935.         /** EOP entries. */
  936.         private final List<EOPEntry> entries;

  938.         /** Indicator for simple interpolation without tidal effects. */
  939.         private final boolean simpleEOP;

  941.         /** Simple constructor.
  942.          * @param conventions IERS conventions to which EOP refers
  943.          * @param interpolationDegree interpolation degree (must be of the form 4k-1)
  944.          * @param entries the EOP data to use
  945.          * @param simpleEOP if true, tidal effects are ignored when interpolating EOP
  946.          * @since 12.0
  947.          */
  948.         DataTransferObject(final IERSConventions conventions,
  949.                            final int interpolationDegree,
  950.                            final List<EOPEntry> entries,
  951.                            final boolean simpleEOP) {
  952.             this.conventions         = conventions;
  953.             this.interpolationDegree = interpolationDegree;
  954.             this.entries             = entries;
  955.             this.simpleEOP           = simpleEOP;
  956.         }

  958.         /** Replace the deserialized data transfer object with a {@link EOPHistory}.
  959.          * @return replacement {@link EOPHistory}
  960.          */
  961.         private Object readResolve() {
  962.             try {
  963.                 return new EOPHistory(conventions, interpolationDegree, entries, simpleEOP);
  964.             } catch (OrekitException oe) {
  965.                 throw new OrekitInternalError(oe);
  966.             }
  967.         }

  969.     }

  971.     /** Internal class for caching tidal correction. */
  972.     private static class TidalCorrectionEntry implements TimeStamped {

  974.         /** Entry date. */
  975.         private final AbsoluteDate date;

  977.         /** Correction. */
  978.         private final double[] correction;

  980.         /** Simple constructor.
  981.          * @param date entry date
  982.          * @param correction correction on the EOP parameters (xp, yp, ut1, lod)
  983.          */
  984.         TidalCorrectionEntry(final AbsoluteDate date, final double[] correction) {
  985.             this.date       = date;
  986.             this.correction = correction;
  987.         }

  989.         /** {@inheritDoc} */
  990.         @Override
  991.         public AbsoluteDate getDate() {
  992.             return date;
  993.         }

  995.     }

  997.     /** Local generator for thread-safe cache. */
  998.     private static class CachedCorrection
  999.         implements TimeVectorFunction, TimeStampedGenerator<TidalCorrectionEntry> {

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

  1004.         /** Step between generated entries. */
  1005.         private final double step;

  1007.         /** Tidal corrections entries cache. */
  1008.         private final TimeStampedCache<TidalCorrectionEntry> cache;

  1010.         /** Simple constructor.
  1011.          * @param tidalCorrection function computing the tidal correction
  1012.          */
  1013.         CachedCorrection(final TimeVectorFunction tidalCorrection) {
  1014.             this.step            = 60 * 60;
  1015.             this.tidalCorrection = tidalCorrection;
  1016.             this.cache           =
  1017.                 new GenericTimeStampedCache<TidalCorrectionEntry>(8,
  1018.                                                                   OrekitConfiguration.getCacheSlotsNumber(),
  1019.                                                                   Constants.JULIAN_DAY * 30,
  1020.                                                                   Constants.JULIAN_DAY,
  1021.                                                                   this);
  1022.         }

  1024.         /** {@inheritDoc} */
  1025.         @Override
  1026.         public double[] value(final AbsoluteDate date) {
  1027.             try {
  1028.                 // set up an interpolator
  1029.                 final HermiteInterpolator interpolator = new HermiteInterpolator();
  1030.                 cache.getNeighbors(date).forEach(entry -> interpolator.addSamplePoint(entry.date.durationFrom(date), entry.correction));

  1032.                 // interpolate to specified date
  1033.                 return interpolator.value(0.0);
  1034.             } catch (TimeStampedCacheException tsce) {
  1035.                 // this should never happen
  1036.                 throw new OrekitInternalError(tsce);
  1037.             }
  1038.         }

  1040.         /** {@inheritDoc} */
  1041.         @Override
  1042.         public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
  1043.             try {

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

  1047.                 final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
  1048.                 final T[] y = MathArrays.buildArray(date.getField(), 4);
  1049.                 final T zero = date.getField().getZero();
  1050.                 final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
  1051.                                                                                            // for example removing derivatives
  1052.                                                                                            // if T was DerivativeStructure
  1053.                 cache.getNeighbors(aDate).forEach(entry -> {
  1054.                     for (int i = 0; i < y.length; ++i) {
  1055.                         y[i] = zero.newInstance(entry.correction[i]);
  1056.                     }
  1057.                     interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
  1058.                 });

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

  1063.             } catch (TimeStampedCacheException tsce) {
  1064.                 // this should never happen
  1065.                 throw new OrekitInternalError(tsce);
  1066.             }
  1067.         }

  1069.         /** {@inheritDoc} */
  1070.         @Override
  1071.         public List<TidalCorrectionEntry> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {

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

  1075.             if (existingDate == null) {

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

  1083.             } else {

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

  1088.                 AbsoluteDate t = existingDate;
  1089.                 if (date.compareTo(t) > 0) {
  1090.                     // forward generation
  1091.                     do {
  1092.                         t = t.shiftedBy(step);
  1093.                         generated.add(new TidalCorrectionEntry(t, tidalCorrection.value(t)));
  1094.                     } while (t.compareTo(date) <= 0);
  1095.                 } else {
  1096.                     // backward generation
  1097.                     do {
  1098.                         t = t.shiftedBy(-step);
  1099.                         generated.add(0, new TidalCorrectionEntry(t, tidalCorrection.value(t)));
  1100.                     } while (t.compareTo(date) >= 0);
  1101.                 }
  1102.             }

  1104.             // return the generated transforms
  1105.             return generated;

  1107.         }
  1108.     }

  1110. }
Created with JaCoCo 0.8.12.202403310830