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