1   /* Copyright 2002-2020 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.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.atomic.AtomicReference;
30  
31  import org.hipparchus.RealFieldElement;
32  import org.hipparchus.util.FastMath;
33  import org.orekit.annotation.DefaultDataContext;
34  import org.orekit.data.AbstractSelfFeedingLoader;
35  import org.orekit.data.DataContext;
36  import org.orekit.data.DataLoader;
37  import org.orekit.data.DataProvidersManager;
38  import org.orekit.errors.OrekitException;
39  import org.orekit.errors.OrekitInternalError;
40  import org.orekit.errors.OrekitMessages;
41  import org.orekit.errors.TimeStampedCacheException;
42  import org.orekit.frames.Frame;
43  import org.orekit.frames.Predefined;
44  import org.orekit.time.AbsoluteDate;
45  import org.orekit.time.ChronologicalComparator;
46  import org.orekit.time.DateComponents;
47  import org.orekit.time.FieldAbsoluteDate;
48  import org.orekit.time.TimeComponents;
49  import org.orekit.time.TimeScale;
50  import org.orekit.time.TimeScales;
51  import org.orekit.utils.Constants;
52  import org.orekit.utils.FieldPVCoordinates;
53  import org.orekit.utils.OrekitConfiguration;
54  import org.orekit.utils.PVCoordinates;
55  import org.orekit.utils.GenericTimeStampedCache;
56  import org.orekit.utils.TimeStampedGenerator;
57  
58  /** Loader for JPL ephemerides binary files (DE 4xx) and similar formats (INPOP 06/08/10).
59   * <p>JPL ephemerides binary files contain ephemerides for all solar system planets.</p>
60   * <p>The JPL ephemerides binary files are recognized thanks to their base names,
61   * which must match the pattern <code>[lu]nx[mp]####.ddd</code> (or
62   * <code>[lu]nx[mp]####.ddd.gz</code> for gzip-compressed files) where # stands for a
63   * digit character and where ddd is an ephemeris type (typically 405 or 406).</p>
64   * <p>The loader supports files encoded in big-endian as well as in little-endian notation.
65   * Usually, big-endian files are named <code>unx[mp]####.ddd</code>, while little-endian files
66   * are named <code>lnx[mp]####.ddd</code>.</p>
67   * <p>The IMCCE ephemerides binary files are recognized thanks to their base names,
68   * which must match the pattern <code>inpop*.dat</code> (or
69   * <code>inpop*.dat.gz</code> for gzip-compressed files) where * stands for any string.</p>
70   * <p>The loader supports files encoded in big-endian as well as in little-endian notation.
71   * Usually, big-endian files contain <code>bigendian</code> in their names, while little-endian files
72   * contain <code>littleendian</code> in their names.</p>
73   * <p>The loader supports files in TDB or TCB time scales.</p>
74   * @author Luc Maisonobe
75   */
76  public class JPLEphemeridesLoader extends AbstractSelfFeedingLoader
77          implements CelestialBodyLoader {
78  
79      /** Default supported files name pattern for JPL DE files. */
80      public static final String DEFAULT_DE_SUPPORTED_NAMES = "^[lu]nx([mp](\\d\\d\\d\\d))+\\.(?:4\\d\\d)$";
81  
82      /** Default supported files name pattern for IMCCE INPOP files. */
83      public static final String DEFAULT_INPOP_SUPPORTED_NAMES = "^inpop.*\\.dat$";
84  
85      /** 50 days in seconds. */
86      private static final double FIFTY_DAYS = 50 * Constants.JULIAN_DAY;
87  
88      /** DE number used by INPOP files. */
89      private static final int INPOP_DE_NUMBER = 100;
90  
91      /** Maximal number of constants in headers. */
92      private static final int CONSTANTS_MAX_NUMBER           = 400;
93  
94      /** Offset of the ephemeris type in first header record. */
95      private static final int HEADER_EPHEMERIS_TYPE_OFFSET   = 2840;
96  
97      /** Offset of the record size (for INPOP files) in first header record. */
98      private static final int HEADER_RECORD_SIZE_OFFSET      = 2856;
99  
100     /** Offset of the start epoch in first header record. */
101     private static final int HEADER_START_EPOCH_OFFSET      = 2652;
102 
103     /** Offset of the end epoch in first header record. */
104     private static final int HEADER_END_EPOCH_OFFSET        = 2660;
105 
106     /** Offset of the astronomical unit in first header record. */
107     private static final int HEADER_ASTRONOMICAL_UNIT_OFFSET = 2680;
108 
109     /** Offset of the Earth-Moon mass ratio in first header record. */
110     private static final int HEADER_EM_RATIO_OFFSET         = 2688;
111 
112     /** Offset of Chebishev coefficients indices in first header record. */
113     private static final int HEADER_CHEBISHEV_INDICES_OFFSET = 2696;
114 
115     /** Offset of libration coefficients indices in first header record. */
116     private static final int HEADER_LIBRATION_INDICES_OFFSET = 2844;
117 
118     /** Offset of chunks duration in first header record. */
119     private static final int HEADER_CHUNK_DURATION_OFFSET    = 2668;
120 
121     /** Offset of the constants names in first header record. */
122     private static final int HEADER_CONSTANTS_NAMES_OFFSET  = 252;
123 
124     /** Offset of the constants values in second header record. */
125     private static final int HEADER_CONSTANTS_VALUES_OFFSET = 0;
126 
127     /** Offset of the range start in the data records. */
128     private static final int DATA_START_RANGE_OFFSET        = 0;
129 
130     /** Offset of the range end in the data records. */
131     private static final int DATE_END_RANGE_OFFSET          = 8;
132 
133     /** The constant name for the astronomical unit. */
134     private static final String CONSTANT_AU = "AU";
135 
136     /** The constant name for the earth-moon mass ratio. */
137     private static final String CONSTANT_EMRAT = "EMRAT";
138 
139     /** List of supported ephemerides types. */
140     public enum EphemerisType {
141 
142         /** Constant for solar system barycenter. */
143         SOLAR_SYSTEM_BARYCENTER,
144 
145         /** Constant for the Sun. */
146         SUN,
147 
148         /** Constant for Mercury. */
149         MERCURY,
150 
151         /** Constant for Venus. */
152         VENUS,
153 
154         /** Constant for the Earth-Moon barycenter. */
155         EARTH_MOON,
156 
157         /** Constant for the Earth. */
158         EARTH,
159 
160         /** Constant for the Moon. */
161         MOON,
162 
163         /** Constant for Mars. */
164         MARS,
165 
166         /** Constant for Jupiter. */
167         JUPITER,
168 
169         /** Constant for Saturn. */
170         SATURN,
171 
172         /** Constant for Uranus. */
173         URANUS,
174 
175         /** Constant for Neptune. */
176         NEPTUNE,
177 
178         /** Constant for Pluto. */
179         PLUTO
180 
181     }
182 
183     /** Interface for raw position-velocity retrieval. */
184     public interface RawPVProvider {
185 
186         /** Get the position-velocity at date.
187          * @param date date at which the position-velocity is desired
188          * @return position-velocity at the specified date
189          */
190         PVCoordinates getRawPV(AbsoluteDate date);
191 
192         /** Get the position-velocity at date.
193          * @param date date at which the position-velocity is desired
194          * @param <T> type of the field elements
195          * @return position-velocity at the specified date
196          */
197         <T extends RealFieldElement<T>> FieldPVCoordinates<T> getRawPV(FieldAbsoluteDate<T> date);
198 
199     }
200 
201     /** Ephemeris for selected body. */
202     private final GenericTimeStampedCache<PosVelChebyshev> ephemerides;
203 
204     /** Constants defined in the file. */
205     private final AtomicReference<Map<String, Double>> constants;
206 
207     /** Ephemeris type to generate. */
208     private final EphemerisType generateType;
209 
210     /** Ephemeris type to load. */
211     private final EphemerisType loadType;
212 
213     /** Time scales to use when loading data. */
214     private final TimeScales timeScales;
215 
216     /** The GCRF implementation. */
217     private final Frame gcrf;
218 
219     /** Current file start epoch. */
220     private AbsoluteDate startEpoch;
221 
222     /** Current file final epoch. */
223     private AbsoluteDate finalEpoch;
224 
225     /** Chunks duration (in seconds). */
226     private double maxChunksDuration;
227 
228     /** Current file chunks duration (in seconds). */
229     private double chunksDuration;
230 
231     /** Index of the first data for selected body. */
232     private int firstIndex;
233 
234     /** Number of coefficients for selected body. */
235     private int coeffs;
236 
237     /** Number of chunks for the selected body. */
238     private int chunks;
239 
240     /** Number of components contained in the file. */
241     private int components;
242 
243     /** Unit of the position coordinates (as a multiple of meters). */
244     private double positionUnit;
245 
246     /** Time scale of the date coordinates. */
247     private TimeScale timeScale;
248 
249     /** Indicator for binary file endianness. */
250     private boolean bigEndian;
251 
252     /** Create a loader for JPL ephemerides binary files. This constructor uses the {@link
253      * DataContext#getDefault() default data context}.
254      *
255      * @param supportedNames regular expression for supported files names
256      * @param generateType ephemeris type to generate
257      * @see #JPLEphemeridesLoader(String, EphemerisType, DataProvidersManager, TimeScales,
258      * Frame)
259      */
260     @DefaultDataContext
261     public JPLEphemeridesLoader(final String supportedNames, final EphemerisType generateType) {
262         this(supportedNames, generateType,
263                 DataContext.getDefault().getDataProvidersManager(),
264                 DataContext.getDefault().getTimeScales(),
265                 DataContext.getDefault().getFrames().getGCRF());
266     }
267 
268     /** Create a loader for JPL ephemerides binary files.
269      * @param supportedNames regular expression for supported files names
270      * @param generateType ephemeris type to generate
271      * @param dataProvidersManager provides access to the ephemeris files.
272      * @param timeScales used to access the TCB and TDB time scales while loading data.
273      * @param gcrf Earth centered frame aligned with ICRF.
274      * @since 10.1
275      */
276     public JPLEphemeridesLoader(final String supportedNames,
277                                 final EphemerisType generateType,
278                                 final DataProvidersManager dataProvidersManager,
279                                 final TimeScales timeScales,
280                                 final Frame gcrf) {
281         super(supportedNames, dataProvidersManager);
282 
283         this.timeScales = timeScales;
284         this.gcrf = gcrf;
285         constants = new AtomicReference<>();
286 
287         this.generateType  = generateType;
288         if (generateType == EphemerisType.SOLAR_SYSTEM_BARYCENTER) {
289             loadType = EphemerisType.EARTH_MOON;
290         } else if (generateType == EphemerisType.EARTH_MOON) {
291             loadType = EphemerisType.MOON;
292         } else {
293             loadType = generateType;
294         }
295 
296         ephemerides = new GenericTimeStampedCache<>(
297                 2, OrekitConfiguration.getCacheSlotsNumber(),
298                 Double.POSITIVE_INFINITY, FIFTY_DAYS,
299                 new EphemerisParser());
300         maxChunksDuration = Double.NaN;
301         chunksDuration    = Double.NaN;
302 
303     }
304 
305     /** Load celestial body.
306      * @param name name of the celestial body
307      * @return loaded celestial body
308      */
309     public CelestialBody loadCelestialBody(final String name) {
310 
311         final double gm       = getLoadedGravitationalCoefficient(generateType);
312         final IAUPole iauPole = PredefinedIAUPoles
313                 .getIAUPole(generateType, timeScales);
314         final double scale;
315         final Frame definingFrameAlignedWithICRF;
316         final RawPVProvider rawPVProvider;
317         String inertialFrameName = null;
318         String bodyOrientedFrameName = null;
319         switch (generateType) {
320             case SOLAR_SYSTEM_BARYCENTER : {
321                 scale = -1.0;
322                 final JPLEphemeridesLoadermeridesLoader">JPLEphemeridesLoader parentLoader = new JPLEphemeridesLoader(
323                         getSupportedNames(),
324                         EphemerisType.EARTH_MOON,
325                         getDataProvidersManager(),
326                         timeScales,
327                         gcrf);
328                 final CelestialBody parentBody =
329                         parentLoader.loadCelestialBody(CelestialBodyFactory.EARTH_MOON);
330                 definingFrameAlignedWithICRF = parentBody.getInertiallyOrientedFrame();
331                 rawPVProvider = new EphemerisRawPVProvider();
332                 inertialFrameName = Predefined.ICRF.getName();
333                 bodyOrientedFrameName = null;
334                 break;
335             }
336             case EARTH_MOON :
337                 scale         = 1.0 / (1.0 + getLoadedEarthMoonMassRatio());
338                 definingFrameAlignedWithICRF = gcrf;
339                 rawPVProvider = new EphemerisRawPVProvider();
340                 break;
341             case EARTH :
342                 scale         = 1.0;
343                 definingFrameAlignedWithICRF = gcrf;
344                 rawPVProvider = new ZeroRawPVProvider();
345                 break;
346             case MOON :
347                 scale         =  1.0;
348                 definingFrameAlignedWithICRF = gcrf;
349                 rawPVProvider = new EphemerisRawPVProvider();
350                 break;
351             default : {
352                 scale = 1.0;
353                 final JPLEphemeridesLoadermeridesLoader">JPLEphemeridesLoader parentLoader = new JPLEphemeridesLoader(
354                         getSupportedNames(),
355                         EphemerisType.SOLAR_SYSTEM_BARYCENTER,
356                         getDataProvidersManager(),
357                         timeScales,
358                         gcrf);
359                 final CelestialBody parentBody =
360                         parentLoader.loadCelestialBody(CelestialBodyFactory.SOLAR_SYSTEM_BARYCENTER);
361                 definingFrameAlignedWithICRF = parentBody.getInertiallyOrientedFrame();
362                 rawPVProvider = new EphemerisRawPVProvider();
363             }
364         }
365 
366         // build the celestial body
367         return new JPLCelestialBody(name, getSupportedNames(), generateType, rawPVProvider,
368                                     gm, scale, iauPole, definingFrameAlignedWithICRF,
369                                     inertialFrameName, bodyOrientedFrameName);
370 
371     }
372 
373     /** Get astronomical unit.
374      * @return astronomical unit in meters
375      */
376     public double getLoadedAstronomicalUnit() {
377         return 1000.0 * getLoadedConstant(CONSTANT_AU);
378     }
379 
380     /** Get Earth/Moon mass ratio.
381      * @return Earth/Moon mass ratio
382      */
383     public double getLoadedEarthMoonMassRatio() {
384         return getLoadedConstant(CONSTANT_EMRAT);
385     }
386 
387     /** Get the gravitational coefficient of a body.
388      * @param body body for which the gravitational coefficient is requested
389      * @return gravitational coefficient in m³/s²
390      */
391     public double getLoadedGravitationalCoefficient(final EphemerisType body) {
392 
393         // coefficient in au³/day²
394         final double rawGM;
395         switch (body) {
396             case SOLAR_SYSTEM_BARYCENTER :
397                 return getLoadedGravitationalCoefficient(EphemerisType.SUN)        +
398                         getLoadedGravitationalCoefficient(EphemerisType.MERCURY)    +
399                         getLoadedGravitationalCoefficient(EphemerisType.VENUS)      +
400                         getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) +
401                         getLoadedGravitationalCoefficient(EphemerisType.MARS)       +
402                         getLoadedGravitationalCoefficient(EphemerisType.JUPITER)    +
403                         getLoadedGravitationalCoefficient(EphemerisType.SATURN)     +
404                         getLoadedGravitationalCoefficient(EphemerisType.URANUS)     +
405                         getLoadedGravitationalCoefficient(EphemerisType.NEPTUNE)    +
406                         getLoadedGravitationalCoefficient(EphemerisType.PLUTO);
407             case SUN :
408                 rawGM = getLoadedConstant("GMS", "GM_Sun");
409                 break;
410             case MERCURY :
411                 rawGM = getLoadedConstant("GM1", "GM_Mer");
412                 break;
413             case VENUS :
414                 rawGM = getLoadedConstant("GM2", "GM_Ven");
415                 break;
416             case EARTH_MOON :
417                 rawGM = getLoadedConstant("GMB", "GM_EMB");
418                 break;
419             case EARTH :
420                 return getLoadedEarthMoonMassRatio() *
421                         getLoadedGravitationalCoefficient(EphemerisType.MOON);
422             case MOON :
423                 return getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) /
424                         (1.0 + getLoadedEarthMoonMassRatio());
425             case MARS :
426                 rawGM = getLoadedConstant("GM4", "GM_Mar");
427                 break;
428             case JUPITER :
429                 rawGM = getLoadedConstant("GM5", "GM_Jup");
430                 break;
431             case SATURN :
432                 rawGM = getLoadedConstant("GM6", "GM_Sat");
433                 break;
434             case URANUS :
435                 rawGM = getLoadedConstant("GM7", "GM_Ura");
436                 break;
437             case NEPTUNE :
438                 rawGM = getLoadedConstant("GM8", "GM_Nep");
439                 break;
440             case PLUTO :
441                 rawGM = getLoadedConstant("GM9", "GM_Plu");
442                 break;
443             default :
444                 throw new OrekitInternalError(null);
445         }
446 
447         final double au    = getLoadedAstronomicalUnit();
448         return rawGM * au * au * au / (Constants.JULIAN_DAY * Constants.JULIAN_DAY);
449 
450     }
451 
452     /** Get a constant defined in the ephemerides headers.
453      * <p>Note that since constants are defined in the JPL headers
454      * files, they are available as soon as one file is available, even
455      * if it doesn't match the desired central date. This is because the
456      * header must be parsed before the dates can be checked.</p>
457      * <p>
458      * There are alternate names for constants since for example JPL names are
459      * different from INPOP names (Sun gravity: GMS or GM_Sun, Mars gravity:
460      * GM4 or GM_Mar...).
461      * </p>
462      * @param names alternate names of the constant
463      * @return value of the constant of NaN if the constant is not defined
464      */
465     public double getLoadedConstant(final String... names) {
466 
467         // lazy loading of constants
468         Map<String, Double> map = constants.get();
469         if (map == null) {
470             final ConstantsParser parser = new ConstantsParser();
471             if (!feed(parser)) {
472                 throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
473             }
474             map = parser.getConstants();
475             constants.compareAndSet(null, map);
476         }
477 
478         for (final String name : names) {
479             if (map.containsKey(name)) {
480                 return map.get(name);
481             }
482         }
483 
484         return Double.NaN;
485 
486     }
487 
488     /** Get the maximal chunks duration.
489      * @return chunks maximal duration in seconds
490      */
491     public double getMaxChunksDuration() {
492         return maxChunksDuration;
493     }
494 
495     /** Parse the first header record.
496      * @param record first header record
497      * @param name name of the file (or zip entry)
498      */
499     private void parseFirstHeaderRecord(final byte[] record, final String name) {
500 
501         // get the ephemerides type
502         final int deNum = extractInt(record, HEADER_EPHEMERIS_TYPE_OFFSET);
503 
504         // as default, 3 polynomial coefficients for the Cartesian coordinates
505         // (x, y, z) are contained in the file, positions are in kilometers
506         // and times are in TDB
507         components   = 3;
508         positionUnit = 1000.0;
509         timeScale    = timeScales.getTDB();
510 
511         if (deNum == INPOP_DE_NUMBER) {
512             // an INPOP file may contain 6 components (including coefficients for the velocity vector)
513             final double format = getLoadedConstant("FORMAT");
514             if (!Double.isNaN(format) && (int) FastMath.IEEEremainder(format, 10) != 1) {
515                 components = 6;
516             }
517 
518             // INPOP files may have their polynomials expressed in AU
519             final double unite = getLoadedConstant("UNITE");
520             if (!Double.isNaN(unite) && (int) unite == 0) {
521                 positionUnit = getLoadedAstronomicalUnit();
522             }
523 
524             // INPOP files may have their times expressed in TCB
525             final double timesc = getLoadedConstant("TIMESC");
526             if (!Double.isNaN(timesc) && (int) timesc == 1) {
527                 timeScale = timeScales.getTCB();
528             }
529 
530         }
531 
532         // extract covered date range
533         startEpoch = extractDate(record, HEADER_START_EPOCH_OFFSET);
534         finalEpoch = extractDate(record, HEADER_END_EPOCH_OFFSET);
535         boolean ok = finalEpoch.compareTo(startEpoch) > 0;
536 
537         // indices of the Chebyshev coefficients for each ephemeris
538         for (int i = 0; i < 12; ++i) {
539             final int row1 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET     + 12 * i);
540             final int row2 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 4 + 12 * i);
541             final int row3 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 8 + 12 * i);
542             ok = ok && (row1 >= 0) && (row2 >= 0) && (row3 >= 0);
543             if (((i ==  0) && (loadType == EphemerisType.MERCURY))    ||
544                     ((i ==  1) && (loadType == EphemerisType.VENUS))      ||
545                     ((i ==  2) && (loadType == EphemerisType.EARTH_MOON)) ||
546                     ((i ==  3) && (loadType == EphemerisType.MARS))       ||
547                     ((i ==  4) && (loadType == EphemerisType.JUPITER))    ||
548                     ((i ==  5) && (loadType == EphemerisType.SATURN))     ||
549                     ((i ==  6) && (loadType == EphemerisType.URANUS))     ||
550                     ((i ==  7) && (loadType == EphemerisType.NEPTUNE))    ||
551                     ((i ==  8) && (loadType == EphemerisType.PLUTO))      ||
552                     ((i ==  9) && (loadType == EphemerisType.MOON))       ||
553                     ((i == 10) && (loadType == EphemerisType.SUN))) {
554                 firstIndex = row1;
555                 coeffs     = row2;
556                 chunks     = row3;
557             }
558         }
559 
560         // compute chunks duration
561         final double timeSpan = extractDouble(record, HEADER_CHUNK_DURATION_OFFSET);
562         ok = ok && (timeSpan > 0) && (timeSpan < 100);
563         chunksDuration = Constants.JULIAN_DAY * (timeSpan / chunks);
564         if (Double.isNaN(maxChunksDuration)) {
565             maxChunksDuration = chunksDuration;
566         } else {
567             maxChunksDuration = FastMath.max(maxChunksDuration, chunksDuration);
568         }
569 
570         // sanity checks
571         if (!ok) {
572             throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
573         }
574 
575     }
576 
577     /** Read first header record.
578      * @param input input stream
579      * @param name name of the file (or zip entry)
580      * @return record record where to put bytes
581      * @exception IOException if a read error occurs
582      */
583     private byte[] readFirstRecord(final InputStream input, final String name)
584         throws IOException {
585 
586         // read first part of record, up to the ephemeris type
587         final byte[] firstPart = new byte[HEADER_RECORD_SIZE_OFFSET + 4];
588         if (!readInRecord(input, firstPart, 0)) {
589             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
590         }
591 
592         // detect the endian format
593         detectEndianess(firstPart);
594 
595         // get the ephemerides type
596         final int deNum = extractInt(firstPart, HEADER_EPHEMERIS_TYPE_OFFSET);
597 
598         // the record size for this file
599         final int recordSize;
600 
601         if (deNum == INPOP_DE_NUMBER) {
602             // INPOP files have an extended DE format, which includes also the record size
603             recordSize = extractInt(firstPart, HEADER_RECORD_SIZE_OFFSET) << 3;
604         } else {
605             // compute the record size for original JPL files
606             recordSize = computeRecordSize(firstPart, name);
607         }
608 
609         if (recordSize <= 0) {
610             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
611         }
612 
613         // build a record with the proper size and finish read of the first complete record
614         final int start = firstPart.length;
615         final byte[] record = new byte[recordSize];
616         System.arraycopy(firstPart, 0, record, 0, firstPart.length);
617         if (!readInRecord(input, record, start)) {
618             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
619         }
620 
621         return record;
622 
623     }
624 
625     /** Parse constants from first two header records.
626      * @param first first header record
627      * @param second second header record
628      * @return map of parsed constants
629      */
630     private Map<String, Double> parseConstants(final byte[] first, final byte[] second) {
631 
632         final Map<String, Double> map = new HashMap<>();
633 
634         for (int i = 0; i < CONSTANTS_MAX_NUMBER; ++i) {
635             // Note: for extracting the strings from the binary file, it makes no difference
636             //       if the file is stored in big-endian or little-endian notation
637             final String constantName = extractString(first, HEADER_CONSTANTS_NAMES_OFFSET + i * 6, 6);
638             if (constantName.length() == 0) {
639                 // no more constants to read
640                 break;
641             }
642             final double constantValue = extractDouble(second, HEADER_CONSTANTS_VALUES_OFFSET + 8 * i);
643             map.put(constantName, constantValue);
644         }
645 
646         // INPOP files do not have constants for AU and EMRAT, thus extract them from
647         // the header record and create a constant for them to be consistent with JPL files
648         if (!map.containsKey(CONSTANT_AU)) {
649             map.put(CONSTANT_AU, extractDouble(first, HEADER_ASTRONOMICAL_UNIT_OFFSET));
650         }
651 
652         if (!map.containsKey(CONSTANT_EMRAT)) {
653             map.put(CONSTANT_EMRAT, extractDouble(first, HEADER_EM_RATIO_OFFSET));
654         }
655 
656         return map;
657 
658     }
659 
660     /** Read bytes into the current record array.
661      * @param input input stream
662      * @param record record where to put bytes
663      * @param start start index where to put bytes
664      * @return true if record has been filled up
665      * @exception IOException if a read error occurs
666      */
667     private boolean readInRecord(final InputStream input, final byte[] record, final int start)
668         throws IOException {
669         int index = start;
670         while (index != record.length) {
671             final int n = input.read(record, index, record.length - index);
672             if (n < 0) {
673                 return false;
674             }
675             index += n;
676         }
677         return true;
678     }
679 
680     /** Detect whether the JPL ephemerides file is stored in big-endian or
681      * little-endian notation.
682      * @param record the array containing the binary JPL header
683      */
684     private void detectEndianess(final byte[] record) {
685 
686         // default to big-endian
687         bigEndian = true;
688 
689         // first try to read the DE number in big-endian format
690         // the number is stored as unsigned int, so we have to convert it properly
691         final long deNum = extractInt(record, HEADER_EPHEMERIS_TYPE_OFFSET) & 0xffffffffL;
692 
693         // simple heuristic: if the read value is larger than half the range of an integer
694         //                   assume the file is in little-endian format
695         if (deNum > (1 << 15)) {
696             bigEndian = false;
697         }
698 
699     }
700 
701     /** Calculate the record size of a JPL ephemerides file.
702      * @param record the byte array containing the header record
703      * @param name the name of the data file
704      * @return the record size for this file
705      */
706     private int computeRecordSize(final byte[] record, final String name) {
707 
708         int recordSize = 0;
709         boolean ok = true;
710         // JPL files always have 3 position components
711         final int nComp = 3;
712 
713         // iterate over the coefficient ptr array and sum up the record size
714         // the coeffPtr array has the dimensions [12][nComp]
715         for (int j = 0; j < 12; j++) {
716             final int nCompCur = (j == 11) ? 2 : nComp;
717 
718             // Note: the array element coeffPtr[j][0] is not needed for the calculation
719             final int idx = HEADER_CHEBISHEV_INDICES_OFFSET + j * nComp * 4;
720             final int coeffPtr1 = extractInt(record, idx + 4);
721             final int coeffPtr2 = extractInt(record, idx + 8);
722 
723             // sanity checks
724             ok = ok && (coeffPtr1 >= 0 || coeffPtr2 >= 0);
725 
726             recordSize += coeffPtr1 * coeffPtr2 * nCompCur;
727         }
728 
729         // the libration ptr array has the dimension [3]
730         // Note: the array element libratPtr[0] is not needed for the calculation
731         final int libratPtr1 = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 4);
732         final int libratPtr2 = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 8);
733 
734         // sanity checks
735         ok = ok && (libratPtr1 >= 0 || libratPtr2 >= 0);
736 
737         recordSize += libratPtr1 * libratPtr2 * nComp + 2;
738         recordSize <<= 3;
739 
740         if (!ok || recordSize <= 0) {
741             throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
742         }
743 
744         return recordSize;
745 
746     }
747 
748     /** Extract a date from a record.
749      * @param record record to parse
750      * @param offset offset of the double within the record
751      * @return extracted date
752      */
753     private AbsoluteDate extractDate(final byte[] record, final int offset) {
754 
755         final double t = extractDouble(record, offset);
756         int    jDay    = (int) FastMath.floor(t);
757         double seconds = (t + 0.5 - jDay) * Constants.JULIAN_DAY;
758         if (seconds >= Constants.JULIAN_DAY) {
759             ++jDay;
760             seconds -= Constants.JULIAN_DAY;
761         }
762         return new AbsoluteDate(new DateComponents(DateComponents.JULIAN_EPOCH, jDay),
763                                 new TimeComponents(seconds), timeScale);
764     }
765 
766     /** Extract a double from a record.
767      * <p>Double numbers are stored according to IEEE 754 standard, with
768      * most significant byte first.</p>
769      * @param record record to parse
770      * @param offset offset of the double within the record
771      * @return extracted double
772      */
773     private double extractDouble(final byte[] record, final int offset) {
774         final long l8 = ((long) record[offset + 0]) & 0xffl;
775         final long l7 = ((long) record[offset + 1]) & 0xffl;
776         final long l6 = ((long) record[offset + 2]) & 0xffl;
777         final long l5 = ((long) record[offset + 3]) & 0xffl;
778         final long l4 = ((long) record[offset + 4]) & 0xffl;
779         final long l3 = ((long) record[offset + 5]) & 0xffl;
780         final long l2 = ((long) record[offset + 6]) & 0xffl;
781         final long l1 = ((long) record[offset + 7]) & 0xffl;
782         final long l;
783         if (bigEndian) {
784             l = (l8 << 56) | (l7 << 48) | (l6 << 40) | (l5 << 32) |
785                 (l4 << 24) | (l3 << 16) | (l2 <<  8) | l1;
786         } else {
787             l = (l1 << 56) | (l2 << 48) | (l3 << 40) | (l4 << 32) |
788                 (l5 << 24) | (l6 << 16) | (l7 <<  8) | l8;
789         }
790         return Double.longBitsToDouble(l);
791     }
792 
793     /** Extract an int from a record.
794      * @param record record to parse
795      * @param offset offset of the double within the record
796      * @return extracted int
797      */
798     private int extractInt(final byte[] record, final int offset) {
799         final int l4 = ((int) record[offset + 0]) & 0xff;
800         final int l3 = ((int) record[offset + 1]) & 0xff;
801         final int l2 = ((int) record[offset + 2]) & 0xff;
802         final int l1 = ((int) record[offset + 3]) & 0xff;
803 
804         if (bigEndian) {
805             return (l4 << 24) | (l3 << 16) | (l2 <<  8) | l1;
806         } else {
807             return (l1 << 24) | (l2 << 16) | (l3 <<  8) | l4;
808         }
809     }
810 
811     /** Extract a String from a record.
812      * @param record record to parse
813      * @param offset offset of the string within the record
814      * @param length maximal length of the string
815      * @return extracted string, with whitespace characters stripped
816      */
817     private String extractString(final byte[] record, final int offset, final int length) {
818         return new String(record, offset, length, StandardCharsets.US_ASCII).trim();
819     }
820 
821     /** Local parser for header constants. */
822     private class ConstantsParser implements DataLoader {
823 
824         /** Local constants map. */
825         private Map<String, Double> localConstants;
826 
827        /** Get the local constants map.
828          * @return local constants map
829          */
830         public Map<String, Double> getConstants() {
831             return localConstants;
832         }
833 
834         /** {@inheritDoc} */
835         public boolean stillAcceptsData() {
836             return localConstants == null;
837         }
838 
839         /** {@inheritDoc} */
840         public void loadData(final InputStream input, final String name)
841             throws IOException, ParseException, OrekitException {
842 
843             // read first header record
844             final byte[] first = readFirstRecord(input, name);
845 
846             // the second record contains the values of the constants used for least-square filtering
847             final byte[] second = new byte[first.length];
848             if (!readInRecord(input, second, 0)) {
849                 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
850             }
851 
852             localConstants = parseConstants(first, second);
853 
854         }
855 
856     }
857 
858     /** Local parser for Chebyshev polynomials. */
859     private class EphemerisParser implements DataLoader, TimeStampedGenerator<PosVelChebyshev> {
860 
861         /** Set of Chebyshev polynomials read. */
862         private final SortedSet<PosVelChebyshev> entries;
863 
864         /** Start of range we are interested in. */
865         private AbsoluteDate start;
866 
867         /** End of range we are interested in. */
868         private AbsoluteDate end;
869 
870         /** Simple constructor.
871          */
872         EphemerisParser() {
873             entries = new TreeSet<>(new ChronologicalComparator());
874         }
875 
876         /** {@inheritDoc} */
877         public List<PosVelChebyshev> generate(final AbsoluteDatete">AbsoluteDate existingDate, final AbsoluteDate date) {
878             try {
879 
880                 // prepare reading
881                 entries.clear();
882                 if (existingDate == null) {
883                     // we want ephemeris data for the first time, set up an arbitrary first range
884                     start = date.shiftedBy(-FIFTY_DAYS);
885                     end   = date.shiftedBy(+FIFTY_DAYS);
886                 } else if (existingDate.compareTo(date) <= 0) {
887                     // we want to extend an existing range towards future dates
888                     start = existingDate;
889                     end   = date;
890                 } else {
891                     // we want to extend an existing range towards past dates
892                     start = date;
893                     end   = existingDate;
894                 }
895 
896                 // get new entries in the specified data range
897                 if (!feed(this)) {
898                     throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
899                 }
900 
901                 return new ArrayList<>(entries);
902 
903             } catch (OrekitException oe) {
904                 throw new TimeStampedCacheException(oe);
905             }
906         }
907 
908         /** {@inheritDoc} */
909         public boolean stillAcceptsData() {
910 
911             // special case for Earth: we do not really load any ephemeris data
912             if (generateType == EphemerisType.EARTH) {
913                 return false;
914             }
915 
916             // we have to look for data in all available ephemerides files as there may be
917             // data overlaps that result in incomplete data
918             if (entries.isEmpty()) {
919                 return true;
920             } else {
921                 // if the requested range is already filled, we do not need to look further
922                 return !(entries.first().getDate().compareTo(start) < 0 &&
923                          entries.last().getDate().compareTo(end)    > 0);
924             }
925 
926         }
927 
928         /** {@inheritDoc} */
929         public void loadData(final InputStream input, final String name)
930             throws IOException {
931 
932             // read first header record
933             final byte[] first = readFirstRecord(input, name);
934 
935             // the second record contains the values of the constants used for least-square filtering
936             final byte[] second = new byte[first.length];
937             if (!readInRecord(input, second, 0)) {
938                 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
939             }
940 
941             if (constants.get() == null) {
942                 constants.compareAndSet(null, parseConstants(first, second));
943             }
944 
945             // check astronomical unit consistency
946             final double au = 1000 * extractDouble(first, HEADER_ASTRONOMICAL_UNIT_OFFSET);
947             if ((au < 1.4e11) || (au > 1.6e11)) {
948                 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
949             }
950             if (FastMath.abs(getLoadedAstronomicalUnit() - au) >= 10.0) {
951                 throw new OrekitException(OrekitMessages.INCONSISTENT_ASTRONOMICAL_UNIT_IN_FILES,
952                                           getLoadedAstronomicalUnit(), au);
953             }
954 
955             // check Earth-Moon mass ratio consistency
956             final double emRat = extractDouble(first, HEADER_EM_RATIO_OFFSET);
957             if ((emRat < 80) || (emRat > 82)) {
958                 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
959             }
960             if (FastMath.abs(getLoadedEarthMoonMassRatio() - emRat) >= 1.0e-5) {
961                 throw new OrekitException(OrekitMessages.INCONSISTENT_EARTH_MOON_RATIO_IN_FILES,
962                                           getLoadedEarthMoonMassRatio(), emRat);
963             }
964 
965             // parse first header record
966             parseFirstHeaderRecord(first, name);
967 
968             if (startEpoch.compareTo(end) < 0 && finalEpoch.compareTo(start) > 0) {
969                 // this file contains data in the range we are looking for, read it
970                 final byte[] record = new byte[first.length];
971                 while (readInRecord(input, record, 0)) {
972                     final AbsoluteDate rangeStart = parseDataRecord(record);
973                     if (rangeStart.compareTo(end) > 0) {
974                         // we have already exceeded the range we were interested in,
975                         // we interrupt parsing here
976                         return;
977                     }
978                 }
979             }
980 
981         }
982 
983         /** Parse regular ephemeris record.
984          * @param record record to parse
985          * @return date of the last parsed chunk
986          */
987         private AbsoluteDate parseDataRecord(final byte[] record) {
988 
989             // extract time range covered by the record
990             final AbsoluteDate rangeStart = extractDate(record, DATA_START_RANGE_OFFSET);
991             if (rangeStart.compareTo(startEpoch) < 0) {
992                 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE, rangeStart, startEpoch, finalEpoch);
993             }
994 
995             final AbsoluteDate rangeEnd   = extractDate(record, DATE_END_RANGE_OFFSET);
996             if (rangeEnd.compareTo(finalEpoch) > 0) {
997                 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE, rangeEnd, startEpoch, finalEpoch);
998             }
999 
1000             if (rangeStart.compareTo(end) > 0 || rangeEnd.compareTo(start) < 0) {
1001                 // we are not interested in this record, don't parse it
1002                 return rangeEnd;
1003             }
1004 
1005             // loop over chunks inside the time range
1006             AbsoluteDate chunkEnd = rangeStart;
1007             final int nbChunks    = chunks;
1008             final int nbCoeffs    = coeffs;
1009             final int first       = firstIndex;
1010             final double duration = chunksDuration;
1011             for (int i = 0; i < nbChunks; ++i) {
1012 
1013                 // set up chunk validity range
1014                 final AbsoluteDate chunkStart = chunkEnd;
1015                 chunkEnd = (i == nbChunks - 1) ? rangeEnd : rangeStart.shiftedBy((i + 1) * duration);
1016 
1017                 // extract Chebyshev coefficients for the selected body
1018                 // and convert them from kilometers to meters
1019                 final double[] xCoeffs = new double[nbCoeffs];
1020                 final double[] yCoeffs = new double[nbCoeffs];
1021                 final double[] zCoeffs = new double[nbCoeffs];
1022 
1023                 for (int k = 0; k < nbCoeffs; ++k) {
1024                     // by now, only use the position components
1025                     // if there are also velocity components contained in the file, ignore them
1026                     final int index = first + components * i * nbCoeffs + k - 1;
1027                     xCoeffs[k] = positionUnit * extractDouble(record, 8 * index);
1028                     yCoeffs[k] = positionUnit * extractDouble(record, 8 * (index +  nbCoeffs));
1029                     zCoeffs[k] = positionUnit * extractDouble(record, 8 * (index + 2 * nbCoeffs));
1030                 }
1031 
1032                 // build the position-velocity model for current chunk
1033                 entries.add(new PosVelChebyshev(chunkStart, timeScale, duration, xCoeffs, yCoeffs, zCoeffs));
1034 
1035             }
1036 
1037             return rangeStart;
1038 
1039         }
1040 
1041     }
1042 
1043     /** Raw position-velocity provider using ephemeris. */
1044     private class EphemerisRawPVProvider implements RawPVProvider {
1045 
1046         /** {@inheritDoc} */
1047         public PVCoordinates getRawPV(final AbsoluteDate date) {
1048 
1049             // get raw PV from Chebyshev polynomials
1050             PosVelChebyshev chebyshev;
1051             try {
1052                 chebyshev = ephemerides.getNeighbors(date).findFirst().get();
1053             } catch (TimeStampedCacheException tce) {
1054                 // we cannot bracket the date, check if the last available chunk covers the specified date
1055                 chebyshev = ephemerides.getLatest();
1056                 if (!chebyshev.inRange(date)) {
1057                     // we were not able to recover from the error, the date is too far
1058                     throw tce;
1059                 }
1060             }
1061 
1062             // evaluate the Chebyshev polynomials
1063             return chebyshev.getPositionVelocityAcceleration(date);
1064 
1065         }
1066 
1067         /** {@inheritDoc} */
1068         public <T extends RealFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {
1069 
1070             // get raw PV from Chebyshev polynomials
1071             PosVelChebyshev chebyshev;
1072             try {
1073                 chebyshev = ephemerides.getNeighbors(date.toAbsoluteDate()).findFirst().get();
1074             } catch (TimeStampedCacheException tce) {
1075                 // we cannot bracket the date, check if the last available chunk covers the specified date
1076                 chebyshev = ephemerides.getLatest();
1077                 if (!chebyshev.inRange(date.toAbsoluteDate())) {
1078                     // we were not able to recover from the error, the date is too far
1079                     throw tce;
1080                 }
1081             }
1082 
1083             // evaluate the Chebyshev polynomials
1084             return chebyshev.getPositionVelocityAcceleration(date);
1085 
1086         }
1087 
1088     }
1089 
1090     /** Raw position-velocity provider providing always zero. */
1091     private static class ZeroRawPVProvider implements RawPVProvider {
1092 
1093         /** {@inheritDoc} */
1094         public PVCoordinates getRawPV(final AbsoluteDate date) {
1095             return PVCoordinates.ZERO;
1096         }
1097 
1098         /** {@inheritDoc} */
1099         public <T extends RealFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {
1100             return FieldPVCoordinates.getZero(date.getField());
1101         }
1102 
1103     }
1104 
1105 }