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 /** Get the EOP data type.
656 * @param date date at which the value is desired
657 * @return data type of the EOP covering the specified date
658 * @since 13.1.1
659 */
660 public EopDataType getEopDataType(final AbsoluteDate date) {
661
662 // check if there is data for date
663 if (!this.hasDataFor(date)) {
664 // no EOP data available for this date, data type is unknown
665 return EopDataType.UNKNOWN;
666 }
667
668 // we have EOP data for date
669 final Optional<EOPEntry> first = getNeighbors(date, 1).findFirst();
670 return first.isPresent() ? first.get().getEopDataType() : EopDataType.UNKNOWN;
671
672 }
673
674 /** Check Earth orientation parameters continuity.
675 * @param maxGap maximal allowed gap between entries (in seconds)
676 */
677 public void checkEOPContinuity(final double maxGap) {
678 TimeStamped preceding = null;
679 for (final TimeStamped current : this.cache.getAll()) {
680
681 // compare the dates of preceding and current entries
682 if (preceding != null && (current.getDate().durationFrom(preceding.getDate())) > maxGap) {
683 throw new OrekitException(OrekitMessages.MISSING_EARTH_ORIENTATION_PARAMETERS_BETWEEN_DATES_GAP,
684 preceding.getDate(), current.getDate(),
685 current.getDate().durationFrom(preceding.getDate()));
686 }
687
688 // prepare next iteration
689 preceding = current;
690
691 }
692 }
693
694 /**
695 * Check if the cache has data for the given date using
696 * {@link #getStartDate()} and {@link #getEndDate()}.
697 *
698 * @param date the requested date
699 * @return true if the {@link #cache} has data for the requested date, false
700 * otherwise.
701 */
702 protected boolean hasDataFor(final AbsoluteDate date) {
703 /*
704 * when there is no EOP data, short circuit getStartDate, which will
705 * throw an exception
706 */
707 return this.hasData && this.getStartDate().compareTo(date) <= 0 &&
708 date.compareTo(this.getEndDate()) <= 0;
709 }
710
711 /** Get a non-modifiable view of the EOP entries.
712 * @return non-modifiable view of the EOP entries
713 */
714 public List<EOPEntry> getEntries() {
715 return cache.getAll();
716 }
717
718 /** Interpolate a single EOP component.
719 * <p>
720 * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
721 * </p>
722 * @param date interpolation date
723 * @param selector selector for EOP entry component
724 * @return interpolated value
725 */
726 private double interpolate(final AbsoluteDate date, final Function<EOPEntry, Double> selector) {
727 try {
728 final HermiteInterpolator interpolator = new HermiteInterpolator();
729 getNeighbors(date, interpolationDegree + 1).
730 forEach(entry -> interpolator.addSamplePoint(entry.getDate().durationFrom(date),
731 new double[] {
732 selector.apply(entry)
733 }));
734 return interpolator.value(0)[0];
735 } catch (TimeStampedCacheException tce) {
736 // this should not happen because of date check performed by caller
737 throw new OrekitInternalError(tce);
738 }
739 }
740
741 /** Interpolate a single EOP component.
742 * <p>
743 * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
744 * </p>
745 * @param date interpolation date
746 * @param aDate interpolation date, as an {@link AbsoluteDate}
747 * @param selector selector for EOP entry component
748 * @param <T> type of the field elements
749 * @return interpolated value
750 */
751 private <T extends CalculusFieldElement<T>> T interpolate(final FieldAbsoluteDate<T> date,
752 final AbsoluteDate aDate,
753 final Function<EOPEntry, Double> selector) {
754 try {
755 final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
756 final T[] y = MathArrays.buildArray(date.getField(), 1);
757 final T zero = date.getField().getZero();
758 final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
759 // for example removing derivatives
760 // if T was DerivativeStructure
761 getNeighbors(aDate, interpolationDegree + 1).
762 forEach(entry -> {
763 y[0] = zero.newInstance(selector.apply(entry));
764 interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
765 });
766 return interpolator.value(date.durationFrom(central))[0]; // here, we introduce derivatives again (in DerivativeStructure case)
767 } catch (TimeStampedCacheException tce) {
768 // this should not happen because of date check performed by caller
769 throw new OrekitInternalError(tce);
770 }
771 }
772
773 /** Interpolate two EOP components.
774 * <p>
775 * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
776 * </p>
777 * @param date interpolation date
778 * @param selector1 selector for first EOP entry component
779 * @param selector2 selector for second EOP entry component
780 * @return interpolated value
781 */
782 private double[] interpolate(final AbsoluteDate date,
783 final Function<EOPEntry, Double> selector1,
784 final Function<EOPEntry, Double> selector2) {
785 try {
786 final HermiteInterpolator interpolator = new HermiteInterpolator();
787 getNeighbors(date, interpolationDegree + 1).
788 forEach(entry -> interpolator.addSamplePoint(entry.getDate().durationFrom(date),
789 new double[] {
790 selector1.apply(entry),
791 selector2.apply(entry)
792 }));
793 return interpolator.value(0);
794 } catch (TimeStampedCacheException tce) {
795 // this should not happen because of date check performed by caller
796 throw new OrekitInternalError(tce);
797 }
798 }
799
800 /** Interpolate two EOP components.
801 * <p>
802 * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
803 * </p>
804 * @param date interpolation date
805 * @param selector1 selector for first EOP entry component
806 * @param selector2 selector for second EOP entry component
807 * @param selector1Rate selector for first EOP entry component rate
808 * @param selector2Rate selector for second EOP entry component rate
809 * @return interpolated value
810 * @since 12.0
811 */
812 private double[] interpolate(final AbsoluteDate date,
813 final Function<EOPEntry, Double> selector1,
814 final Function<EOPEntry, Double> selector2,
815 final Function<EOPEntry, Double> selector1Rate,
816 final Function<EOPEntry, Double> selector2Rate) {
817 try {
818 final HermiteInterpolator interpolator = new HermiteInterpolator();
819 getNeighbors(date, (interpolationDegree + 1) / 2).
820 forEach(entry -> interpolator.addSamplePoint(entry.getDate().durationFrom(date),
821 new double[] {
822 selector1.apply(entry),
823 selector2.apply(entry)
824 },
825 new double[] {
826 selector1Rate.apply(entry),
827 selector2Rate.apply(entry)
828 }));
829 return interpolator.value(0);
830 } catch (TimeStampedCacheException tce) {
831 // this should not happen because of date check performed by caller
832 throw new OrekitInternalError(tce);
833 }
834 }
835
836 /** Interpolate two EOP components.
837 * <p>
838 * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
839 * </p>
840 * @param date interpolation date
841 * @param aDate interpolation date, as an {@link AbsoluteDate}
842 * @param selector1 selector for first EOP entry component
843 * @param selector2 selector for second EOP entry component
844 * @param <T> type of the field elements
845 * @return interpolated value
846 */
847 private <T extends CalculusFieldElement<T>> T[] interpolate(final FieldAbsoluteDate<T> date,
848 final AbsoluteDate aDate,
849 final Function<EOPEntry, Double> selector1,
850 final Function<EOPEntry, Double> selector2) {
851 try {
852 final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
853 final T[] y = MathArrays.buildArray(date.getField(), 2);
854 final T zero = date.getField().getZero();
855 final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
856 // for example removing derivatives
857 // if T was DerivativeStructure
858 getNeighbors(aDate, interpolationDegree + 1).
859 forEach(entry -> {
860 y[0] = zero.newInstance(selector1.apply(entry));
861 y[1] = zero.newInstance(selector2.apply(entry));
862 interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
863 });
864 return interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
865 } catch (TimeStampedCacheException tce) {
866 // this should not happen because of date check performed by caller
867 throw new OrekitInternalError(tce);
868 }
869 }
870
871 /** Check if some derivatives are missing.
872 * @param raw raw EOP history
873 * @return true if some derivatives are missing
874 * @since 12.0
875 */
876 private boolean missSomeDerivatives(final Collection<? extends EOPEntry> raw) {
877 for (final EOPEntry entry : raw) {
878 if (Double.isNaN(entry.getLOD() + entry.getXRate() + entry.getYRate())) {
879 return true;
880 }
881 }
882 return false;
883 }
884
885 /** Fix missing derivatives.
886 * @param entry EOP entry to fix
887 * @param rawCache raw EOP history cache
888 * @return fixed entry
889 * @since 12.0
890 */
891 private EOPEntry fixDerivatives(final EOPEntry entry, final ImmutableTimeStampedCache<EOPEntry> rawCache) {
892
893 // helper function to differentiate some EOP parameters
894 final BiFunction<EOPEntry, Function<EOPEntry, Double>, Double> differentiator =
895 (e, selector) -> {
896 final HermiteInterpolator interpolator = new HermiteInterpolator();
897 rawCache.getNeighbors(e.getDate()).
898 forEach(f -> interpolator.addSamplePoint(f.getDate().durationFrom(e.getDate()),
899 new double[] {
900 selector.apply(f)
901 }));
902 return interpolator.derivatives(0.0, 1)[1][0];
903 };
904
905 if (Double.isNaN(entry.getLOD() + entry.getXRate() + entry.getYRate())) {
906 final double lod = Double.isNaN(entry.getLOD()) ?
907 -differentiator.apply(entry, EOPEntry::getUT1MinusUTC) :
908 entry.getLOD();
909 final double xRate = Double.isNaN(entry.getXRate()) ?
910 differentiator.apply(entry, EOPEntry::getX) :
911 entry.getXRate();
912 final double yRate = Double.isNaN(entry.getYRate()) ?
913 differentiator.apply(entry, EOPEntry::getY) :
914 entry.getYRate();
915 return new EOPEntry(entry.getMjd(),
916 entry.getUT1MinusUTC(), lod,
917 entry.getX(), entry.getY(), xRate, yRate,
918 entry.getDdPsi(), entry.getDdEps(),
919 entry.getDx(), entry.getDy(),
920 entry.getITRFType(), entry.getDate(), entry.getEopDataType());
921 } else {
922 // the entry already has all derivatives
923 return entry;
924 }
925 }
926
927 /** Internal class for caching tidal correction. */
928 private static class TidalCorrectionEntry implements TimeStamped {
929
930 /** Entry date. */
931 private final AbsoluteDate date;
932
933 /** Correction. */
934 private final double[] correction;
935
936 /** Simple constructor.
937 * @param date entry date
938 * @param correction correction on the EOP parameters (xp, yp, ut1, lod)
939 */
940 TidalCorrectionEntry(final AbsoluteDate date, final double[] correction) {
941 this.date = date;
942 this.correction = correction;
943 }
944
945 /** {@inheritDoc} */
946 @Override
947 public AbsoluteDate getDate() {
948 return date;
949 }
950
951 }
952
953 /** Local generator for thread-safe cache. */
954 private static class CachedCorrection
955 implements TimeVectorFunction, TimeStampedGenerator<TidalCorrectionEntry> {
956
957 /** Correction to apply to EOP (may be null). */
958 private final TimeVectorFunction tidalCorrection;
959
960 /** Step between generated entries. */
961 private final double step;
962
963 /** Tidal corrections entries cache. */
964 private final TimeStampedCache<TidalCorrectionEntry> cache;
965
966 /** Simple constructor.
967 * @param tidalCorrection function computing the tidal correction
968 */
969 CachedCorrection(final TimeVectorFunction tidalCorrection) {
970 this.step = 60 * 60;
971 this.tidalCorrection = tidalCorrection;
972 this.cache =
973 new GenericTimeStampedCache<>(8,
974 OrekitConfiguration.getCacheSlotsNumber(),
975 Constants.JULIAN_DAY * 30,
976 Constants.JULIAN_DAY,
977 this);
978 }
979
980 /** {@inheritDoc} */
981 @Override
982 public double[] value(final AbsoluteDate date) {
983 try {
984 // set up an interpolator
985 final HermiteInterpolator interpolator = new HermiteInterpolator();
986 cache.getNeighbors(date).forEach(entry -> interpolator.addSamplePoint(entry.date.durationFrom(date), entry.correction));
987
988 // interpolate to specified date
989 return interpolator.value(0.0);
990 } catch (TimeStampedCacheException tsce) {
991 // this should never happen
992 throw new OrekitInternalError(tsce);
993 }
994 }
995
996 /** {@inheritDoc} */
997 @Override
998 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
999 try {
1000
1001 final AbsoluteDate aDate = date.toAbsoluteDate();
1002
1003 final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
1004 final T[] y = MathArrays.buildArray(date.getField(), 4);
1005 final T zero = date.getField().getZero();
1006 final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
1007 // for example removing derivatives
1008 // if T was DerivativeStructure
1009 cache.getNeighbors(aDate).forEach(entry -> {
1010 for (int i = 0; i < y.length; ++i) {
1011 y[i] = zero.newInstance(entry.correction[i]);
1012 }
1013 interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
1014 });
1015
1016 // interpolate to specified date
1017 return interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
1018
1019 } catch (TimeStampedCacheException tsce) {
1020 // this should never happen
1021 throw new OrekitInternalError(tsce);
1022 }
1023 }
1024
1025 /** {@inheritDoc} */
1026 @Override
1027 public List<TidalCorrectionEntry> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {
1028
1029 final List<TidalCorrectionEntry> generated = new ArrayList<>();
1030
1031 if (existingDate == null) {
1032
1033 // no prior existing entries, just generate a first set
1034 for (int i = -cache.getMaxNeighborsSize() / 2; generated.size() < cache.getMaxNeighborsSize(); ++i) {
1035 final AbsoluteDate t = date.shiftedBy(i * step);
1036 generated.add(new TidalCorrectionEntry(t, tidalCorrection.value(t)));
1037 }
1038
1039 } else {
1040
1041 // some entries have already been generated
1042 // add the missing ones up to specified date
1043
1044 AbsoluteDate t = existingDate;
1045 if (date.compareTo(t) > 0) {
1046 // forward generation
1047 do {
1048 t = t.shiftedBy(step);
1049 generated.add(new TidalCorrectionEntry(t, tidalCorrection.value(t)));
1050 } while (t.compareTo(date) <= 0);
1051 } else {
1052 // backward generation
1053 do {
1054 t = t.shiftedBy(-step);
1055 generated.add(0, new TidalCorrectionEntry(t, tidalCorrection.value(t)));
1056 } while (t.compareTo(date) >= 0);
1057 }
1058 }
1059
1060 // return the generated transforms
1061 return generated;
1062
1063 }
1064 }
1065
1066 }