1   /* Copyright 2002-2026 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.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  /** Loader for JPL ephemerides binary files (DE 4xx) and similar formats (INPOP 06/08/10).
63   * <p>JPL ephemerides binary files contain ephemerides for all solar system planets.</p>
64   * <p>The JPL ephemerides binary files are recognized thanks to their base names,
65   * which must match the pattern <code>[lu]nx[mp]####.ddd</code> (or
66   * <code>[lu]nx[mp]####.ddd.gz</code> for gzip-compressed files) where # stands for a
67   * digit character and where ddd is an ephemeris type (typically 405 or 406).</p>
68   * <p>The loader supports files encoded in big-endian as well as in little-endian notation.
69   * Usually, big-endian files are named <code>unx[mp]####.ddd</code>, while little-endian files
70   * are named <code>lnx[mp]####.ddd</code>.</p>
71   * <p>The IMCCE ephemerides binary files are recognized thanks to their base names,
72   * which must match the pattern <code>inpop*.dat</code> (or
73   * <code>inpop*.dat.gz</code> for gzip-compressed files) where * stands for any string.</p>
74   * <p>The loader supports files encoded in big-endian as well as in little-endian notation.
75   * Usually, big-endian files contain <code>bigendian</code> in their names, while little-endian files
76   * contain <code>littleendian</code> in their names.</p>
77   * <p>The loader supports files in TDB or TCB time scales.</p>
78   * @author Luc Maisonobe
79   */
80  public class JPLEphemeridesLoader extends AbstractSelfFeedingLoader
81          implements CelestialBodyLoader {
82  
83      /** Default supported files name pattern for JPL DE files. */
84      public static final String DEFAULT_DE_SUPPORTED_NAMES = "^[lu]nx([mp](\\d{4,5}))+\\.(?:4\\d\\d)$";
85  
86      /** Default supported files name pattern for JPL DE files. */
87      public static final String DEFAULT_DE_2021_SUPPORTED_NAMES = "^linux_([mp](\\d{4,5}))+\\.(?:4\\d\\d)$";
88  
89      /** Default supported files name pattern for IMCCE INPOP files. */
90      public static final String DEFAULT_INPOP_SUPPORTED_NAMES = "^inpop.*\\.dat$";
91  
92      /** 50 days. */
93      private static final TimeOffset FIFTY_DAYS =
94              new TimeOffset(50, TimeUnit.DAYS);
95  
96      /** DE number used by INPOP files. */
97      private static final int INPOP_DE_NUMBER = 100;
98  
99      /** Maximal number of constants in headers. */
100     private static final int CONSTANTS_MAX_NUMBER           = 400;
101 
102     /** Offset of the ephemeris type in first header record. */
103     private static final int HEADER_EPHEMERIS_TYPE_OFFSET   = 2840;
104 
105     /** Offset of the record size (for INPOP files) in first header record. */
106     private static final int HEADER_RECORD_SIZE_OFFSET      = 2856;
107 
108     /** Offset of the start epoch in first header record. */
109     private static final int HEADER_START_EPOCH_OFFSET      = 2652;
110 
111     /** Offset of the end epoch in first header record. */
112     private static final int HEADER_END_EPOCH_OFFSET        = 2660;
113 
114     /** Offset of the astronomical unit in first header record. */
115     private static final int HEADER_ASTRONOMICAL_UNIT_OFFSET = 2680;
116 
117     /** Offset of the Earth-Moon mass ratio in first header record. */
118     private static final int HEADER_EM_RATIO_OFFSET         = 2688;
119 
120     /** Offset of Chebishev coefficients indices in first header record. */
121     private static final int HEADER_CHEBISHEV_INDICES_OFFSET = 2696;
122 
123     /** Offset of libration coefficients indices in first header record. */
124     private static final int HEADER_LIBRATION_INDICES_OFFSET = 2844;
125 
126     /** Offset of chunks duration in first header record. */
127     private static final int HEADER_CHUNK_DURATION_OFFSET    = 2668;
128 
129     /** Offset of the constants names in first header record. */
130     private static final int HEADER_CONSTANTS_NAMES_OFFSET  = 252;
131 
132     /** Offset of the constants values in second header record. */
133     private static final int HEADER_CONSTANTS_VALUES_OFFSET = 0;
134 
135     /** Offset of the range start in the data records. */
136     private static final int DATA_START_RANGE_OFFSET        = 0;
137 
138     /** Offset of the range end in the data records. */
139     private static final int DATE_END_RANGE_OFFSET          = 8;
140 
141     /** The constant name for the astronomical unit. */
142     private static final String CONSTANT_AU = "AU";
143 
144     /** The constant name for the earth-moon mass ratio. */
145     private static final String CONSTANT_EMRAT = "EMRAT";
146 
147     /** List of supported ephemerides types. */
148     public enum EphemerisType {
149 
150         /** Constant for solar system barycenter. */
151         SOLAR_SYSTEM_BARYCENTER,
152 
153         /** Constant for the Sun. */
154         SUN,
155 
156         /** Constant for Mercury. */
157         MERCURY,
158 
159         /** Constant for Venus. */
160         VENUS,
161 
162         /** Constant for the Earth-Moon barycenter. */
163         EARTH_MOON,
164 
165         /** Constant for the Earth. */
166         EARTH,
167 
168         /** Constant for the Moon. */
169         MOON,
170 
171         /** Constant for Mars. */
172         MARS,
173 
174         /** Constant for Jupiter. */
175         JUPITER,
176 
177         /** Constant for Saturn. */
178         SATURN,
179 
180         /** Constant for Uranus. */
181         URANUS,
182 
183         /** Constant for Neptune. */
184         NEPTUNE,
185 
186         /** Constant for Pluto. */
187         PLUTO,
188 
189         /** Constant for Lunar librations. */
190         LIBRATION
191 
192     }
193 
194     /** Interface for raw position-velocity retrieval. */
195     public interface RawPVProvider {
196 
197         /** Get the position-velocity at date.
198          * @param date date at which the position-velocity is desired
199          * @return position-velocity at the specified date
200          */
201         PVCoordinates getRawPV(AbsoluteDate date);
202 
203         /** Get the position at date.
204          * @param date date at which the position is desired
205          * @return position at the specified date
206          */
207         default Vector3D getRawPosition(final AbsoluteDate date) {
208             return getRawPV(date).getPosition();
209         }
210 
211         /** Get the position-velocity at date.
212          * @param date date at which the position-velocity is desired
213          * @param <T> type of the field elements
214          * @return position-velocity at the specified date
215          */
216         <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(FieldAbsoluteDate<T> date);
217 
218         /** Get the position at date.
219          * @param date date at which the position is desired
220          * @param <T> type of the field elements
221          * @return position at the specified date
222          */
223         default <T extends CalculusFieldElement<T>> FieldVector3D<T> getRawPosition(final FieldAbsoluteDate<T> date) {
224             return getRawPV(date).getPosition();
225         }
226 
227     }
228 
229     /** Ephemeris for selected body. */
230     private final GenericTimeStampedCache<PosVelChebyshev> ephemerides;
231 
232     /** Constants defined in the file. */
233     private final AtomicReference<Map<String, Double>> constants;
234 
235     /** Ephemeris type to generate. */
236     private final EphemerisType generateType;
237 
238     /** Ephemeris type to load. */
239     private final EphemerisType loadType;
240 
241     /** Time scales to use when loading data. */
242     private final TimeScales timeScales;
243 
244     /** The GCRF implementation. */
245     private final Frame gcrf;
246 
247     /** Used to retrieve ICRF and EMB frames for building other bodies. */
248     private final CelestialBodies celestialBodies;
249 
250     /** Current file start epoch. */
251     private AbsoluteDate startEpoch;
252 
253     /** Current file final epoch. */
254     private AbsoluteDate finalEpoch;
255 
256     /** Chunks duration (in seconds). */
257     private double maxChunksDuration;
258 
259     /** Current file chunks duration (in seconds). */
260     private TimeOffset chunksDuration;
261 
262     /** Index of the first data for selected body. */
263     private int firstIndex;
264 
265     /** Number of coefficients for selected body. */
266     private int coeffs;
267 
268     /** Number of chunks for the selected body. */
269     private int chunks;
270 
271     /** Number of components contained in the file. */
272     private int components;
273 
274     /** Unit of the position coordinates (as a multiple of meters). */
275     private double positionUnit;
276 
277     /** Time scale of the date coordinates. */
278     private TimeScale timeScale;
279 
280     /** Indicator for binary file endianness. */
281     private boolean bigEndian;
282 
283     /** Create a loader for JPL ephemerides binary files. This constructor uses the {@link
284      * DataContext#getDefault() default data context}.
285      *
286      * @param supportedNames regular expression for supported files names
287      * @param generateType ephemeris type to generate
288      * @see #JPLEphemeridesLoader(String, EphemerisType, DataProvidersManager, TimeScales,
289      * Frame, CelestialBodies)
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      * Create a loader for JPL ephemerides binary files.
302      *
303      * <p>Requiring {@code celestialBodies} to be passed to this constructor is
304      * not a great design when this class is used to load the data for
305      * {@code celestialBodies}. But this loader needs to know the ICRF and Earth
306      * Moon Barycenter frames used by {@code celestialBodies} to avoid creating
307      * duplicate frames, as was done in Orekit 13 and before.
308      *
309      * @param supportedNames       regular expression for supported files names
310      * @param generateType         ephemeris type to generate
311      * @param dataProvidersManager provides access to the ephemeris files.
312      * @param timeScales           used to access the TCB and TDB time scales
313      *                             while loading data.
314      * @param gcrf                 Earth centered frame aligned with ICRF.
315      * @param celestialBodies      used to load the ICRF and Earth-Moon
316      *                             Barycenter frames.
317      * @since 14.0
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     /** Load celestial body.
351      * @param name name of the celestial body
352      * @return loaded celestial body
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         // build the celestial body
397         return new JPLCelestialBody(name, getSupportedNames(), generateType, rawPVProvider,
398                                     gm, scale, iauPole, definingFrameAlignedWithICRF);
399 
400     }
401 
402     /** Load libration.
403      * @return loaded libration
404      * @since 14.0
405      */
406     public JPLLibration loadLibration() {
407         final RawPVProvider rawPVProvider = new EphemerisRawPVProvider();
408 
409         // build the libration
410         return new JPLLibration(rawPVProvider);
411     }
412 
413     /** Get astronomical unit.
414      * @return astronomical unit in meters
415      */
416     public double getLoadedAstronomicalUnit() {
417         return UnitsConverter.KILOMETRES_TO_METRES.convert(getLoadedConstant(CONSTANT_AU));
418     }
419 
420     /** Get Earth/Moon mass ratio.
421      * @return Earth/Moon mass ratio
422      */
423     public double getLoadedEarthMoonMassRatio() {
424         return getLoadedConstant(CONSTANT_EMRAT);
425     }
426 
427     /** Get the gravitational coefficient of a body.
428      * @param body body for which the gravitational coefficient is requested
429      * @return gravitational coefficient in m³/s²
430      */
431     public double getLoadedGravitationalCoefficient(final EphemerisType body) {
432 
433         // coefficient in au³/day²
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     /** Get a constant defined in the ephemerides headers.
493      * <p>Note that since constants are defined in the JPL headers
494      * files, they are available as soon as one file is available, even
495      * if it doesn't match the desired central date. This is because the
496      * header must be parsed before the dates can be checked.</p>
497      * <p>
498      * There are alternate names for constants since for example JPL names are
499      * different from INPOP names (Sun gravity: GMS or GM_Sun, Mars gravity:
500      * GM4 or GM_Mar...).
501      * </p>
502      * @param names alternate names of the constant
503      * @return value of the constant of NaN if the constant is not defined
504      */
505     public double getLoadedConstant(final String... names) {
506 
507         // lazy loading of constants
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     /** Get the maximal chunks duration.
529      * @return chunks maximal duration in seconds
530      */
531     public double getMaxChunksDuration() {
532         return maxChunksDuration;
533     }
534 
535     /** Parse the first header record.
536      * @param record first header record
537      * @param name name of the file (or zip entry)
538      */
539     private void parseFirstHeaderRecord(final byte[] record, final String name) {
540 
541         // get the ephemerides type
542         final int deNum = extractInt(record, HEADER_EPHEMERIS_TYPE_OFFSET);
543 
544         // as default, 3 polynomial coefficients for the Cartesian coordinates
545         // (x, y, z) are contained in the file, positions are in kilometers
546         // and times are in TDB
547         components   = 3;
548         positionUnit = UnitsConverter.KILOMETRES_TO_METRES.convert(1.0);
549         timeScale    = timeScales.getTDB();
550 
551         if (deNum == INPOP_DE_NUMBER) {
552             // an INPOP file may contain 6 components (including coefficients for the velocity vector)
553             final double format = getLoadedConstant("FORMAT");
554             if (!Double.isNaN(format) && (int) FastMath.IEEEremainder(format, 10) != 1) {
555                 components = 6;
556             }
557 
558             // INPOP files may have their polynomials expressed in AU
559             final double unite = getLoadedConstant("UNITE");
560             if (!Double.isNaN(unite) && (int) unite == 0) {
561                 positionUnit = getLoadedAstronomicalUnit();
562             }
563 
564             // INPOP files may have their times expressed in TCB
565             final double timesc = getLoadedConstant("TIMESC");
566             if (!Double.isNaN(timesc) && (int) timesc == 1) {
567                 timeScale = timeScales.getTCB();
568             }
569 
570         }
571 
572         // extract covered date range
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             // indices of the Chebyshev coefficients for lunar librations
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             // indices of the Chebyshev coefficients for each ephemeris
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         // compute chunks duration
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         // sanity checks
622         if (!ok) {
623             throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
624         }
625 
626     }
627 
628     /** Read first header record.
629      * @param input input stream
630      * @param name name of the file (or zip entry)
631      * @return record record where to put bytes
632      * @exception IOException if a read error occurs
633      */
634     private byte[] readFirstRecord(final InputStream input, final String name)
635         throws IOException {
636 
637         // read first part of record, up to the ephemeris type
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         // detect the endian format
644         detectEndianess(firstPart);
645 
646         // get the ephemerides type
647         final int deNum = extractInt(firstPart, HEADER_EPHEMERIS_TYPE_OFFSET);
648 
649         // the record size for this file
650         final int recordSize;
651 
652         if (deNum == INPOP_DE_NUMBER) {
653             // INPOP files have an extended DE format, which includes also the record size
654             recordSize = extractInt(firstPart, HEADER_RECORD_SIZE_OFFSET) << 3;
655         } else {
656             // compute the record size for original JPL files
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         // build a record with the proper size and finish read of the first complete record
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     /** Parse constants from first two header records.
677      * @param first first header record
678      * @param second second header record
679      * @return map of parsed constants
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             // Note: for extracting the strings from the binary file, it makes no difference
687             //       if the file is stored in big-endian or little-endian notation
688             final String constantName = extractString(first, HEADER_CONSTANTS_NAMES_OFFSET + i * 6, 6);
689             if (constantName.length() == 0) {
690                 // no more constants to read
691                 break;
692             }
693             final double constantValue = extractDouble(second, HEADER_CONSTANTS_VALUES_OFFSET + 8 * i);
694             map.put(constantName, constantValue);
695         }
696 
697         // INPOP files do not have constants for AU and EMRAT, thus extract them from
698         // the header record and create a constant for them to be consistent with JPL files
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     /** Read bytes into the current record array.
712      * @param input input stream
713      * @param record record where to put bytes
714      * @param start start index where to put bytes
715      * @return true if record has been filled up
716      * @exception IOException if a read error occurs
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     /** Detect whether the JPL ephemerides file is stored in big-endian or
732      * little-endian notation.
733      * @param record the array containing the binary JPL header
734      */
735     private void detectEndianess(final byte[] record) {
736 
737         // default to big-endian
738         bigEndian = true;
739 
740         // first try to read the DE number in big-endian format
741         // the number is stored as unsigned int, so we have to convert it properly
742         final long deNum = extractInt(record, HEADER_EPHEMERIS_TYPE_OFFSET) & 0xffffffffL;
743 
744         // simple heuristic: if the read value is larger than half the range of an integer
745         //                   assume the file is in little-endian format
746         if (deNum > (1 << 15)) {
747             bigEndian = false;
748         }
749 
750     }
751 
752     /** Calculate the record size of a JPL ephemerides file.
753      * @param record the byte array containing the header record
754      * @param name the name of the data file
755      * @return the record size for this file
756      */
757     private int computeRecordSize(final byte[] record, final String name) {
758 
759         int recordSize = 0;
760         boolean ok = true;
761         // JPL files always have 3 position components
762         final int nComp = 3;
763 
764         // iterate over the coefficient ptr array and sum up the record size
765         // the coeffPtr array has the dimensions [12][nComp]
766         for (int j = 0; j < 12; j++) {
767             final int nCompCur = (j == 11) ? 2 : nComp;
768 
769             // Note: the array element coeffPtr[j][0] is not needed for the calculation
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             // sanity checks
775             ok = ok && (coeffPtr1 >= 0 || coeffPtr2 >= 0);
776 
777             recordSize += coeffPtr1 * coeffPtr2 * nCompCur;
778         }
779 
780         // the libration ptr array has the dimension [3]
781         // Note: the array element libratPtr[0] is not needed for the calculation
782         final int libratPtr1 = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 4);
783         final int libratPtr2 = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 8);
784 
785         // sanity checks
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     /** Extract a date from a record.
800      * @param record record to parse
801      * @param offset offset of the double within the record
802      * @return extracted date
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     /** Extract a double from a record.
818      * <p>Double numbers are stored according to IEEE 754 standard, with
819      * most significant byte first.</p>
820      * @param record record to parse
821      * @param offset offset of the double within the record
822      * @return extracted double
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     /** Extract an int from a record.
845      * @param record record to parse
846      * @param offset offset of the double within the record
847      * @return extracted int
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     /** Extract a String from a record.
863      * @param record record to parse
864      * @param offset offset of the string within the record
865      * @param length maximal length of the string
866      * @return extracted string, with whitespace characters stripped
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     /** Local parser for header constants. */
873     private class ConstantsParser implements DataLoader {
874 
875         /** Local constants map. */
876         private Map<String, Double> localConstants;
877 
878        /** Get the local constants map.
879          * @return local constants map
880          */
881         public Map<String, Double> getConstants() {
882             return localConstants;
883         }
884 
885         /** {@inheritDoc} */
886         public boolean stillAcceptsData() {
887             return localConstants == null;
888         }
889 
890         /** {@inheritDoc} */
891         public void loadData(final InputStream input, final String name)
892             throws IOException, ParseException, OrekitException {
893 
894             // read first header record
895             final byte[] first = readFirstRecord(input, name);
896 
897             // the second record contains the values of the constants used for least-square filtering
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     /** Local parser for Chebyshev polynomials. */
910     private class EphemerisParser implements DataLoader, TimeStampedGenerator<PosVelChebyshev> {
911 
912         /** Set of Chebyshev polynomials read. */
913         private final SortedSet<PosVelChebyshev> entries;
914 
915         /** Start of range we are interested in. */
916         private AbsoluteDate start;
917 
918         /** End of range we are interested in. */
919         private AbsoluteDate end;
920 
921         /** Simple constructor.
922          */
923         EphemerisParser() {
924             entries = new TreeSet<>(new ChronologicalComparator());
925         }
926 
927         /** {@inheritDoc} */
928         public List<PosVelChebyshev> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {
929             try {
930 
931                 // prepare reading
932                 entries.clear();
933                 if (existingDate == null) {
934                     // we want ephemeris data for the first time, set up an arbitrary first range
935                     start = date.shiftedBy(FIFTY_DAYS.negate());
936                     end   = date.shiftedBy(FIFTY_DAYS);
937                 } else if (existingDate.compareTo(date) <= 0) {
938                     // we want to extend an existing range towards future dates
939                     start = existingDate;
940                     end   = date;
941                 } else {
942                     // we want to extend an existing range towards past dates
943                     start = date;
944                     end   = existingDate;
945                 }
946 
947                 // get new entries in the specified data range
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         /** {@inheritDoc} */
960         public boolean stillAcceptsData() {
961 
962             // special case for Earth: we do not really load any ephemeris data
963             if (generateType == EphemerisType.EARTH) {
964                 return false;
965             }
966 
967             // we have to look for data in all available ephemerides files as there may be
968             // data overlaps that result in incomplete data
969             if (entries.isEmpty()) {
970                 return true;
971             } else {
972                 // if the requested range is already filled, we do not need to look further
973                 return !(entries.getFirst().getDate().compareTo(start) < 0 &&
974                          entries.getLast().getDate().compareTo(end)    > 0);
975             }
976 
977         }
978 
979         /** {@inheritDoc} */
980         public void loadData(final InputStream input, final String name)
981             throws IOException {
982 
983             // read first header record
984             final byte[] first = readFirstRecord(input, name);
985 
986             // the second record contains the values of the constants used for least-square filtering
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             // check astronomical unit consistency
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             // check Earth-Moon mass ratio consistency
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             // parse first header record
1017             parseFirstHeaderRecord(first, name);
1018 
1019             if (startEpoch.compareTo(end) < 0 && finalEpoch.compareTo(start) > 0) {
1020                 // this file contains data in the range we are looking for, read it
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                         // we have already exceeded the range we were interested in,
1026                         // we interrupt parsing here
1027                         return;
1028                     }
1029                 }
1030             }
1031 
1032         }
1033 
1034         /** Parse regular ephemeris record.
1035          * @param record record to parse
1036          * @return date of the last parsed chunk
1037          */
1038         private AbsoluteDate parseDataRecord(final byte[] record) {
1039 
1040             // extract time range covered by the record
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                 // we are not interested in this record, don't parse it
1055                 return rangeEnd;
1056             }
1057 
1058             // loop over chunks inside the time range
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                 // set up chunk validity range
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                 // extract Chebyshev coefficients for the selected body
1078                 // and convert them from kilometers to meters
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                     // by now, only use the position components
1085                     // if there are also velocity components contained in the file, ignore them
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                 // build the position-velocity model for current chunk
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     /** Raw position-velocity provider using ephemeris. */
1105     private class EphemerisRawPVProvider implements RawPVProvider {
1106 
1107 
1108         /** Retrieve correct Chebyshev polynomial for given date.
1109          * @param date date
1110          * @return PosVelChebyshev Chebyshev polynomial class for position-velocity-acceleration
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                 // we cannot bracket the date, check if the last available chunk covers the specified date
1118                 chebyshev = ephemerides.getLatest();
1119                 if (!chebyshev.inRange(date)) {
1120                     // we were not able to recover from the error, the date is too far
1121                     throw tce;
1122                 }
1123             }
1124             return chebyshev;
1125         }
1126 
1127         /** {@inheritDoc} */
1128         public PVCoordinates getRawPV(final AbsoluteDate date) {
1129 
1130             // get raw PV from Chebyshev polynomials
1131             final PosVelChebyshev chebyshev = getChebyshev(date);
1132 
1133             // evaluate the Chebyshev polynomials
1134             return chebyshev.getPositionVelocityAcceleration(date);
1135 
1136         }
1137 
1138         /** {@inheritDoc} */
1139         public <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {
1140 
1141             // get raw PV from Chebyshev polynomials
1142             final PosVelChebyshev chebyshev = getChebyshev(date.toAbsoluteDate());
1143 
1144             // evaluate the Chebyshev polynomials
1145             return chebyshev.getPositionVelocityAcceleration(date);
1146 
1147         }
1148 
1149         /** {@inheritDoc} */
1150         @Override
1151         public Vector3D getRawPosition(final AbsoluteDate date) {
1152             // get raw PV from Chebyshev polynomials
1153             final PosVelChebyshev chebyshev = getChebyshev(date);
1154 
1155             // evaluate the Chebyshev polynomials
1156             return chebyshev.getPosition(date);
1157         }
1158 
1159         /** {@inheritDoc} */
1160         @Override
1161         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getRawPosition(final FieldAbsoluteDate<T> date) {
1162             // get raw PV from Chebyshev polynomials
1163             final PosVelChebyshev chebyshev = getChebyshev(date.toAbsoluteDate());
1164 
1165             // evaluate the Chebyshev polynomials
1166             return chebyshev.getPosition(date);
1167         }
1168 
1169     }
1170 
1171     /** Raw position-velocity provider providing always zero. */
1172     static class ZeroRawPVProvider implements RawPVProvider {
1173 
1174         /** {@inheritDoc} */
1175         public PVCoordinates getRawPV(final AbsoluteDate date) {
1176             return PVCoordinates.ZERO;
1177         }
1178 
1179         /** {@inheritDoc} */
1180         public <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {
1181             return FieldPVCoordinates.getZero(date.getField());
1182         }
1183 
1184     }
1185 
1186 }