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;
18  
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;
27  
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;
51  
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 {
57  
58      /** Default interpolation degree.
59       * @since 12.0
60       */
61      public static final int DEFAULT_INTERPOLATION_DEGREE = 3;
62  
63      /** Interpolation degree.
64       * @since 12.0
65       */
66      private final int interpolationDegree;
67  
68      /**
69       * If this history has any EOP data.
70       *
71       * @see #hasDataFor(AbsoluteDate)
72       */
73      private final boolean hasData;
74  
75      /** EOP history entries. */
76      private final transient ImmutableTimeStampedCache<EOPEntry> cache;
77  
78      /** IERS conventions to which EOP refers. */
79      private final IERSConventions conventions;
80  
81      /** Correction to apply to EOP (may be null). */
82      private final transient TimeVectorFunction tidalCorrection;
83  
84      /** Time scales to use when computing corrections. */
85      private final transient TimeScales timeScales;
86  
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     }
104 
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     }
122 
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) {
136 
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         }
142 
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     }
168 
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     }
177 
178     /** Get interpolation degree.
179      * @return interpolation degree
180      * @since 12.0
181      */
182     public int getInterpolationDegree() {
183         return interpolationDegree;
184     }
185 
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     }
195 
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     }
205 
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     }
213 
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     }
220 
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     }
227 
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     }
234 
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) {
241 
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         }
247 
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         }
261 
262     }
263 
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) {
272 
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         }
279 
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         }
293 
294     }
295 
296     /** Local class for DUT1 interpolation, crossing leaps safely. */
297     private static class DUT1Interpolator implements Consumer<EOPEntry> {
298 
299         /** DUT at first entry. */
300         private double firstDUT;
301 
302         /** Indicator for dates just before a leap occurring during the interpolation sample. */
303         private boolean beforeLeap;
304 
305         /** Interpolator to use. */
306         private final HermiteInterpolator interpolator;
307 
308         /** Interpolation date. */
309         private AbsoluteDate date;
310 
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         }
320 
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         }
346 
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         }
354 
355     }
356 
357     /** Local class for DUT1 interpolation, crossing leaps safely. */
358     private static class FieldDUT1Interpolator<T extends CalculusFieldElement<T>> implements Consumer<EOPEntry> {
359 
360         /** DUT at first entry. */
361         private double firstDUT;
362 
363         /** Indicator for dates just before a leap occurring during the interpolation sample. */
364         private boolean beforeLeap;
365 
366         /** Interpolator to use. */
367         private final FieldHermiteInterpolator<T> interpolator;
368 
369         /** Interpolation date. */
370         private FieldAbsoluteDate<T> date;
371 
372         /** Interpolation date. */
373         private AbsoluteDate absDate;
374 
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         }
386 
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         }
408 
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         }
416 
417     }
418 
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     }
433 
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) {
440 
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         }
446 
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;
453 
454     }
455 
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) {
464 
465         final AbsoluteDate aDate = date.toAbsoluteDate();
466 
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         }
472 
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         }
478 
479         return interpolated;
480 
481     }
482 
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) {
490 
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         }
501 
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]);
512 
513     }
514 
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) {
523 
524         final AbsoluteDate aDate = date.toAbsoluteDate();
525 
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         }
536 
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]);
545 
546     }
547 
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) {
555 
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         }
561 
562         // we have EOP data for date -> interpolate correction
563         return interpolate(date, EOPEntry::getDdPsi, EOPEntry::getDdEps);
564 
565     }
566 
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) {
575 
576         final AbsoluteDate aDate = date.toAbsoluteDate();
577 
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         }
583 
584         // we have EOP data for date -> interpolate correction
585         return interpolate(date, aDate, EOPEntry::getDdPsi, EOPEntry::getDdEps);
586 
587     }
588 
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) {
596 
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         }
602 
603         // we have EOP data for date -> interpolate correction
604         return interpolate(date, EOPEntry::getDx, EOPEntry::getDy);
605 
606     }
607 
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) {
616 
617         final AbsoluteDate aDate = date.toAbsoluteDate();
618 
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         }
624 
625         // we have EOP data for date -> interpolate correction
626         return interpolate(date, aDate, EOPEntry::getDx, EOPEntry::getDy);
627 
628     }
629 
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) {
636 
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         }
642 
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;
647 
648         } catch (TimeStampedCacheException tce) {
649             // this should not happen because of date check performed at start
650             throw new OrekitInternalError(tce);
651         }
652 
653     }
654 
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()) {
661 
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             }
668 
669             // prepare next iteration
670             preceding = current;
671 
672         }
673     }
674 
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     }
691 
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     }
698 
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     }
721 
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     }
753 
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     }
780 
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     }
816 
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     }
851 
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     }
865 
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) {
873 
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                         };
885 
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     }
907 
908     /** Internal class for caching tidal correction. */
909     private static class TidalCorrectionEntry implements TimeStamped {
910 
911         /** Entry date. */
912         private final AbsoluteDate date;
913 
914         /** Correction. */
915         private final double[] correction;
916 
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         }
925 
926         /** {@inheritDoc} */
927         @Override
928         public AbsoluteDate getDate() {
929             return date;
930         }
931 
932     }
933 
934     /** Local generator for thread-safe cache. */
935     private static class CachedCorrection
936         implements TimeVectorFunction, TimeStampedGenerator<TidalCorrectionEntry> {
937 
938         /** Correction to apply to EOP (may be null). */
939         private final TimeVectorFunction tidalCorrection;
940 
941         /** Step between generated entries. */
942         private final double step;
943 
944         /** Tidal corrections entries cache. */
945         private final TimeStampedCache<TidalCorrectionEntry> cache;
946 
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         }
960 
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));
968 
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         }
976 
977         /** {@inheritDoc} */
978         @Override
979         public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
980             try {
981 
982                 final AbsoluteDate aDate = date.toAbsoluteDate();
983 
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                 });
996 
997                 // interpolate to specified date
998                 return interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
999 
1000             } catch (TimeStampedCacheException tsce) {
1001                 // this should never happen
1002                 throw new OrekitInternalError(tsce);
1003             }
1004         }
1005 
1006         /** {@inheritDoc} */
1007         @Override
1008         public List<TidalCorrectionEntry> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {
1009 
1010             final List<TidalCorrectionEntry> generated = new ArrayList<>();
1011 
1012             if (existingDate == null) {
1013 
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                 }
1019 
1020             } else {
1021 
1022                 // some entries have already been generated
1023                 // add the missing ones up to specified date
1024 
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             }
1040 
1041             // return the generated transforms
1042             return generated;
1043 
1044         }
1045     }
1046 
1047 }