1 /* Copyright 2002-2026 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.ChronologicalComparator;
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;
52
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 {
58
59 /** Default interpolation degree.
60 * @since 12.0
61 */
62 public static final int DEFAULT_INTERPOLATION_DEGREE = 3;
63
64 /** Interpolation degree.
65 * @since 12.0
66 */
67 private final int interpolationDegree;
68
69 /**
70 * If this history has any EOP data.
71 *
72 * @see #hasDataFor(AbsoluteDate)
73 */
74 private final boolean hasData;
75
76 /** EOP history entries. */
77 private final ImmutableTimeStampedCache<EOPEntry> cache;
78
79 /** IERS conventions to which EOP refers. */
80 private final IERSConventions conventions;
81
82 /** Correction to apply to EOP (may be null). */
83 private final TimeVectorFunction tidalCorrection;
84
85 /** Time scales to use when computing corrections. */
86 private final TimeScales timeScales;
87
88 /** Simple constructor.
89 * <p>
90 * This method uses the {@link DataContext#getDefault() default data context}.
91 * </p>
92 * <p>
93 * Since 14.0, the {@code data} provided to this constructor may be scattered,
94 * unsorted and incomplete as it may come from different files. This is typically
95 * true with Bulletin A, where rapid data pole motion is published weekly and
96 * nutation data is published monthly, several weeks later. The data is therefore
97 * sorted, deduplicated and merged internally. For this reason, it is <em>not</em>
98 * recommended to use {@link java.util.SortedSet} for providing the data, as this
99 * would preserve only one entry for each date and eliminate arbitrarily the other
100 * entries for the same date.
101 * </p>
102 * @param conventions IERS conventions to which EOP refers
103 * @param interpolationDegree interpolation degree (must be of the form 4k-1)
104 * @param data the EOP data to use
105 * (it is <em>not</em> recommended to use {@link java.util.SortedSet} here, see above)
106 * @param simpleEOP if true, tidal effects are ignored when interpolating EOP
107 * @see #EOPHistory(IERSConventions, int, Collection, boolean, TimeScales)
108 */
109 @DefaultDataContext
110 protected EOPHistory(final IERSConventions conventions,
111 final int interpolationDegree,
112 final Collection<? extends EOPEntry> data,
113 final boolean simpleEOP) {
114 this(conventions, interpolationDegree, data, simpleEOP, DataContext.getDefault().getTimeScales());
115 }
116
117 /** Simple constructor.
118 * <p>
119 * Since 14.0, the {@code data} provided to this constructor may be scattered,
120 * unsorted and incomplete as it may come from different files. This is typically
121 * true with Bulletin A, where rapid data pole motion is published weekly and
122 * nutation data is published monthly, several weeks later. The data is therefore
123 * sorted, deduplicated and merged internally. For this reason, it is <em>not</em>
124 * recommended to use {@link java.util.SortedSet} for providing the data, as this
125 * would preserve only one entry for each date and eliminate arbitrarily the other
126 * entries for the same date.
127 * </p>
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 * (it is <em>not</em> recommended to use {@link java.util.SortedSet} here, see above)
132 * @param simpleEOP if true, tidal effects are ignored when interpolating EOP
133 * @param timeScales to use when computing EOP corrections.
134 * @since 10.1
135 */
136 public EOPHistory(final IERSConventions conventions,
137 final int interpolationDegree,
138 final Collection<? extends EOPEntry> data,
139 final boolean simpleEOP,
140 final TimeScales timeScales) {
141 this(conventions, interpolationDegree, data,
142 simpleEOP ? null : new CachedCorrection(conventions.getEOPTidalCorrection(timeScales)),
143 timeScales);
144 }
145
146 /** Simple constructor.
147 * @param conventions IERS conventions to which EOP refers
148 * @param interpolationDegree interpolation degree (must be of the form 4k-1)
149 * @param data the EOP data to use
150 * @param tidalCorrection correction to apply to EOP
151 * @param timeScales to use when computing EOP corrections
152 * @since 12
153 */
154 private EOPHistory(final IERSConventions conventions,
155 final int interpolationDegree,
156 final Collection<? extends EOPEntry> data,
157 final TimeVectorFunction tidalCorrection,
158 final TimeScales timeScales) {
159
160 // check interpolation degree is 4k-1
161 final int k = (interpolationDegree + 1) / 4;
162 if (interpolationDegree != 4 * k - 1) {
163 throw new OrekitException(OrekitMessages.WRONG_EOP_INTERPOLATION_DEGREE, interpolationDegree);
164 }
165
166 // deduplicate data
167 final List<EOPEntry> deduplicated = deduplicate(data);
168
169 this.conventions = conventions;
170 this.interpolationDegree = interpolationDegree;
171 this.tidalCorrection = tidalCorrection;
172 this.timeScales = timeScales;
173 if (!deduplicated.isEmpty()) {
174 // enough data to interpolate
175 if (missSomeDerivatives(deduplicated)) {
176 // we need to estimate the missing derivatives
177 final ImmutableTimeStampedCache<EOPEntry> rawCache =
178 new ImmutableTimeStampedCache<>(FastMath.min(interpolationDegree + 1,
179 deduplicated.size()), deduplicated);
180 final List<EOPEntry> fixedData = new ArrayList<>();
181 for (final EOPEntry entry : rawCache.getAll()) {
182 fixedData.add(fixDerivatives(entry, rawCache));
183 }
184 cache = new ImmutableTimeStampedCache<>(FastMath.min(interpolationDegree + 1,
185 fixedData.size()), fixedData);
186 } else {
187 cache = new ImmutableTimeStampedCache<>(FastMath.min(interpolationDegree + 1,
188 deduplicated.size()), deduplicated);
189 }
190 hasData = true;
191 } else {
192 // not enough data to interpolate -> always use null correction
193 cache = ImmutableTimeStampedCache.emptyCache();
194 hasData = false;
195 }
196 }
197
198 /** Deduplicate entries.
199 * <p>
200 * Some sources provide EOP data with scattered and incomplete entries. A typical
201 * example is Bulletin A, as the xp, yp and UT1-UTC data are published as rapid data
202 * on a weekly basis whereas pole offsets Δδψ/Δδε and x/y are published on a monthly
203 * basis. This implies data for one day is scattered among several files published
204 * weeks apart and need to be merged. This method performs this merge, combining
205 * entries that correspond to a single date. The result is a sorted list with only
206 * one entry for each date.
207 * </p>
208 * @param rawData raw data, that may contain several partial entries for some dates
209 * @return deduplicated data, with only one entry for each date, with merged fields
210 * @since 14.0
211 */
212 private List<EOPEntry> deduplicate(final Collection<? extends EOPEntry> rawData) {
213
214 // copy data into an independent list
215 final List<EOPEntry> deduplicated = new ArrayList<>(rawData);
216
217 // sort chronologically
218 deduplicated.sort(new ChronologicalComparator());
219
220 // process each entry, combining the ones that correspond to the same date
221 int i = 0;
222 while (i < deduplicated.size() - 1) {
223
224 // beware we must call get(i) and get(i + 1) at each iteration because
225 // we mess up with the both the list and the index within the loop
226 final EOPEntry current = deduplicated.get(i);
227 final EOPEntry next = deduplicated.get(i + 1);
228
229 if (next.getDate().equals(current.getDate())) {
230 // the two entries are at the same date, we combine them together
231 deduplicated.set(i, new EOPEntry(current, next));
232 deduplicated.remove(i + 1);
233 } else {
234 // this entry has been completed, we can go to the next one
235 ++i;
236 }
237
238 }
239
240 return deduplicated;
241
242 }
243
244 /**
245 * Determine if this history uses simplified EOP corrections.
246 *
247 * @return {@code true} if tidal corrections are ignored, {@code false} otherwise.
248 */
249 public boolean isSimpleEop() {
250 return tidalCorrection == null;
251 }
252
253 /** Get interpolation degree.
254 * @return interpolation degree
255 * @since 12.0
256 */
257 public int getInterpolationDegree() {
258 return interpolationDegree;
259 }
260
261 /**
262 * Get the time scales used in computing EOP corrections.
263 *
264 * @return set of time scales.
265 * @since 10.1
266 */
267 public TimeScales getTimeScales() {
268 return timeScales;
269 }
270
271 /** Get version of the instance that does not cache tidal correction.
272 * @return version of the instance that does not cache tidal correction
273 * @since 12.0
274 */
275 public EOPHistory getEOPHistoryWithoutCachedTidalCorrection() {
276 return new EOPHistory(conventions, interpolationDegree, getEntries(),
277 conventions.getEOPTidalCorrection(timeScales),
278 timeScales);
279 }
280
281 /** Check if the instance caches tidal corrections.
282 * @return true if the instance caches tidal corrections
283 * @since 12.0
284 */
285 public boolean cachesTidalCorrection() {
286 return tidalCorrection instanceof CachedCorrection;
287 }
288
289 /** Get the IERS conventions to which these EOP apply.
290 * @return IERS conventions to which these EOP apply
291 */
292 public IERSConventions getConventions() {
293 return conventions;
294 }
295
296 /** Get the date of the first available Earth Orientation Parameters.
297 * @return the start date of the available data
298 */
299 public AbsoluteDate getStartDate() {
300 return this.cache.getEarliest().getDate();
301 }
302
303 /** Get the date of the last available Earth Orientation Parameters.
304 * @return the end date of the available data
305 */
306 public AbsoluteDate getEndDate() {
307 return this.cache.getLatest().getDate();
308 }
309
310 /** Get the UT1-UTC value.
311 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
312 * @param date date at which the value is desired
313 * @return UT1-UTC in seconds (0 if date is outside covered range)
314 */
315 public double getUT1MinusUTC(final AbsoluteDate date) {
316
317 //check if there is data for date
318 if (!this.hasDataFor(date)) {
319 // no EOP data available for this date, we use a default 0.0 offset
320 return (tidalCorrection == null) ? 0.0 : tidalCorrection.value(date)[2];
321 }
322
323 // we have EOP data -> interpolate offset
324 try {
325 final DUT1Interpolator interpolator = new DUT1Interpolator(date);
326 getNeighbors(date, interpolationDegree + 1).forEach(interpolator);
327 double interpolated = interpolator.getInterpolated();
328 if (tidalCorrection != null) {
329 interpolated += tidalCorrection.value(date)[2];
330 }
331 return interpolated;
332 } catch (TimeStampedCacheException tce) {
333 //this should not happen because of date check above
334 throw new OrekitInternalError(tce);
335 }
336
337 }
338
339 /** Get the UT1-UTC value.
340 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
341 * @param date date at which the value is desired
342 * @param <T> type of the field elements
343 * @return UT1-UTC in seconds (0 if date is outside covered range)
344 * @since 9.0
345 */
346 public <T extends CalculusFieldElement<T>> T getUT1MinusUTC(final FieldAbsoluteDate<T> date) {
347
348 //check if there is data for date
349 final AbsoluteDate absDate = date.toAbsoluteDate();
350 if (!this.hasDataFor(absDate)) {
351 // no EOP data available for this date, we use a default 0.0 offset
352 return (tidalCorrection == null) ? date.getField().getZero() : tidalCorrection.value(date)[2];
353 }
354
355 // we have EOP data -> interpolate offset
356 try {
357 final FieldDUT1Interpolator<T> interpolator = new FieldDUT1Interpolator<>(date, absDate);
358 getNeighbors(absDate, interpolationDegree + 1).forEach(interpolator);
359 T interpolated = interpolator.getInterpolated();
360 if (tidalCorrection != null) {
361 interpolated = interpolated.add(tidalCorrection.value(date)[2]);
362 }
363 return interpolated;
364 } catch (TimeStampedCacheException tce) {
365 //this should not happen because of date check above
366 throw new OrekitInternalError(tce);
367 }
368
369 }
370
371 /** Local class for DUT1 interpolation, crossing leaps safely. */
372 private static class DUT1Interpolator implements Consumer<EOPEntry> {
373
374 /** DUT at first entry. */
375 private double firstDUT;
376
377 /** Indicator for dates just before a leap occurring during the interpolation sample. */
378 private boolean beforeLeap;
379
380 /** Interpolator to use. */
381 private final HermiteInterpolator interpolator;
382
383 /** Interpolation date. */
384 private final AbsoluteDate date;
385
386 /** Simple constructor.
387 * @param date interpolation date
388 */
389 DUT1Interpolator(final AbsoluteDate date) {
390 this.firstDUT = Double.NaN;
391 this.beforeLeap = true;
392 this.interpolator = new HermiteInterpolator();
393 this.date = date;
394 }
395
396 /** {@inheritDoc} */
397 @Override
398 public void accept(final EOPEntry neighbor) {
399 if (Double.isNaN(firstDUT)) {
400 firstDUT = neighbor.getUT1MinusUTC();
401 }
402 final double dut;
403 if (neighbor.getUT1MinusUTC() - firstDUT > 0.9) {
404 // there was a leap second between the entries
405 dut = neighbor.getUT1MinusUTC() - 1.0;
406 // UTCScale considers the discontinuity to occur at the start of the leap
407 // second so this code must use the same convention. EOP entries are time
408 // stamped at midnight UTC so 1 second before is the start of the leap
409 // second.
410 if (neighbor.getDate().shiftedBy(-1).compareTo(date) <= 0) {
411 beforeLeap = false;
412 }
413 } else {
414 dut = neighbor.getUT1MinusUTC();
415 }
416 interpolator.addSamplePoint(neighbor.getDate().durationFrom(date),
417 new double[] {
418 dut
419 });
420 }
421
422 /** Get the interpolated value.
423 * @return interpolated value
424 */
425 public double getInterpolated() {
426 final double interpolated = interpolator.value(0)[0];
427 return beforeLeap ? interpolated : interpolated + 1.0;
428 }
429
430 }
431
432 /** Local class for DUT1 interpolation, crossing leaps safely. */
433 private static class FieldDUT1Interpolator<T extends CalculusFieldElement<T>> implements Consumer<EOPEntry> {
434
435 /** DUT at first entry. */
436 private double firstDUT;
437
438 /** Indicator for dates just before a leap occurring during the interpolation sample. */
439 private boolean beforeLeap;
440
441 /** Interpolator to use. */
442 private final FieldHermiteInterpolator<T> interpolator;
443
444 /** Interpolation date. */
445 private final FieldAbsoluteDate<T> date;
446
447 /** Interpolation date. */
448 private final AbsoluteDate absDate;
449
450 /** Simple constructor.
451 * @param date interpolation date
452 * @param absDate interpolation date
453 */
454 FieldDUT1Interpolator(final FieldAbsoluteDate<T> date, final AbsoluteDate absDate) {
455 this.firstDUT = Double.NaN;
456 this.beforeLeap = true;
457 this.interpolator = new FieldHermiteInterpolator<>();
458 this.date = date;
459 this.absDate = absDate;
460 }
461
462 /** {@inheritDoc} */
463 @Override
464 public void accept(final EOPEntry neighbor) {
465 if (Double.isNaN(firstDUT)) {
466 firstDUT = neighbor.getUT1MinusUTC();
467 }
468 final double dut;
469 if (neighbor.getUT1MinusUTC() - firstDUT > 0.9) {
470 // there was a leap second between the entries
471 dut = neighbor.getUT1MinusUTC() - 1.0;
472 if (neighbor.getDate().compareTo(absDate) <= 0) {
473 beforeLeap = false;
474 }
475 } else {
476 dut = neighbor.getUT1MinusUTC();
477 }
478 final T[] array = MathArrays.buildArray(date.getField(), 1);
479 array[0] = date.getField().getZero().newInstance(dut);
480 interpolator.addSamplePoint(date.durationFrom(neighbor.getDate()).negate(),
481 array);
482 }
483
484 /** Get the interpolated value.
485 * @return interpolated value
486 */
487 public T getInterpolated() {
488 final T interpolated = interpolator.value(date.getField().getZero())[0];
489 return beforeLeap ? interpolated : interpolated.add(1.0);
490 }
491
492 }
493
494 /**
495 * Get the entries surrounding a central date.
496 * <p>
497 * See {@link #hasDataFor(AbsoluteDate)} to determine if the cache has data
498 * for {@code central} without throwing an exception.
499 *
500 * @param central central date
501 * @param n number of neighbors
502 * @return array of cached entries surrounding specified date
503 * @since 12.0
504 */
505 protected Stream<EOPEntry> getNeighbors(final AbsoluteDate central, final int n) {
506 return cache.getNeighbors(central, n);
507 }
508
509 /** Get the LoD (Length of Day) value.
510 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
511 * @param date date at which the value is desired
512 * @return LoD in seconds (0 if date is outside covered range)
513 */
514 public double getLOD(final AbsoluteDate date) {
515
516 // check if there is data for date
517 if (!this.hasDataFor(date)) {
518 // no EOP data available for this date, we use a default null correction
519 return (tidalCorrection == null) ? 0.0 : tidalCorrection.value(date)[3];
520 }
521
522 // we have EOP data for date -> interpolate correction
523 double interpolated = interpolate(date, EOPEntry::getLOD);
524 if (tidalCorrection != null) {
525 interpolated += tidalCorrection.value(date)[3];
526 }
527 return interpolated;
528
529 }
530
531 /** Get the LoD (Length of Day) value.
532 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
533 * @param date date at which the value is desired
534 * @param <T> type of the filed elements
535 * @return LoD in seconds (0 if date is outside covered range)
536 * @since 9.0
537 */
538 public <T extends CalculusFieldElement<T>> T getLOD(final FieldAbsoluteDate<T> date) {
539
540 final AbsoluteDate aDate = date.toAbsoluteDate();
541
542 // check if there is data for date
543 if (!this.hasDataFor(aDate)) {
544 // no EOP data available for this date, we use a default null correction
545 return (tidalCorrection == null) ? date.getField().getZero() : tidalCorrection.value(date)[3];
546 }
547
548 // we have EOP data for date -> interpolate correction
549 T interpolated = interpolate(date, aDate, EOPEntry::getLOD);
550 if (tidalCorrection != null) {
551 interpolated = interpolated.add(tidalCorrection.value(date)[3]);
552 }
553
554 return interpolated;
555
556 }
557
558 /** Get the pole IERS Reference Pole correction.
559 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
560 * @param date date at which the correction is desired
561 * @return pole correction ({@link PoleCorrection#NULL_CORRECTION
562 * PoleCorrection.NULL_CORRECTION} if date is outside covered range)
563 */
564 public PoleCorrection getPoleCorrection(final AbsoluteDate date) {
565
566 // check if there is data for date
567 if (!this.hasDataFor(date)) {
568 // no EOP data available for this date, we use a default null correction
569 if (tidalCorrection == null) {
570 return PoleCorrection.NULL_CORRECTION;
571 } else {
572 final double[] correction = tidalCorrection.value(date);
573 return new PoleCorrection(correction[0], correction[1]);
574 }
575 }
576
577 // we have EOP data for date -> interpolate correction
578 final double[] interpolated = interpolate(date,
579 EOPEntry::getX, EOPEntry::getY,
580 EOPEntry::getXRate, EOPEntry::getYRate);
581 if (tidalCorrection != null) {
582 final double[] correction = tidalCorrection.value(date);
583 interpolated[0] += correction[0];
584 interpolated[1] += correction[1];
585 }
586 return new PoleCorrection(interpolated[0], interpolated[1]);
587
588 }
589
590 /** Get the pole IERS Reference Pole correction.
591 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
592 * @param date date at which the correction is desired
593 * @param <T> type of the field elements
594 * @return pole correction ({@link PoleCorrection#NULL_CORRECTION
595 * PoleCorrection.NULL_CORRECTION} if date is outside covered range)
596 */
597 public <T extends CalculusFieldElement<T>> FieldPoleCorrection<T> getPoleCorrection(final FieldAbsoluteDate<T> date) {
598
599 final AbsoluteDate aDate = date.toAbsoluteDate();
600
601 // check if there is data for date
602 if (!this.hasDataFor(aDate)) {
603 // no EOP data available for this date, we use a default null correction
604 if (tidalCorrection == null) {
605 return new FieldPoleCorrection<>(date.getField().getZero(), date.getField().getZero());
606 } else {
607 final T[] correction = tidalCorrection.value(date);
608 return new FieldPoleCorrection<>(correction[0], correction[1]);
609 }
610 }
611
612 // we have EOP data for date -> interpolate correction
613 final T[] interpolated = interpolate(date, aDate, EOPEntry::getX, EOPEntry::getY);
614 if (tidalCorrection != null) {
615 final T[] correction = tidalCorrection.value(date);
616 interpolated[0] = interpolated[0].add(correction[0]);
617 interpolated[1] = interpolated[1].add(correction[1]);
618 }
619 return new FieldPoleCorrection<>(interpolated[0], interpolated[1]);
620
621 }
622
623 /** Get the correction to the nutation parameters for equinox-based paradigm.
624 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
625 * @param date date at which the correction is desired
626 * @return nutation correction in longitude ΔΨ and in obliquity Δε
627 * (zero if date is outside covered range)
628 */
629 public double[] getEquinoxNutationCorrection(final AbsoluteDate date) {
630
631 // check if there is data for date
632 if (!this.hasDataFor(date)) {
633 // no EOP data available for this date, we use a default null correction
634 return new double[2];
635 }
636
637 // we have EOP data for date -> interpolate correction
638 return interpolate(date, EOPEntry::getDdPsi, EOPEntry::getDdEps);
639
640 }
641
642 /** Get the correction to the nutation parameters for equinox-based paradigm.
643 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
644 * @param date date at which the correction is desired
645 * @param <T> type of the field elements
646 * @return nutation correction in longitude ΔΨ and in obliquity Δε
647 * (zero if date is outside covered range)
648 */
649 public <T extends CalculusFieldElement<T>> T[] getEquinoxNutationCorrection(final FieldAbsoluteDate<T> date) {
650
651 final AbsoluteDate aDate = date.toAbsoluteDate();
652
653 // check if there is data for date
654 if (!this.hasDataFor(aDate)) {
655 // no EOP data available for this date, we use a default null correction
656 return MathArrays.buildArray(date.getField(), 2);
657 }
658
659 // we have EOP data for date -> interpolate correction
660 return interpolate(date, aDate, EOPEntry::getDdPsi, EOPEntry::getDdEps);
661
662 }
663
664 /** Get the correction to the nutation parameters for Non-Rotating Origin paradigm.
665 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
666 * @param date date at which the correction is desired
667 * @return nutation correction in Celestial Intermediate Pole coordinates
668 * δX and δY (zero if date is outside covered range)
669 */
670 public double[] getNonRotatingOriginNutationCorrection(final AbsoluteDate date) {
671
672 // check if there is data for date
673 if (!this.hasDataFor(date)) {
674 // no EOP data available for this date, we use a default null correction
675 return new double[2];
676 }
677
678 // we have EOP data for date -> interpolate correction
679 return interpolate(date, EOPEntry::getDx, EOPEntry::getDy);
680
681 }
682
683 /** Get the correction to the nutation parameters for Non-Rotating Origin paradigm.
684 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
685 * @param date date at which the correction is desired
686 * @param <T> type of the filed elements
687 * @return nutation correction in Celestial Intermediate Pole coordinates
688 * δX and δY (zero if date is outside covered range)
689 */
690 public <T extends CalculusFieldElement<T>> T[] getNonRotatingOriginNutationCorrection(final FieldAbsoluteDate<T> date) {
691
692 final AbsoluteDate aDate = date.toAbsoluteDate();
693
694 // check if there is data for date
695 if (!this.hasDataFor(aDate)) {
696 // no EOP data available for this date, we use a default null correction
697 return MathArrays.buildArray(date.getField(), 2);
698 }
699
700 // we have EOP data for date -> interpolate correction
701 return interpolate(date, aDate, EOPEntry::getDx, EOPEntry::getDy);
702
703 }
704
705 /** Get the ITRF version.
706 * @param date date at which the value is desired
707 * @return ITRF version of the EOP covering the specified date
708 * @since 9.2
709 */
710 public ITRFVersion getITRFVersion(final AbsoluteDate date) {
711
712 // check if there is data for date
713 if (!this.hasDataFor(date)) {
714 // no EOP data available for this date, we use a default ITRF 2014
715 return ITRFVersion.ITRF_2014;
716 }
717
718 try {
719 // we have EOP data for date
720 final Optional<EOPEntry> first = getNeighbors(date, 1).findFirst();
721 return first.isPresent() ? first.get().getITRFType() : ITRFVersion.ITRF_2014;
722
723 } catch (TimeStampedCacheException tce) {
724 // this should not happen because of date check performed at start
725 throw new OrekitInternalError(tce);
726 }
727
728 }
729
730 /** Get the EOP data type.
731 * @param date date at which the value is desired
732 * @return data type of the EOP covering the specified date
733 * @since 13.1.1
734 */
735 public EopDataType getEopDataType(final AbsoluteDate date) {
736
737 // check if there is data for date
738 if (!this.hasDataFor(date)) {
739 // no EOP data available for this date, data type is unknown
740 return EopDataType.UNKNOWN;
741 }
742
743 // we have EOP data for date
744 final Optional<EOPEntry> first = getNeighbors(date, 1).findFirst();
745 return first.isPresent() ? first.get().getEopDataType() : EopDataType.UNKNOWN;
746
747 }
748
749 /** Check Earth orientation parameters continuity.
750 * @param maxGap maximal allowed gap between entries (in seconds)
751 */
752 public void checkEOPContinuity(final double maxGap) {
753 TimeStamped preceding = null;
754 for (final TimeStamped current : this.cache.getAll()) {
755
756 // compare the dates of preceding and current entries
757 if (preceding != null && (current.getDate().durationFrom(preceding.getDate())) > maxGap) {
758 throw new OrekitException(OrekitMessages.MISSING_EARTH_ORIENTATION_PARAMETERS_BETWEEN_DATES_GAP,
759 preceding.getDate(), current.getDate(),
760 current.getDate().durationFrom(preceding.getDate()));
761 }
762
763 // prepare next iteration
764 preceding = current;
765
766 }
767 }
768
769 /**
770 * Check if the cache has data for the given date using
771 * {@link #getStartDate()} and {@link #getEndDate()}.
772 *
773 * @param date the requested date
774 * @return true if the {@link #cache} has data for the requested date, false
775 * otherwise.
776 */
777 protected boolean hasDataFor(final AbsoluteDate date) {
778 /*
779 * when there is no EOP data, short circuit getStartDate, which will
780 * throw an exception
781 */
782 return this.hasData && this.getStartDate().compareTo(date) <= 0 &&
783 date.compareTo(this.getEndDate()) <= 0;
784 }
785
786 /** Get a non-modifiable view of the EOP entries.
787 * @return non-modifiable view of the EOP entries
788 */
789 public List<EOPEntry> getEntries() {
790 return cache.getAll();
791 }
792
793 /** Interpolate a single EOP component.
794 * <p>
795 * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
796 * </p>
797 * @param date interpolation date
798 * @param selector selector for EOP entry component
799 * @return interpolated value
800 */
801 private double interpolate(final AbsoluteDate date, final Function<EOPEntry, Double> selector) {
802 try {
803 final HermiteInterpolator interpolator = new HermiteInterpolator();
804 getNeighbors(date, interpolationDegree + 1).
805 forEach(entry -> interpolator.addSamplePoint(entry.getDate().durationFrom(date),
806 new double[] {
807 selector.apply(entry)
808 }));
809 return interpolator.value(0)[0];
810 } catch (TimeStampedCacheException tce) {
811 // this should not happen because of date check performed by caller
812 throw new OrekitInternalError(tce);
813 }
814 }
815
816 /** Interpolate a single EOP component.
817 * <p>
818 * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
819 * </p>
820 * @param date interpolation date
821 * @param aDate interpolation date, as an {@link AbsoluteDate}
822 * @param selector selector for EOP entry component
823 * @param <T> type of the field elements
824 * @return interpolated value
825 */
826 private <T extends CalculusFieldElement<T>> T interpolate(final FieldAbsoluteDate<T> date,
827 final AbsoluteDate aDate,
828 final Function<EOPEntry, Double> selector) {
829 try {
830 final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
831 final T[] y = MathArrays.buildArray(date.getField(), 1);
832 final T zero = date.getField().getZero();
833 final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
834 // for example removing derivatives
835 // if T was DerivativeStructure
836 getNeighbors(aDate, interpolationDegree + 1).
837 forEach(entry -> {
838 y[0] = zero.newInstance(selector.apply(entry));
839 interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
840 });
841 return interpolator.value(date.durationFrom(central))[0]; // here, we introduce derivatives again (in DerivativeStructure case)
842 } catch (TimeStampedCacheException tce) {
843 // this should not happen because of date check performed by caller
844 throw new OrekitInternalError(tce);
845 }
846 }
847
848 /** Interpolate two EOP components.
849 * <p>
850 * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
851 * </p>
852 * @param date interpolation date
853 * @param selector1 selector for first EOP entry component
854 * @param selector2 selector for second EOP entry component
855 * @return interpolated value
856 */
857 private double[] interpolate(final AbsoluteDate date,
858 final Function<EOPEntry, Double> selector1,
859 final Function<EOPEntry, Double> selector2) {
860 try {
861 final HermiteInterpolator interpolator = new HermiteInterpolator();
862 getNeighbors(date, interpolationDegree + 1).
863 forEach(entry -> interpolator.addSamplePoint(entry.getDate().durationFrom(date),
864 new double[] {
865 selector1.apply(entry),
866 selector2.apply(entry)
867 }));
868 return interpolator.value(0);
869 } catch (TimeStampedCacheException tce) {
870 // this should not happen because of date check performed by caller
871 throw new OrekitInternalError(tce);
872 }
873 }
874
875 /** Interpolate two EOP components.
876 * <p>
877 * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
878 * </p>
879 * @param date interpolation date
880 * @param selector1 selector for first EOP entry component
881 * @param selector2 selector for second EOP entry component
882 * @param selector1Rate selector for first EOP entry component rate
883 * @param selector2Rate selector for second EOP entry component rate
884 * @return interpolated value
885 * @since 12.0
886 */
887 private double[] interpolate(final AbsoluteDate date,
888 final Function<EOPEntry, Double> selector1,
889 final Function<EOPEntry, Double> selector2,
890 final Function<EOPEntry, Double> selector1Rate,
891 final Function<EOPEntry, Double> selector2Rate) {
892 try {
893 final HermiteInterpolator interpolator = new HermiteInterpolator();
894 getNeighbors(date, (interpolationDegree + 1) / 2).
895 forEach(entry -> interpolator.addSamplePoint(entry.getDate().durationFrom(date),
896 new double[] {
897 selector1.apply(entry),
898 selector2.apply(entry)
899 },
900 new double[] {
901 selector1Rate.apply(entry),
902 selector2Rate.apply(entry)
903 }));
904 return interpolator.value(0);
905 } catch (TimeStampedCacheException tce) {
906 // this should not happen because of date check performed by caller
907 throw new OrekitInternalError(tce);
908 }
909 }
910
911 /** Interpolate two EOP components.
912 * <p>
913 * This method should be called <em>only</em> when {@link #hasDataFor(AbsoluteDate)} returns true.
914 * </p>
915 * @param date interpolation date
916 * @param aDate interpolation date, as an {@link AbsoluteDate}
917 * @param selector1 selector for first EOP entry component
918 * @param selector2 selector for second EOP entry component
919 * @param <T> type of the field elements
920 * @return interpolated value
921 */
922 private <T extends CalculusFieldElement<T>> T[] interpolate(final FieldAbsoluteDate<T> date,
923 final AbsoluteDate aDate,
924 final Function<EOPEntry, Double> selector1,
925 final Function<EOPEntry, Double> selector2) {
926 try {
927 final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
928 final T[] y = MathArrays.buildArray(date.getField(), 2);
929 final T zero = date.getField().getZero();
930 final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
931 // for example removing derivatives
932 // if T was DerivativeStructure
933 getNeighbors(aDate, interpolationDegree + 1).
934 forEach(entry -> {
935 y[0] = zero.newInstance(selector1.apply(entry));
936 y[1] = zero.newInstance(selector2.apply(entry));
937 interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
938 });
939 return interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
940 } catch (TimeStampedCacheException tce) {
941 // this should not happen because of date check performed by caller
942 throw new OrekitInternalError(tce);
943 }
944 }
945
946 /** Check if some derivatives are missing.
947 * @param raw raw EOP history
948 * @return true if some derivatives are missing
949 * @since 12.0
950 */
951 private boolean missSomeDerivatives(final Collection<? extends EOPEntry> raw) {
952 for (final EOPEntry entry : raw) {
953 if (Double.isNaN(entry.getLOD() + entry.getXRate() + entry.getYRate())) {
954 return true;
955 }
956 }
957 return false;
958 }
959
960 /** Fix missing derivatives.
961 * @param entry EOP entry to fix
962 * @param rawCache raw EOP history cache
963 * @return fixed entry
964 * @since 12.0
965 */
966 private EOPEntry fixDerivatives(final EOPEntry entry, final ImmutableTimeStampedCache<EOPEntry> rawCache) {
967
968 // helper function to differentiate some EOP parameters
969 final BiFunction<EOPEntry, Function<EOPEntry, Double>, Double> differentiator =
970 (e, selector) -> {
971 final HermiteInterpolator interpolator = new HermiteInterpolator();
972 rawCache.getNeighbors(e.getDate()).
973 forEach(f -> interpolator.addSamplePoint(f.getDate().durationFrom(e.getDate()),
974 new double[] {
975 selector.apply(f)
976 }));
977 return interpolator.derivatives(0.0, 1)[1][0];
978 };
979
980 if (Double.isNaN(entry.getLOD() + entry.getXRate() + entry.getYRate())) {
981 final double lod = Double.isNaN(entry.getLOD()) ?
982 -Constants.JULIAN_DAY * differentiator.apply(entry, EOPEntry::getUT1MinusUTC) :
983 entry.getLOD();
984 final double xRate = Double.isNaN(entry.getXRate()) ?
985 differentiator.apply(entry, EOPEntry::getX) :
986 entry.getXRate();
987 final double yRate = Double.isNaN(entry.getYRate()) ?
988 differentiator.apply(entry, EOPEntry::getY) :
989 entry.getYRate();
990 return new EOPEntry(entry.getMjd(),
991 entry.getUT1MinusUTC(), lod,
992 entry.getX(), entry.getY(), xRate, yRate,
993 entry.getDdPsi(), entry.getDdEps(),
994 entry.getDx(), entry.getDy(),
995 entry.getITRFType(), entry.getDate(), entry.getEopDataType(),
996 entry.getDtOrigin(), entry.getNutOrigin(), entry.getCipOrigin());
997 } else {
998 // the entry already has all derivatives
999 return entry;
1000 }
1001 }
1002
1003 /** Internal class for caching tidal correction. */
1004 private static class TidalCorrectionEntry implements TimeStamped {
1005
1006 /** Entry date. */
1007 private final AbsoluteDate date;
1008
1009 /** Correction. */
1010 private final double[] correction;
1011
1012 /** Simple constructor.
1013 * @param date entry date
1014 * @param correction correction on the EOP parameters (xp, yp, ut1, lod)
1015 */
1016 TidalCorrectionEntry(final AbsoluteDate date, final double[] correction) {
1017 this.date = date;
1018 this.correction = correction;
1019 }
1020
1021 /** {@inheritDoc} */
1022 @Override
1023 public AbsoluteDate getDate() {
1024 return date;
1025 }
1026
1027 }
1028
1029 /** Local generator for thread-safe cache. */
1030 private static class CachedCorrection
1031 implements TimeVectorFunction, TimeStampedGenerator<TidalCorrectionEntry> {
1032
1033 /** Correction to apply to EOP (may be null). */
1034 private final TimeVectorFunction tidalCorrection;
1035
1036 /** Step between generated entries. */
1037 private final double step;
1038
1039 /** Tidal corrections entries cache. */
1040 private final TimeStampedCache<TidalCorrectionEntry> cache;
1041
1042 /** Simple constructor.
1043 * @param tidalCorrection function computing the tidal correction
1044 */
1045 CachedCorrection(final TimeVectorFunction tidalCorrection) {
1046 this.step = 60. * 60;
1047 this.tidalCorrection = tidalCorrection;
1048 this.cache =
1049 new GenericTimeStampedCache<>(8,
1050 OrekitConfiguration.getCacheSlotsNumber(),
1051 Constants.JULIAN_DAY * 30,
1052 Constants.JULIAN_DAY,
1053 this);
1054 }
1055
1056 /** {@inheritDoc} */
1057 @Override
1058 public double[] value(final AbsoluteDate date) {
1059 try {
1060 // set up an interpolator
1061 final HermiteInterpolator interpolator = new HermiteInterpolator();
1062 cache.getNeighbors(date).forEach(entry -> interpolator.addSamplePoint(entry.date.durationFrom(date), entry.correction));
1063
1064 // interpolate to specified date
1065 return interpolator.value(0.0);
1066 } catch (TimeStampedCacheException tsce) {
1067 // this should never happen
1068 throw new OrekitInternalError(tsce);
1069 }
1070 }
1071
1072 /** {@inheritDoc} */
1073 @Override
1074 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
1075 try {
1076
1077 final AbsoluteDate aDate = date.toAbsoluteDate();
1078
1079 final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
1080 final T[] y = MathArrays.buildArray(date.getField(), 4);
1081 final T zero = date.getField().getZero();
1082 final FieldAbsoluteDate<T> central = new FieldAbsoluteDate<>(aDate, zero); // here, we attempt to get a constant date,
1083 // for example removing derivatives
1084 // if T was DerivativeStructure
1085 cache.getNeighbors(aDate).forEach(entry -> {
1086 for (int i = 0; i < y.length; ++i) {
1087 y[i] = zero.newInstance(entry.correction[i]);
1088 }
1089 interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
1090 });
1091
1092 // interpolate to specified date
1093 return interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
1094
1095 } catch (TimeStampedCacheException tsce) {
1096 // this should never happen
1097 throw new OrekitInternalError(tsce);
1098 }
1099 }
1100
1101 /** {@inheritDoc} */
1102 @Override
1103 public List<TidalCorrectionEntry> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {
1104
1105 final List<TidalCorrectionEntry> generated = new ArrayList<>();
1106
1107 if (existingDate == null) {
1108
1109 // no prior existing entries, just generate a first set
1110 for (int i = -cache.getMaxNeighborsSize() / 2; generated.size() < cache.getMaxNeighborsSize(); ++i) {
1111 final AbsoluteDate t = date.shiftedBy(i * step);
1112 generated.add(new TidalCorrectionEntry(t, tidalCorrection.value(t)));
1113 }
1114
1115 } else {
1116
1117 // some entries have already been generated
1118 // add the missing ones up to specified date
1119
1120 AbsoluteDate t = existingDate;
1121 if (date.compareTo(t) > 0) {
1122 // forward generation
1123 do {
1124 t = t.shiftedBy(step);
1125 generated.add(new TidalCorrectionEntry(t, tidalCorrection.value(t)));
1126 } while (t.compareTo(date) <= 0);
1127 } else {
1128 // backward generation
1129 do {
1130 t = t.shiftedBy(-step);
1131 generated.addFirst(new TidalCorrectionEntry(t, tidalCorrection.value(t)));
1132 } while (t.compareTo(date) >= 0);
1133 }
1134 }
1135
1136 // return the generated transforms
1137 return generated;
1138
1139 }
1140 }
1141
1142 }