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