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