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