1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.bodies;
18
19 import java.io.IOException;
20 import java.io.InputStream;
21 import java.nio.charset.StandardCharsets;
22 import java.text.ParseException;
23 import java.util.ArrayList;
24 import java.util.HashMap;
25 import java.util.List;
26 import java.util.Map;
27 import java.util.SortedSet;
28 import java.util.TreeSet;
29 import java.util.concurrent.TimeUnit;
30 import java.util.concurrent.atomic.AtomicReference;
31
32 import org.hipparchus.CalculusFieldElement;
33 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
34 import org.hipparchus.geometry.euclidean.threed.Vector3D;
35 import org.hipparchus.util.FastMath;
36 import org.orekit.annotation.DefaultDataContext;
37 import org.orekit.data.AbstractSelfFeedingLoader;
38 import org.orekit.data.DataContext;
39 import org.orekit.data.DataLoader;
40 import org.orekit.data.DataProvidersManager;
41 import org.orekit.errors.OrekitException;
42 import org.orekit.errors.OrekitInternalError;
43 import org.orekit.errors.OrekitMessages;
44 import org.orekit.errors.TimeStampedCacheException;
45 import org.orekit.frames.Frame;
46 import org.orekit.time.AbsoluteDate;
47 import org.orekit.time.ChronologicalComparator;
48 import org.orekit.time.DateComponents;
49 import org.orekit.time.FieldAbsoluteDate;
50 import org.orekit.time.TimeComponents;
51 import org.orekit.time.TimeOffset;
52 import org.orekit.time.TimeScale;
53 import org.orekit.time.TimeScales;
54 import org.orekit.utils.Constants;
55 import org.orekit.utils.FieldPVCoordinates;
56 import org.orekit.utils.OrekitConfiguration;
57 import org.orekit.utils.PVCoordinates;
58 import org.orekit.utils.GenericTimeStampedCache;
59 import org.orekit.utils.TimeStampedGenerator;
60 import org.orekit.utils.units.UnitsConverter;
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80 public class JPLEphemeridesLoader extends AbstractSelfFeedingLoader
81 implements CelestialBodyLoader {
82
83
84 public static final String DEFAULT_DE_SUPPORTED_NAMES = "^[lu]nx([mp](\\d{4,5}))+\\.(?:4\\d\\d)$";
85
86
87 public static final String DEFAULT_DE_2021_SUPPORTED_NAMES = "^linux_([mp](\\d{4,5}))+\\.(?:4\\d\\d)$";
88
89
90 public static final String DEFAULT_INPOP_SUPPORTED_NAMES = "^inpop.*\\.dat$";
91
92
93 private static final TimeOffset FIFTY_DAYS =
94 new TimeOffset(50, TimeUnit.DAYS);
95
96
97 private static final int INPOP_DE_NUMBER = 100;
98
99
100 private static final int CONSTANTS_MAX_NUMBER = 400;
101
102
103 private static final int HEADER_EPHEMERIS_TYPE_OFFSET = 2840;
104
105
106 private static final int HEADER_RECORD_SIZE_OFFSET = 2856;
107
108
109 private static final int HEADER_START_EPOCH_OFFSET = 2652;
110
111
112 private static final int HEADER_END_EPOCH_OFFSET = 2660;
113
114
115 private static final int HEADER_ASTRONOMICAL_UNIT_OFFSET = 2680;
116
117
118 private static final int HEADER_EM_RATIO_OFFSET = 2688;
119
120
121 private static final int HEADER_CHEBISHEV_INDICES_OFFSET = 2696;
122
123
124 private static final int HEADER_LIBRATION_INDICES_OFFSET = 2844;
125
126
127 private static final int HEADER_CHUNK_DURATION_OFFSET = 2668;
128
129
130 private static final int HEADER_CONSTANTS_NAMES_OFFSET = 252;
131
132
133 private static final int HEADER_CONSTANTS_VALUES_OFFSET = 0;
134
135
136 private static final int DATA_START_RANGE_OFFSET = 0;
137
138
139 private static final int DATE_END_RANGE_OFFSET = 8;
140
141
142 private static final String CONSTANT_AU = "AU";
143
144
145 private static final String CONSTANT_EMRAT = "EMRAT";
146
147
148 public enum EphemerisType {
149
150
151 SOLAR_SYSTEM_BARYCENTER,
152
153
154 SUN,
155
156
157 MERCURY,
158
159
160 VENUS,
161
162
163 EARTH_MOON,
164
165
166 EARTH,
167
168
169 MOON,
170
171
172 MARS,
173
174
175 JUPITER,
176
177
178 SATURN,
179
180
181 URANUS,
182
183
184 NEPTUNE,
185
186
187 PLUTO,
188
189
190 LIBRATION
191
192 }
193
194
195 public interface RawPVProvider {
196
197
198
199
200
201 PVCoordinates getRawPV(AbsoluteDate date);
202
203
204
205
206
207 default Vector3D getRawPosition(final AbsoluteDate date) {
208 return getRawPV(date).getPosition();
209 }
210
211
212
213
214
215
216 <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(FieldAbsoluteDate<T> date);
217
218
219
220
221
222
223 default <T extends CalculusFieldElement<T>> FieldVector3D<T> getRawPosition(final FieldAbsoluteDate<T> date) {
224 return getRawPV(date).getPosition();
225 }
226
227 }
228
229
230 private final GenericTimeStampedCache<PosVelChebyshev> ephemerides;
231
232
233 private final AtomicReference<Map<String, Double>> constants;
234
235
236 private final EphemerisType generateType;
237
238
239 private final EphemerisType loadType;
240
241
242 private final TimeScales timeScales;
243
244
245 private final Frame gcrf;
246
247
248 private final CelestialBodies celestialBodies;
249
250
251 private AbsoluteDate startEpoch;
252
253
254 private AbsoluteDate finalEpoch;
255
256
257 private double maxChunksDuration;
258
259
260 private TimeOffset chunksDuration;
261
262
263 private int firstIndex;
264
265
266 private int coeffs;
267
268
269 private int chunks;
270
271
272 private int components;
273
274
275 private double positionUnit;
276
277
278 private TimeScale timeScale;
279
280
281 private boolean bigEndian;
282
283
284
285
286
287
288
289
290
291 @DefaultDataContext
292 public JPLEphemeridesLoader(final String supportedNames, final EphemerisType generateType) {
293 this(supportedNames, generateType,
294 DataContext.getDefault().getDataProvidersManager(),
295 DataContext.getDefault().getTimeScales(),
296 DataContext.getDefault().getFrames().getGCRF(),
297 DataContext.getDefault().getCelestialBodies());
298 }
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319 public JPLEphemeridesLoader(final String supportedNames,
320 final EphemerisType generateType,
321 final DataProvidersManager dataProvidersManager,
322 final TimeScales timeScales,
323 final Frame gcrf,
324 final CelestialBodies celestialBodies) {
325 super(supportedNames, dataProvidersManager);
326
327 this.timeScales = timeScales;
328 this.gcrf = gcrf;
329 this.celestialBodies = celestialBodies;
330 constants = new AtomicReference<>();
331
332 this.generateType = generateType;
333 if (generateType == EphemerisType.SOLAR_SYSTEM_BARYCENTER) {
334 loadType = EphemerisType.EARTH_MOON;
335 } else if (generateType == EphemerisType.EARTH_MOON) {
336 loadType = EphemerisType.MOON;
337 } else {
338 loadType = generateType;
339 }
340
341 ephemerides = new GenericTimeStampedCache<>(
342 2, OrekitConfiguration.getCacheSlotsNumber(),
343 Double.POSITIVE_INFINITY, FIFTY_DAYS.toDouble(),
344 new EphemerisParser());
345 maxChunksDuration = Double.NaN;
346 chunksDuration = null;
347
348 }
349
350
351
352
353
354 @Override
355 public CelestialBody loadCelestialBody(final String name) {
356
357 final double gm = getLoadedGravitationalCoefficient(generateType);
358 final IAUPole iauPole = PredefinedIAUPoles
359 .getIAUPole(generateType, timeScales);
360 final double scale;
361 final Frame definingFrameAlignedWithICRF;
362 final RawPVProvider rawPVProvider;
363 switch (generateType) {
364 case SOLAR_SYSTEM_BARYCENTER : {
365 scale = -1.0;
366 definingFrameAlignedWithICRF = celestialBodies
367 .getBody(CelestialBodyFactory.EARTH_MOON)
368 .getIcrfAlignedFrame();
369 rawPVProvider = new EphemerisRawPVProvider();
370 break;
371 }
372 case EARTH_MOON :
373 scale = 1.0 / (1.0 + getLoadedEarthMoonMassRatio());
374 definingFrameAlignedWithICRF = gcrf;
375 rawPVProvider = new EphemerisRawPVProvider();
376 break;
377 case EARTH :
378 scale = 1.0;
379 definingFrameAlignedWithICRF = gcrf;
380 rawPVProvider = new ZeroRawPVProvider();
381 break;
382 case MOON :
383 scale = 1.0;
384 definingFrameAlignedWithICRF = gcrf;
385 rawPVProvider = new EphemerisRawPVProvider();
386 break;
387 default : {
388 scale = 1.0;
389 definingFrameAlignedWithICRF = celestialBodies
390 .getBody(CelestialBodyFactory.SOLAR_SYSTEM_BARYCENTER)
391 .getIcrfAlignedFrame();
392 rawPVProvider = new EphemerisRawPVProvider();
393 }
394 }
395
396
397 return new JPLCelestialBody(name, getSupportedNames(), generateType, rawPVProvider,
398 gm, scale, iauPole, definingFrameAlignedWithICRF);
399
400 }
401
402
403
404
405
406 public JPLLibration loadLibration() {
407 final RawPVProvider rawPVProvider = new EphemerisRawPVProvider();
408
409
410 return new JPLLibration(rawPVProvider);
411 }
412
413
414
415
416 public double getLoadedAstronomicalUnit() {
417 return UnitsConverter.KILOMETRES_TO_METRES.convert(getLoadedConstant(CONSTANT_AU));
418 }
419
420
421
422
423 public double getLoadedEarthMoonMassRatio() {
424 return getLoadedConstant(CONSTANT_EMRAT);
425 }
426
427
428
429
430
431 public double getLoadedGravitationalCoefficient(final EphemerisType body) {
432
433
434 final double rawGM;
435 switch (body) {
436 case SOLAR_SYSTEM_BARYCENTER :
437 return getLoadedGravitationalCoefficient(EphemerisType.SUN) +
438 getLoadedGravitationalCoefficient(EphemerisType.MERCURY) +
439 getLoadedGravitationalCoefficient(EphemerisType.VENUS) +
440 getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) +
441 getLoadedGravitationalCoefficient(EphemerisType.MARS) +
442 getLoadedGravitationalCoefficient(EphemerisType.JUPITER) +
443 getLoadedGravitationalCoefficient(EphemerisType.SATURN) +
444 getLoadedGravitationalCoefficient(EphemerisType.URANUS) +
445 getLoadedGravitationalCoefficient(EphemerisType.NEPTUNE) +
446 getLoadedGravitationalCoefficient(EphemerisType.PLUTO);
447 case SUN :
448 rawGM = getLoadedConstant("GMS", "GM_Sun");
449 break;
450 case MERCURY :
451 rawGM = getLoadedConstant("GM1", "GM_Mer");
452 break;
453 case VENUS :
454 rawGM = getLoadedConstant("GM2", "GM_Ven");
455 break;
456 case EARTH_MOON :
457 rawGM = getLoadedConstant("GMB", "GM_EMB");
458 break;
459 case EARTH :
460 return getLoadedEarthMoonMassRatio() *
461 getLoadedGravitationalCoefficient(EphemerisType.MOON);
462 case MOON :
463 return getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) /
464 (1.0 + getLoadedEarthMoonMassRatio());
465 case MARS :
466 rawGM = getLoadedConstant("GM4", "GM_Mar");
467 break;
468 case JUPITER :
469 rawGM = getLoadedConstant("GM5", "GM_Jup");
470 break;
471 case SATURN :
472 rawGM = getLoadedConstant("GM6", "GM_Sat");
473 break;
474 case URANUS :
475 rawGM = getLoadedConstant("GM7", "GM_Ura");
476 break;
477 case NEPTUNE :
478 rawGM = getLoadedConstant("GM8", "GM_Nep");
479 break;
480 case PLUTO :
481 rawGM = getLoadedConstant("GM9", "GM_Plu");
482 break;
483 default :
484 throw new OrekitInternalError(null);
485 }
486
487 final double au = getLoadedAstronomicalUnit();
488 return rawGM * au * au * au / (Constants.JULIAN_DAY * Constants.JULIAN_DAY);
489
490 }
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505 public double getLoadedConstant(final String... names) {
506
507
508 Map<String, Double> map = constants.get();
509 if (map == null) {
510 final ConstantsParser parser = new ConstantsParser();
511 if (!feed(parser)) {
512 throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
513 }
514 map = parser.getConstants();
515 constants.compareAndSet(null, map);
516 }
517
518 for (final String name : names) {
519 if (map.containsKey(name)) {
520 return map.get(name);
521 }
522 }
523
524 return Double.NaN;
525
526 }
527
528
529
530
531 public double getMaxChunksDuration() {
532 return maxChunksDuration;
533 }
534
535
536
537
538
539 private void parseFirstHeaderRecord(final byte[] record, final String name) {
540
541
542 final int deNum = extractInt(record, HEADER_EPHEMERIS_TYPE_OFFSET);
543
544
545
546
547 components = 3;
548 positionUnit = UnitsConverter.KILOMETRES_TO_METRES.convert(1.0);
549 timeScale = timeScales.getTDB();
550
551 if (deNum == INPOP_DE_NUMBER) {
552
553 final double format = getLoadedConstant("FORMAT");
554 if (!Double.isNaN(format) && (int) FastMath.IEEEremainder(format, 10) != 1) {
555 components = 6;
556 }
557
558
559 final double unite = getLoadedConstant("UNITE");
560 if (!Double.isNaN(unite) && (int) unite == 0) {
561 positionUnit = getLoadedAstronomicalUnit();
562 }
563
564
565 final double timesc = getLoadedConstant("TIMESC");
566 if (!Double.isNaN(timesc) && (int) timesc == 1) {
567 timeScale = timeScales.getTCB();
568 }
569
570 }
571
572
573 startEpoch = extractDate(record, HEADER_START_EPOCH_OFFSET);
574 finalEpoch = extractDate(record, HEADER_END_EPOCH_OFFSET);
575 boolean ok = finalEpoch.compareTo(startEpoch) > 0;
576
577 if (loadType == EphemerisType.LIBRATION) {
578
579 firstIndex = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET);
580 coeffs = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 4);
581 chunks = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 8);
582 ok = ok && firstIndex >= 0 && coeffs >= 0 && chunks >= 0;
583 positionUnit = 1.0;
584 }
585 else {
586
587 for (int i = 0; i < 12; ++i) {
588 final int row1 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 12 * i);
589 final int row2 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 4 + 12 * i);
590 final int row3 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 8 + 12 * i);
591 ok = ok && row1 >= 0 && row2 >= 0 && row3 >= 0;
592 if (i == 0 && loadType == EphemerisType.MERCURY ||
593 i == 1 && loadType == EphemerisType.VENUS ||
594 i == 2 && loadType == EphemerisType.EARTH_MOON ||
595 i == 3 && loadType == EphemerisType.MARS ||
596 i == 4 && loadType == EphemerisType.JUPITER ||
597 i == 5 && loadType == EphemerisType.SATURN ||
598 i == 6 && loadType == EphemerisType.URANUS ||
599 i == 7 && loadType == EphemerisType.NEPTUNE ||
600 i == 8 && loadType == EphemerisType.PLUTO ||
601 i == 9 && loadType == EphemerisType.MOON ||
602 i == 10 && loadType == EphemerisType.SUN) {
603 firstIndex = row1;
604 coeffs = row2;
605 chunks = row3;
606 }
607 }
608 }
609
610
611 final double timeSpan = extractDouble(record, HEADER_CHUNK_DURATION_OFFSET);
612 ok = ok && timeSpan > 0 && timeSpan < 100;
613 chunksDuration = new TimeOffset(timeSpan).divide(chunks).multiply(86400L);
614 if (Double.isNaN(maxChunksDuration)) {
615 maxChunksDuration = chunksDuration.toDouble();
616 } else {
617 maxChunksDuration = FastMath
618 .max(maxChunksDuration, chunksDuration.toDouble());
619 }
620
621
622 if (!ok) {
623 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
624 }
625
626 }
627
628
629
630
631
632
633
634 private byte[] readFirstRecord(final InputStream input, final String name)
635 throws IOException {
636
637
638 final byte[] firstPart = new byte[HEADER_RECORD_SIZE_OFFSET + 4];
639 if (!readInRecord(input, firstPart, 0)) {
640 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
641 }
642
643
644 detectEndianess(firstPart);
645
646
647 final int deNum = extractInt(firstPart, HEADER_EPHEMERIS_TYPE_OFFSET);
648
649
650 final int recordSize;
651
652 if (deNum == INPOP_DE_NUMBER) {
653
654 recordSize = extractInt(firstPart, HEADER_RECORD_SIZE_OFFSET) << 3;
655 } else {
656
657 recordSize = computeRecordSize(firstPart, name);
658 }
659
660 if (recordSize <= 0) {
661 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
662 }
663
664
665 final int start = firstPart.length;
666 final byte[] record = new byte[recordSize];
667 System.arraycopy(firstPart, 0, record, 0, firstPart.length);
668 if (!readInRecord(input, record, start)) {
669 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
670 }
671
672 return record;
673
674 }
675
676
677
678
679
680
681 private Map<String, Double> parseConstants(final byte[] first, final byte[] second) {
682
683 final Map<String, Double> map = new HashMap<>();
684
685 for (int i = 0; i < CONSTANTS_MAX_NUMBER; ++i) {
686
687
688 final String constantName = extractString(first, HEADER_CONSTANTS_NAMES_OFFSET + i * 6, 6);
689 if (constantName.length() == 0) {
690
691 break;
692 }
693 final double constantValue = extractDouble(second, HEADER_CONSTANTS_VALUES_OFFSET + 8 * i);
694 map.put(constantName, constantValue);
695 }
696
697
698
699 if (!map.containsKey(CONSTANT_AU)) {
700 map.put(CONSTANT_AU, extractDouble(first, HEADER_ASTRONOMICAL_UNIT_OFFSET));
701 }
702
703 if (!map.containsKey(CONSTANT_EMRAT)) {
704 map.put(CONSTANT_EMRAT, extractDouble(first, HEADER_EM_RATIO_OFFSET));
705 }
706
707 return map;
708
709 }
710
711
712
713
714
715
716
717
718 private boolean readInRecord(final InputStream input, final byte[] record, final int start)
719 throws IOException {
720 int index = start;
721 while (index != record.length) {
722 final int n = input.read(record, index, record.length - index);
723 if (n < 0) {
724 return false;
725 }
726 index += n;
727 }
728 return true;
729 }
730
731
732
733
734
735 private void detectEndianess(final byte[] record) {
736
737
738 bigEndian = true;
739
740
741
742 final long deNum = extractInt(record, HEADER_EPHEMERIS_TYPE_OFFSET) & 0xffffffffL;
743
744
745
746 if (deNum > (1 << 15)) {
747 bigEndian = false;
748 }
749
750 }
751
752
753
754
755
756
757 private int computeRecordSize(final byte[] record, final String name) {
758
759 int recordSize = 0;
760 boolean ok = true;
761
762 final int nComp = 3;
763
764
765
766 for (int j = 0; j < 12; j++) {
767 final int nCompCur = (j == 11) ? 2 : nComp;
768
769
770 final int idx = HEADER_CHEBISHEV_INDICES_OFFSET + j * nComp * 4;
771 final int coeffPtr1 = extractInt(record, idx + 4);
772 final int coeffPtr2 = extractInt(record, idx + 8);
773
774
775 ok = ok && (coeffPtr1 >= 0 || coeffPtr2 >= 0);
776
777 recordSize += coeffPtr1 * coeffPtr2 * nCompCur;
778 }
779
780
781
782 final int libratPtr1 = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 4);
783 final int libratPtr2 = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 8);
784
785
786 ok = ok && (libratPtr1 >= 0 || libratPtr2 >= 0);
787
788 recordSize += libratPtr1 * libratPtr2 * nComp + 2;
789 recordSize <<= 3;
790
791 if (!ok || recordSize <= 0) {
792 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
793 }
794
795 return recordSize;
796
797 }
798
799
800
801
802
803
804 private AbsoluteDate extractDate(final byte[] record, final int offset) {
805
806 final double t = extractDouble(record, offset);
807 int jDay = (int) FastMath.floor(t);
808 double seconds = (t + 0.5 - jDay) * Constants.JULIAN_DAY;
809 if (seconds >= Constants.JULIAN_DAY) {
810 ++jDay;
811 seconds -= Constants.JULIAN_DAY;
812 }
813 return new AbsoluteDate(new DateComponents(DateComponents.JULIAN_EPOCH, jDay),
814 new TimeComponents(seconds), timeScale);
815 }
816
817
818
819
820
821
822
823
824 private double extractDouble(final byte[] record, final int offset) {
825 final long l8 = ((long) record[offset + 0]) & 0xffl;
826 final long l7 = ((long) record[offset + 1]) & 0xffl;
827 final long l6 = ((long) record[offset + 2]) & 0xffl;
828 final long l5 = ((long) record[offset + 3]) & 0xffl;
829 final long l4 = ((long) record[offset + 4]) & 0xffl;
830 final long l3 = ((long) record[offset + 5]) & 0xffl;
831 final long l2 = ((long) record[offset + 6]) & 0xffl;
832 final long l1 = ((long) record[offset + 7]) & 0xffl;
833 final long l;
834 if (bigEndian) {
835 l = (l8 << 56) | (l7 << 48) | (l6 << 40) | (l5 << 32) |
836 (l4 << 24) | (l3 << 16) | (l2 << 8) | l1;
837 } else {
838 l = (l1 << 56) | (l2 << 48) | (l3 << 40) | (l4 << 32) |
839 (l5 << 24) | (l6 << 16) | (l7 << 8) | l8;
840 }
841 return Double.longBitsToDouble(l);
842 }
843
844
845
846
847
848
849 private int extractInt(final byte[] record, final int offset) {
850 final int l4 = ((int) record[offset + 0]) & 0xff;
851 final int l3 = ((int) record[offset + 1]) & 0xff;
852 final int l2 = ((int) record[offset + 2]) & 0xff;
853 final int l1 = ((int) record[offset + 3]) & 0xff;
854
855 if (bigEndian) {
856 return (l4 << 24) | (l3 << 16) | (l2 << 8) | l1;
857 } else {
858 return (l1 << 24) | (l2 << 16) | (l3 << 8) | l4;
859 }
860 }
861
862
863
864
865
866
867
868 private String extractString(final byte[] record, final int offset, final int length) {
869 return new String(record, offset, length, StandardCharsets.US_ASCII).trim();
870 }
871
872
873 private class ConstantsParser implements DataLoader {
874
875
876 private Map<String, Double> localConstants;
877
878
879
880
881 public Map<String, Double> getConstants() {
882 return localConstants;
883 }
884
885
886 public boolean stillAcceptsData() {
887 return localConstants == null;
888 }
889
890
891 public void loadData(final InputStream input, final String name)
892 throws IOException, ParseException, OrekitException {
893
894
895 final byte[] first = readFirstRecord(input, name);
896
897
898 final byte[] second = new byte[first.length];
899 if (!readInRecord(input, second, 0)) {
900 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
901 }
902
903 localConstants = parseConstants(first, second);
904
905 }
906
907 }
908
909
910 private class EphemerisParser implements DataLoader, TimeStampedGenerator<PosVelChebyshev> {
911
912
913 private final SortedSet<PosVelChebyshev> entries;
914
915
916 private AbsoluteDate start;
917
918
919 private AbsoluteDate end;
920
921
922
923 EphemerisParser() {
924 entries = new TreeSet<>(new ChronologicalComparator());
925 }
926
927
928 public List<PosVelChebyshev> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {
929 try {
930
931
932 entries.clear();
933 if (existingDate == null) {
934
935 start = date.shiftedBy(FIFTY_DAYS.negate());
936 end = date.shiftedBy(FIFTY_DAYS);
937 } else if (existingDate.compareTo(date) <= 0) {
938
939 start = existingDate;
940 end = date;
941 } else {
942
943 start = date;
944 end = existingDate;
945 }
946
947
948 if (!feed(this)) {
949 throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
950 }
951
952 return new ArrayList<>(entries);
953
954 } catch (OrekitException oe) {
955 throw new TimeStampedCacheException(oe);
956 }
957 }
958
959
960 public boolean stillAcceptsData() {
961
962
963 if (generateType == EphemerisType.EARTH) {
964 return false;
965 }
966
967
968
969 if (entries.isEmpty()) {
970 return true;
971 } else {
972
973 return !(entries.getFirst().getDate().compareTo(start) < 0 &&
974 entries.getLast().getDate().compareTo(end) > 0);
975 }
976
977 }
978
979
980 public void loadData(final InputStream input, final String name)
981 throws IOException {
982
983
984 final byte[] first = readFirstRecord(input, name);
985
986
987 final byte[] second = new byte[first.length];
988 if (!readInRecord(input, second, 0)) {
989 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
990 }
991
992 if (constants.get() == null) {
993 constants.compareAndSet(null, parseConstants(first, second));
994 }
995
996
997 final double au = 1000 * extractDouble(first, HEADER_ASTRONOMICAL_UNIT_OFFSET);
998 if (au < 1.4e11 || au > 1.6e11) {
999 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
1000 }
1001 if (FastMath.abs(getLoadedAstronomicalUnit() - au) >= 10.0) {
1002 throw new OrekitException(OrekitMessages.INCONSISTENT_ASTRONOMICAL_UNIT_IN_FILES,
1003 getLoadedAstronomicalUnit(), au);
1004 }
1005
1006
1007 final double emRat = extractDouble(first, HEADER_EM_RATIO_OFFSET);
1008 if (emRat < 80 || emRat > 82) {
1009 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
1010 }
1011 if (FastMath.abs(getLoadedEarthMoonMassRatio() - emRat) >= 1.0e-5) {
1012 throw new OrekitException(OrekitMessages.INCONSISTENT_EARTH_MOON_RATIO_IN_FILES,
1013 getLoadedEarthMoonMassRatio(), emRat);
1014 }
1015
1016
1017 parseFirstHeaderRecord(first, name);
1018
1019 if (startEpoch.compareTo(end) < 0 && finalEpoch.compareTo(start) > 0) {
1020
1021 final byte[] record = new byte[first.length];
1022 while (readInRecord(input, record, 0)) {
1023 final AbsoluteDate rangeStart = parseDataRecord(record);
1024 if (rangeStart.compareTo(end) > 0) {
1025
1026
1027 return;
1028 }
1029 }
1030 }
1031
1032 }
1033
1034
1035
1036
1037
1038 private AbsoluteDate parseDataRecord(final byte[] record) {
1039
1040
1041 final AbsoluteDate rangeStart = extractDate(record, DATA_START_RANGE_OFFSET);
1042 if (rangeStart.compareTo(startEpoch) < 0) {
1043 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE,
1044 rangeStart, startEpoch, finalEpoch, startEpoch.durationFrom(rangeStart));
1045 }
1046
1047 final AbsoluteDate rangeEnd = extractDate(record, DATE_END_RANGE_OFFSET);
1048 if (rangeEnd.compareTo(finalEpoch) > 0) {
1049 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_AFTER,
1050 rangeEnd, startEpoch, finalEpoch, rangeEnd.durationFrom(finalEpoch));
1051 }
1052
1053 if (rangeStart.compareTo(end) > 0 || rangeEnd.compareTo(start) < 0) {
1054
1055 return rangeEnd;
1056 }
1057
1058
1059 AbsoluteDate chunkEnd = rangeStart;
1060 final int nbChunks = chunks;
1061 final int nbCoeffs = coeffs;
1062 final int first = firstIndex;
1063 final TimeOffset duration = chunksDuration;
1064 for (int i = 0; i < nbChunks; ++i) {
1065
1066
1067 final AbsoluteDate chunkStart = chunkEnd;
1068 if (i == nbChunks - 1) {
1069 chunkEnd = rangeEnd;
1070 } else {
1071 chunkEnd = new AbsoluteDate(
1072 rangeStart,
1073 duration.multiply(i + 1),
1074 timeScale);
1075 }
1076
1077
1078
1079 final double[] xCoeffs = new double[nbCoeffs];
1080 final double[] yCoeffs = new double[nbCoeffs];
1081 final double[] zCoeffs = new double[nbCoeffs];
1082
1083 for (int k = 0; k < nbCoeffs; ++k) {
1084
1085
1086 final int index = first + components * i * nbCoeffs + k - 1;
1087 xCoeffs[k] = positionUnit * extractDouble(record, 8 * index);
1088 yCoeffs[k] = positionUnit * extractDouble(record, 8 * (index + nbCoeffs));
1089 zCoeffs[k] = positionUnit * extractDouble(record, 8 * (index + 2 * nbCoeffs));
1090 }
1091
1092
1093 entries.add(new PosVelChebyshev(chunkStart, timeScale,
1094 duration.toDouble(), xCoeffs, yCoeffs, zCoeffs));
1095
1096 }
1097
1098 return rangeStart;
1099
1100 }
1101
1102 }
1103
1104
1105 private class EphemerisRawPVProvider implements RawPVProvider {
1106
1107
1108
1109
1110
1111
1112 private PosVelChebyshev getChebyshev(final AbsoluteDate date) {
1113 PosVelChebyshev chebyshev;
1114 try {
1115 chebyshev = ephemerides.getNeighbors(date).findFirst().get();
1116 } catch (TimeStampedCacheException tce) {
1117
1118 chebyshev = ephemerides.getLatest();
1119 if (!chebyshev.inRange(date)) {
1120
1121 throw tce;
1122 }
1123 }
1124 return chebyshev;
1125 }
1126
1127
1128 public PVCoordinates getRawPV(final AbsoluteDate date) {
1129
1130
1131 final PosVelChebyshev chebyshev = getChebyshev(date);
1132
1133
1134 return chebyshev.getPositionVelocityAcceleration(date);
1135
1136 }
1137
1138
1139 public <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {
1140
1141
1142 final PosVelChebyshev chebyshev = getChebyshev(date.toAbsoluteDate());
1143
1144
1145 return chebyshev.getPositionVelocityAcceleration(date);
1146
1147 }
1148
1149
1150 @Override
1151 public Vector3D getRawPosition(final AbsoluteDate date) {
1152
1153 final PosVelChebyshev chebyshev = getChebyshev(date);
1154
1155
1156 return chebyshev.getPosition(date);
1157 }
1158
1159
1160 @Override
1161 public <T extends CalculusFieldElement<T>> FieldVector3D<T> getRawPosition(final FieldAbsoluteDate<T> date) {
1162
1163 final PosVelChebyshev chebyshev = getChebyshev(date.toAbsoluteDate());
1164
1165
1166 return chebyshev.getPosition(date);
1167 }
1168
1169 }
1170
1171
1172 static class ZeroRawPVProvider implements RawPVProvider {
1173
1174
1175 public PVCoordinates getRawPV(final AbsoluteDate date) {
1176 return PVCoordinates.ZERO;
1177 }
1178
1179
1180 public <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {
1181 return FieldPVCoordinates.getZero(date.getField());
1182 }
1183
1184 }
1185
1186 }