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