JPLEphemeridesLoader.java

  1. /* Copyright 2002-2024 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. import java.io.IOException;
  19. import java.io.InputStream;
  20. import java.nio.charset.StandardCharsets;
  21. import java.text.ParseException;
  22. import java.util.ArrayList;
  23. import java.util.HashMap;
  24. import java.util.List;
  25. import java.util.Map;
  26. import java.util.SortedSet;
  27. import java.util.TreeSet;
  28. import java.util.concurrent.atomic.AtomicReference;

  29. import org.hipparchus.CalculusFieldElement;
  30. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  31. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  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. import org.orekit.utils.units.UnitsConverter;

  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.     /** Default supported files name pattern for JPL DE files. */
  79.     public static final String DEFAULT_DE_SUPPORTED_NAMES = "^[lu]nx([mp](\\d\\d\\d\\d))+\\.(?:4\\d\\d)$";

  80.     /** Default supported files name pattern for IMCCE INPOP files. */
  81.     public static final String DEFAULT_INPOP_SUPPORTED_NAMES = "^inpop.*\\.dat$";

  82.     /** 50 days in seconds. */
  83.     private static final double FIFTY_DAYS = 50 * Constants.JULIAN_DAY;

  84.     /** DE number used by INPOP files. */
  85.     private static final int INPOP_DE_NUMBER = 100;

  86.     /** Maximal number of constants in headers. */
  87.     private static final int CONSTANTS_MAX_NUMBER           = 400;

  88.     /** Offset of the ephemeris type in first header record. */
  89.     private static final int HEADER_EPHEMERIS_TYPE_OFFSET   = 2840;

  90.     /** Offset of the record size (for INPOP files) in first header record. */
  91.     private static final int HEADER_RECORD_SIZE_OFFSET      = 2856;

  92.     /** Offset of the start epoch in first header record. */
  93.     private static final int HEADER_START_EPOCH_OFFSET      = 2652;

  94.     /** Offset of the end epoch in first header record. */
  95.     private static final int HEADER_END_EPOCH_OFFSET        = 2660;

  96.     /** Offset of the astronomical unit in first header record. */
  97.     private static final int HEADER_ASTRONOMICAL_UNIT_OFFSET = 2680;

  98.     /** Offset of the Earth-Moon mass ratio in first header record. */
  99.     private static final int HEADER_EM_RATIO_OFFSET         = 2688;

  100.     /** Offset of Chebishev coefficients indices in first header record. */
  101.     private static final int HEADER_CHEBISHEV_INDICES_OFFSET = 2696;

  102.     /** Offset of libration coefficients indices in first header record. */
  103.     private static final int HEADER_LIBRATION_INDICES_OFFSET = 2844;

  104.     /** Offset of chunks duration in first header record. */
  105.     private static final int HEADER_CHUNK_DURATION_OFFSET    = 2668;

  106.     /** Offset of the constants names in first header record. */
  107.     private static final int HEADER_CONSTANTS_NAMES_OFFSET  = 252;

  108.     /** Offset of the constants values in second header record. */
  109.     private static final int HEADER_CONSTANTS_VALUES_OFFSET = 0;

  110.     /** Offset of the range start in the data records. */
  111.     private static final int DATA_START_RANGE_OFFSET        = 0;

  112.     /** Offset of the range end in the data records. */
  113.     private static final int DATE_END_RANGE_OFFSET          = 8;

  114.     /** The constant name for the astronomical unit. */
  115.     private static final String CONSTANT_AU = "AU";

  116.     /** The constant name for the earth-moon mass ratio. */
  117.     private static final String CONSTANT_EMRAT = "EMRAT";

  118.     /** List of supported ephemerides types. */
  119.     public enum EphemerisType {

  120.         /** Constant for solar system barycenter. */
  121.         SOLAR_SYSTEM_BARYCENTER,

  122.         /** Constant for the Sun. */
  123.         SUN,

  124.         /** Constant for Mercury. */
  125.         MERCURY,

  126.         /** Constant for Venus. */
  127.         VENUS,

  128.         /** Constant for the Earth-Moon barycenter. */
  129.         EARTH_MOON,

  130.         /** Constant for the Earth. */
  131.         EARTH,

  132.         /** Constant for the Moon. */
  133.         MOON,

  134.         /** Constant for Mars. */
  135.         MARS,

  136.         /** Constant for Jupiter. */
  137.         JUPITER,

  138.         /** Constant for Saturn. */
  139.         SATURN,

  140.         /** Constant for Uranus. */
  141.         URANUS,

  142.         /** Constant for Neptune. */
  143.         NEPTUNE,

  144.         /** Constant for Pluto. */
  145.         PLUTO

  146.     }

  147.     /** Interface for raw position-velocity retrieval. */
  148.     public interface RawPVProvider {

  149.         /** Get the position-velocity at date.
  150.          * @param date date at which the position-velocity is desired
  151.          * @return position-velocity at the specified date
  152.          */
  153.         PVCoordinates getRawPV(AbsoluteDate date);

  154.         /** Get the position at date.
  155.          * @param date date at which the position is desired
  156.          * @return position at the specified date
  157.          */
  158.         default Vector3D getRawPosition(final AbsoluteDate date) {
  159.             return getRawPV(date).getPosition();
  160.         }

  161.         /** Get the position-velocity at date.
  162.          * @param date date at which the position-velocity is desired
  163.          * @param <T> type of the field elements
  164.          * @return position-velocity at the specified date
  165.          */
  166.         <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(FieldAbsoluteDate<T> date);

  167.         /** Get the position at date.
  168.          * @param date date at which the position is desired
  169.          * @param <T> type of the field elements
  170.          * @return position at the specified date
  171.          */
  172.         default <T extends CalculusFieldElement<T>> FieldVector3D<T> getRawPosition(final FieldAbsoluteDate<T> date) {
  173.             return getRawPV(date).getPosition();
  174.         }

  175.     }

  176.     /** Ephemeris for selected body. */
  177.     private final GenericTimeStampedCache<PosVelChebyshev> ephemerides;

  178.     /** Constants defined in the file. */
  179.     private final AtomicReference<Map<String, Double>> constants;

  180.     /** Ephemeris type to generate. */
  181.     private final EphemerisType generateType;

  182.     /** Ephemeris type to load. */
  183.     private final EphemerisType loadType;

  184.     /** Time scales to use when loading data. */
  185.     private final TimeScales timeScales;

  186.     /** The GCRF implementation. */
  187.     private final Frame gcrf;

  188.     /** Current file start epoch. */
  189.     private AbsoluteDate startEpoch;

  190.     /** Current file final epoch. */
  191.     private AbsoluteDate finalEpoch;

  192.     /** Chunks duration (in seconds). */
  193.     private double maxChunksDuration;

  194.     /** Current file chunks duration (in seconds). */
  195.     private double chunksDuration;

  196.     /** Index of the first data for selected body. */
  197.     private int firstIndex;

  198.     /** Number of coefficients for selected body. */
  199.     private int coeffs;

  200.     /** Number of chunks for the selected body. */
  201.     private int chunks;

  202.     /** Number of components contained in the file. */
  203.     private int components;

  204.     /** Unit of the position coordinates (as a multiple of meters). */
  205.     private double positionUnit;

  206.     /** Time scale of the date coordinates. */
  207.     private TimeScale timeScale;

  208.     /** Indicator for binary file endianness. */
  209.     private boolean bigEndian;

  210.     /** Create a loader for JPL ephemerides binary files. This constructor uses the {@link
  211.      * DataContext#getDefault() default data context}.
  212.      *
  213.      * @param supportedNames regular expression for supported files names
  214.      * @param generateType ephemeris type to generate
  215.      * @see #JPLEphemeridesLoader(String, EphemerisType, DataProvidersManager, TimeScales,
  216.      * Frame)
  217.      */
  218.     @DefaultDataContext
  219.     public JPLEphemeridesLoader(final String supportedNames, final EphemerisType generateType) {
  220.         this(supportedNames, generateType,
  221.                 DataContext.getDefault().getDataProvidersManager(),
  222.                 DataContext.getDefault().getTimeScales(),
  223.                 DataContext.getDefault().getFrames().getGCRF());
  224.     }

  225.     /** Create a loader for JPL ephemerides binary files.
  226.      * @param supportedNames regular expression for supported files names
  227.      * @param generateType ephemeris type to generate
  228.      * @param dataProvidersManager provides access to the ephemeris files.
  229.      * @param timeScales used to access the TCB and TDB time scales while loading data.
  230.      * @param gcrf Earth centered frame aligned with ICRF.
  231.      * @since 10.1
  232.      */
  233.     public JPLEphemeridesLoader(final String supportedNames,
  234.                                 final EphemerisType generateType,
  235.                                 final DataProvidersManager dataProvidersManager,
  236.                                 final TimeScales timeScales,
  237.                                 final Frame gcrf) {
  238.         super(supportedNames, dataProvidersManager);

  239.         this.timeScales = timeScales;
  240.         this.gcrf = gcrf;
  241.         constants = new AtomicReference<>();

  242.         this.generateType  = generateType;
  243.         if (generateType == EphemerisType.SOLAR_SYSTEM_BARYCENTER) {
  244.             loadType = EphemerisType.EARTH_MOON;
  245.         } else if (generateType == EphemerisType.EARTH_MOON) {
  246.             loadType = EphemerisType.MOON;
  247.         } else {
  248.             loadType = generateType;
  249.         }

  250.         ephemerides = new GenericTimeStampedCache<>(
  251.                 2, OrekitConfiguration.getCacheSlotsNumber(),
  252.                 Double.POSITIVE_INFINITY, FIFTY_DAYS,
  253.                 new EphemerisParser());
  254.         maxChunksDuration = Double.NaN;
  255.         chunksDuration    = Double.NaN;

  256.     }

  257.     /** Load celestial body.
  258.      * @param name name of the celestial body
  259.      * @return loaded celestial body
  260.      */
  261.     public CelestialBody loadCelestialBody(final String name) {

  262.         final double gm       = getLoadedGravitationalCoefficient(generateType);
  263.         final IAUPole iauPole = PredefinedIAUPoles
  264.                 .getIAUPole(generateType, timeScales);
  265.         final double scale;
  266.         final Frame definingFrameAlignedWithICRF;
  267.         final RawPVProvider rawPVProvider;
  268.         String inertialFrameName = null;
  269.         String bodyOrientedFrameName = null;
  270.         switch (generateType) {
  271.             case SOLAR_SYSTEM_BARYCENTER : {
  272.                 scale = -1.0;
  273.                 final JPLEphemeridesLoader parentLoader = new JPLEphemeridesLoader(
  274.                         getSupportedNames(),
  275.                         EphemerisType.EARTH_MOON,
  276.                         getDataProvidersManager(),
  277.                         timeScales,
  278.                         gcrf);
  279.                 final CelestialBody parentBody =
  280.                         parentLoader.loadCelestialBody(CelestialBodyFactory.EARTH_MOON);
  281.                 definingFrameAlignedWithICRF = parentBody.getInertiallyOrientedFrame();
  282.                 rawPVProvider = new EphemerisRawPVProvider();
  283.                 inertialFrameName = Predefined.ICRF.getName();
  284.                 bodyOrientedFrameName = null;
  285.                 break;
  286.             }
  287.             case EARTH_MOON :
  288.                 scale         = 1.0 / (1.0 + getLoadedEarthMoonMassRatio());
  289.                 definingFrameAlignedWithICRF = gcrf;
  290.                 rawPVProvider = new EphemerisRawPVProvider();
  291.                 break;
  292.             case EARTH :
  293.                 scale         = 1.0;
  294.                 definingFrameAlignedWithICRF = gcrf;
  295.                 rawPVProvider = new ZeroRawPVProvider();
  296.                 break;
  297.             case MOON :
  298.                 scale         =  1.0;
  299.                 definingFrameAlignedWithICRF = gcrf;
  300.                 rawPVProvider = new EphemerisRawPVProvider();
  301.                 break;
  302.             default : {
  303.                 scale = 1.0;
  304.                 final JPLEphemeridesLoader parentLoader = new JPLEphemeridesLoader(
  305.                         getSupportedNames(),
  306.                         EphemerisType.SOLAR_SYSTEM_BARYCENTER,
  307.                         getDataProvidersManager(),
  308.                         timeScales,
  309.                         gcrf);
  310.                 final CelestialBody parentBody =
  311.                         parentLoader.loadCelestialBody(CelestialBodyFactory.SOLAR_SYSTEM_BARYCENTER);
  312.                 definingFrameAlignedWithICRF = parentBody.getInertiallyOrientedFrame();
  313.                 rawPVProvider = new EphemerisRawPVProvider();
  314.             }
  315.         }

  316.         // build the celestial body
  317.         return new JPLCelestialBody(name, getSupportedNames(), generateType, rawPVProvider,
  318.                                     gm, scale, iauPole, definingFrameAlignedWithICRF,
  319.                                     inertialFrameName, bodyOrientedFrameName);

  320.     }

  321.     /** Get astronomical unit.
  322.      * @return astronomical unit in meters
  323.      */
  324.     public double getLoadedAstronomicalUnit() {
  325.         return UnitsConverter.KILOMETRES_TO_METRES.convert(getLoadedConstant(CONSTANT_AU));
  326.     }

  327.     /** Get Earth/Moon mass ratio.
  328.      * @return Earth/Moon mass ratio
  329.      */
  330.     public double getLoadedEarthMoonMassRatio() {
  331.         return getLoadedConstant(CONSTANT_EMRAT);
  332.     }

  333.     /** Get the gravitational coefficient of a body.
  334.      * @param body body for which the gravitational coefficient is requested
  335.      * @return gravitational coefficient in m³/s²
  336.      */
  337.     public double getLoadedGravitationalCoefficient(final EphemerisType body) {

  338.         // coefficient in au³/day²
  339.         final double rawGM;
  340.         switch (body) {
  341.             case SOLAR_SYSTEM_BARYCENTER :
  342.                 return getLoadedGravitationalCoefficient(EphemerisType.SUN)        +
  343.                         getLoadedGravitationalCoefficient(EphemerisType.MERCURY)    +
  344.                         getLoadedGravitationalCoefficient(EphemerisType.VENUS)      +
  345.                         getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) +
  346.                         getLoadedGravitationalCoefficient(EphemerisType.MARS)       +
  347.                         getLoadedGravitationalCoefficient(EphemerisType.JUPITER)    +
  348.                         getLoadedGravitationalCoefficient(EphemerisType.SATURN)     +
  349.                         getLoadedGravitationalCoefficient(EphemerisType.URANUS)     +
  350.                         getLoadedGravitationalCoefficient(EphemerisType.NEPTUNE)    +
  351.                         getLoadedGravitationalCoefficient(EphemerisType.PLUTO);
  352.             case SUN :
  353.                 rawGM = getLoadedConstant("GMS", "GM_Sun");
  354.                 break;
  355.             case MERCURY :
  356.                 rawGM = getLoadedConstant("GM1", "GM_Mer");
  357.                 break;
  358.             case VENUS :
  359.                 rawGM = getLoadedConstant("GM2", "GM_Ven");
  360.                 break;
  361.             case EARTH_MOON :
  362.                 rawGM = getLoadedConstant("GMB", "GM_EMB");
  363.                 break;
  364.             case EARTH :
  365.                 return getLoadedEarthMoonMassRatio() *
  366.                         getLoadedGravitationalCoefficient(EphemerisType.MOON);
  367.             case MOON :
  368.                 return getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) /
  369.                         (1.0 + getLoadedEarthMoonMassRatio());
  370.             case MARS :
  371.                 rawGM = getLoadedConstant("GM4", "GM_Mar");
  372.                 break;
  373.             case JUPITER :
  374.                 rawGM = getLoadedConstant("GM5", "GM_Jup");
  375.                 break;
  376.             case SATURN :
  377.                 rawGM = getLoadedConstant("GM6", "GM_Sat");
  378.                 break;
  379.             case URANUS :
  380.                 rawGM = getLoadedConstant("GM7", "GM_Ura");
  381.                 break;
  382.             case NEPTUNE :
  383.                 rawGM = getLoadedConstant("GM8", "GM_Nep");
  384.                 break;
  385.             case PLUTO :
  386.                 rawGM = getLoadedConstant("GM9", "GM_Plu");
  387.                 break;
  388.             default :
  389.                 throw new OrekitInternalError(null);
  390.         }

  391.         final double au    = getLoadedAstronomicalUnit();
  392.         return rawGM * au * au * au / (Constants.JULIAN_DAY * Constants.JULIAN_DAY);

  393.     }

  394.     /** Get a constant defined in the ephemerides headers.
  395.      * <p>Note that since constants are defined in the JPL headers
  396.      * files, they are available as soon as one file is available, even
  397.      * if it doesn't match the desired central date. This is because the
  398.      * header must be parsed before the dates can be checked.</p>
  399.      * <p>
  400.      * There are alternate names for constants since for example JPL names are
  401.      * different from INPOP names (Sun gravity: GMS or GM_Sun, Mars gravity:
  402.      * GM4 or GM_Mar...).
  403.      * </p>
  404.      * @param names alternate names of the constant
  405.      * @return value of the constant of NaN if the constant is not defined
  406.      */
  407.     public double getLoadedConstant(final String... names) {

  408.         // lazy loading of constants
  409.         Map<String, Double> map = constants.get();
  410.         if (map == null) {
  411.             final ConstantsParser parser = new ConstantsParser();
  412.             if (!feed(parser)) {
  413.                 throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
  414.             }
  415.             map = parser.getConstants();
  416.             constants.compareAndSet(null, map);
  417.         }

  418.         for (final String name : names) {
  419.             if (map.containsKey(name)) {
  420.                 return map.get(name);
  421.             }
  422.         }

  423.         return Double.NaN;

  424.     }

  425.     /** Get the maximal chunks duration.
  426.      * @return chunks maximal duration in seconds
  427.      */
  428.     public double getMaxChunksDuration() {
  429.         return maxChunksDuration;
  430.     }

  431.     /** Parse the first header record.
  432.      * @param record first header record
  433.      * @param name name of the file (or zip entry)
  434.      */
  435.     private void parseFirstHeaderRecord(final byte[] record, final String name) {

  436.         // get the ephemerides type
  437.         final int deNum = extractInt(record, HEADER_EPHEMERIS_TYPE_OFFSET);

  438.         // as default, 3 polynomial coefficients for the Cartesian coordinates
  439.         // (x, y, z) are contained in the file, positions are in kilometers
  440.         // and times are in TDB
  441.         components   = 3;
  442.         positionUnit = UnitsConverter.KILOMETRES_TO_METRES.convert(1.0);
  443.         timeScale    = timeScales.getTDB();

  444.         if (deNum == INPOP_DE_NUMBER) {
  445.             // an INPOP file may contain 6 components (including coefficients for the velocity vector)
  446.             final double format = getLoadedConstant("FORMAT");
  447.             if (!Double.isNaN(format) && (int) FastMath.IEEEremainder(format, 10) != 1) {
  448.                 components = 6;
  449.             }

  450.             // INPOP files may have their polynomials expressed in AU
  451.             final double unite = getLoadedConstant("UNITE");
  452.             if (!Double.isNaN(unite) && (int) unite == 0) {
  453.                 positionUnit = getLoadedAstronomicalUnit();
  454.             }

  455.             // INPOP files may have their times expressed in TCB
  456.             final double timesc = getLoadedConstant("TIMESC");
  457.             if (!Double.isNaN(timesc) && (int) timesc == 1) {
  458.                 timeScale = timeScales.getTCB();
  459.             }

  460.         }

  461.         // extract covered date range
  462.         startEpoch = extractDate(record, HEADER_START_EPOCH_OFFSET);
  463.         finalEpoch = extractDate(record, HEADER_END_EPOCH_OFFSET);
  464.         boolean ok = finalEpoch.compareTo(startEpoch) > 0;

  465.         // indices of the Chebyshev coefficients for each ephemeris
  466.         for (int i = 0; i < 12; ++i) {
  467.             final int row1 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET     + 12 * i);
  468.             final int row2 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 4 + 12 * i);
  469.             final int row3 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 8 + 12 * i);
  470.             ok = ok && row1 >= 0 && row2 >= 0 && row3 >= 0;
  471.             if (i ==  0 && loadType == EphemerisType.MERCURY    ||
  472.                     i ==  1 && loadType == EphemerisType.VENUS      ||
  473.                     i ==  2 && loadType == EphemerisType.EARTH_MOON ||
  474.                     i ==  3 && loadType == EphemerisType.MARS       ||
  475.                     i ==  4 && loadType == EphemerisType.JUPITER    ||
  476.                     i ==  5 && loadType == EphemerisType.SATURN     ||
  477.                     i ==  6 && loadType == EphemerisType.URANUS     ||
  478.                     i ==  7 && loadType == EphemerisType.NEPTUNE    ||
  479.                     i ==  8 && loadType == EphemerisType.PLUTO      ||
  480.                     i ==  9 && loadType == EphemerisType.MOON       ||
  481.                     i == 10 && loadType == EphemerisType.SUN) {
  482.                 firstIndex = row1;
  483.                 coeffs     = row2;
  484.                 chunks     = row3;
  485.             }
  486.         }

  487.         // compute chunks duration
  488.         final double timeSpan = extractDouble(record, HEADER_CHUNK_DURATION_OFFSET);
  489.         ok = ok && timeSpan > 0 && timeSpan < 100;
  490.         chunksDuration = Constants.JULIAN_DAY * (timeSpan / chunks);
  491.         if (Double.isNaN(maxChunksDuration)) {
  492.             maxChunksDuration = chunksDuration;
  493.         } else {
  494.             maxChunksDuration = FastMath.max(maxChunksDuration, chunksDuration);
  495.         }

  496.         // sanity checks
  497.         if (!ok) {
  498.             throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  499.         }

  500.     }

  501.     /** Read first header record.
  502.      * @param input input stream
  503.      * @param name name of the file (or zip entry)
  504.      * @return record record where to put bytes
  505.      * @exception IOException if a read error occurs
  506.      */
  507.     private byte[] readFirstRecord(final InputStream input, final String name)
  508.         throws IOException {

  509.         // read first part of record, up to the ephemeris type
  510.         final byte[] firstPart = new byte[HEADER_RECORD_SIZE_OFFSET + 4];
  511.         if (!readInRecord(input, firstPart, 0)) {
  512.             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  513.         }

  514.         // detect the endian format
  515.         detectEndianess(firstPart);

  516.         // get the ephemerides type
  517.         final int deNum = extractInt(firstPart, HEADER_EPHEMERIS_TYPE_OFFSET);

  518.         // the record size for this file
  519.         final int recordSize;

  520.         if (deNum == INPOP_DE_NUMBER) {
  521.             // INPOP files have an extended DE format, which includes also the record size
  522.             recordSize = extractInt(firstPart, HEADER_RECORD_SIZE_OFFSET) << 3;
  523.         } else {
  524.             // compute the record size for original JPL files
  525.             recordSize = computeRecordSize(firstPart, name);
  526.         }

  527.         if (recordSize <= 0) {
  528.             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  529.         }

  530.         // build a record with the proper size and finish read of the first complete record
  531.         final int start = firstPart.length;
  532.         final byte[] record = new byte[recordSize];
  533.         System.arraycopy(firstPart, 0, record, 0, firstPart.length);
  534.         if (!readInRecord(input, record, start)) {
  535.             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  536.         }

  537.         return record;

  538.     }

  539.     /** Parse constants from first two header records.
  540.      * @param first first header record
  541.      * @param second second header record
  542.      * @return map of parsed constants
  543.      */
  544.     private Map<String, Double> parseConstants(final byte[] first, final byte[] second) {

  545.         final Map<String, Double> map = new HashMap<>();

  546.         for (int i = 0; i < CONSTANTS_MAX_NUMBER; ++i) {
  547.             // Note: for extracting the strings from the binary file, it makes no difference
  548.             //       if the file is stored in big-endian or little-endian notation
  549.             final String constantName = extractString(first, HEADER_CONSTANTS_NAMES_OFFSET + i * 6, 6);
  550.             if (constantName.length() == 0) {
  551.                 // no more constants to read
  552.                 break;
  553.             }
  554.             final double constantValue = extractDouble(second, HEADER_CONSTANTS_VALUES_OFFSET + 8 * i);
  555.             map.put(constantName, constantValue);
  556.         }

  557.         // INPOP files do not have constants for AU and EMRAT, thus extract them from
  558.         // the header record and create a constant for them to be consistent with JPL files
  559.         if (!map.containsKey(CONSTANT_AU)) {
  560.             map.put(CONSTANT_AU, extractDouble(first, HEADER_ASTRONOMICAL_UNIT_OFFSET));
  561.         }

  562.         if (!map.containsKey(CONSTANT_EMRAT)) {
  563.             map.put(CONSTANT_EMRAT, extractDouble(first, HEADER_EM_RATIO_OFFSET));
  564.         }

  565.         return map;

  566.     }

  567.     /** Read bytes into the current record array.
  568.      * @param input input stream
  569.      * @param record record where to put bytes
  570.      * @param start start index where to put bytes
  571.      * @return true if record has been filled up
  572.      * @exception IOException if a read error occurs
  573.      */
  574.     private boolean readInRecord(final InputStream input, final byte[] record, final int start)
  575.         throws IOException {
  576.         int index = start;
  577.         while (index != record.length) {
  578.             final int n = input.read(record, index, record.length - index);
  579.             if (n < 0) {
  580.                 return false;
  581.             }
  582.             index += n;
  583.         }
  584.         return true;
  585.     }

  586.     /** Detect whether the JPL ephemerides file is stored in big-endian or
  587.      * little-endian notation.
  588.      * @param record the array containing the binary JPL header
  589.      */
  590.     private void detectEndianess(final byte[] record) {

  591.         // default to big-endian
  592.         bigEndian = true;

  593.         // first try to read the DE number in big-endian format
  594.         // the number is stored as unsigned int, so we have to convert it properly
  595.         final long deNum = extractInt(record, HEADER_EPHEMERIS_TYPE_OFFSET) & 0xffffffffL;

  596.         // simple heuristic: if the read value is larger than half the range of an integer
  597.         //                   assume the file is in little-endian format
  598.         if (deNum > (1 << 15)) {
  599.             bigEndian = false;
  600.         }

  601.     }

  602.     /** Calculate the record size of a JPL ephemerides file.
  603.      * @param record the byte array containing the header record
  604.      * @param name the name of the data file
  605.      * @return the record size for this file
  606.      */
  607.     private int computeRecordSize(final byte[] record, final String name) {

  608.         int recordSize = 0;
  609.         boolean ok = true;
  610.         // JPL files always have 3 position components
  611.         final int nComp = 3;

  612.         // iterate over the coefficient ptr array and sum up the record size
  613.         // the coeffPtr array has the dimensions [12][nComp]
  614.         for (int j = 0; j < 12; j++) {
  615.             final int nCompCur = (j == 11) ? 2 : nComp;

  616.             // Note: the array element coeffPtr[j][0] is not needed for the calculation
  617.             final int idx = HEADER_CHEBISHEV_INDICES_OFFSET + j * nComp * 4;
  618.             final int coeffPtr1 = extractInt(record, idx + 4);
  619.             final int coeffPtr2 = extractInt(record, idx + 8);

  620.             // sanity checks
  621.             ok = ok && (coeffPtr1 >= 0 || coeffPtr2 >= 0);

  622.             recordSize += coeffPtr1 * coeffPtr2 * nCompCur;
  623.         }

  624.         // the libration ptr array has the dimension [3]
  625.         // Note: the array element libratPtr[0] is not needed for the calculation
  626.         final int libratPtr1 = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 4);
  627.         final int libratPtr2 = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 8);

  628.         // sanity checks
  629.         ok = ok && (libratPtr1 >= 0 || libratPtr2 >= 0);

  630.         recordSize += libratPtr1 * libratPtr2 * nComp + 2;
  631.         recordSize <<= 3;

  632.         if (!ok || recordSize <= 0) {
  633.             throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  634.         }

  635.         return recordSize;

  636.     }

  637.     /** Extract a date from a record.
  638.      * @param record record to parse
  639.      * @param offset offset of the double within the record
  640.      * @return extracted date
  641.      */
  642.     private AbsoluteDate extractDate(final byte[] record, final int offset) {

  643.         final double t = extractDouble(record, offset);
  644.         int    jDay    = (int) FastMath.floor(t);
  645.         double seconds = (t + 0.5 - jDay) * Constants.JULIAN_DAY;
  646.         if (seconds >= Constants.JULIAN_DAY) {
  647.             ++jDay;
  648.             seconds -= Constants.JULIAN_DAY;
  649.         }
  650.         return new AbsoluteDate(new DateComponents(DateComponents.JULIAN_EPOCH, jDay),
  651.                                 new TimeComponents(seconds), timeScale);
  652.     }

  653.     /** Extract a double from a record.
  654.      * <p>Double numbers are stored according to IEEE 754 standard, with
  655.      * most significant byte first.</p>
  656.      * @param record record to parse
  657.      * @param offset offset of the double within the record
  658.      * @return extracted double
  659.      */
  660.     private double extractDouble(final byte[] record, final int offset) {
  661.         final long l8 = ((long) record[offset + 0]) & 0xffl;
  662.         final long l7 = ((long) record[offset + 1]) & 0xffl;
  663.         final long l6 = ((long) record[offset + 2]) & 0xffl;
  664.         final long l5 = ((long) record[offset + 3]) & 0xffl;
  665.         final long l4 = ((long) record[offset + 4]) & 0xffl;
  666.         final long l3 = ((long) record[offset + 5]) & 0xffl;
  667.         final long l2 = ((long) record[offset + 6]) & 0xffl;
  668.         final long l1 = ((long) record[offset + 7]) & 0xffl;
  669.         final long l;
  670.         if (bigEndian) {
  671.             l = (l8 << 56) | (l7 << 48) | (l6 << 40) | (l5 << 32) |
  672.                 (l4 << 24) | (l3 << 16) | (l2 <<  8) | l1;
  673.         } else {
  674.             l = (l1 << 56) | (l2 << 48) | (l3 << 40) | (l4 << 32) |
  675.                 (l5 << 24) | (l6 << 16) | (l7 <<  8) | l8;
  676.         }
  677.         return Double.longBitsToDouble(l);
  678.     }

  679.     /** Extract an int from a record.
  680.      * @param record record to parse
  681.      * @param offset offset of the double within the record
  682.      * @return extracted int
  683.      */
  684.     private int extractInt(final byte[] record, final int offset) {
  685.         final int l4 = ((int) record[offset + 0]) & 0xff;
  686.         final int l3 = ((int) record[offset + 1]) & 0xff;
  687.         final int l2 = ((int) record[offset + 2]) & 0xff;
  688.         final int l1 = ((int) record[offset + 3]) & 0xff;

  689.         if (bigEndian) {
  690.             return (l4 << 24) | (l3 << 16) | (l2 <<  8) | l1;
  691.         } else {
  692.             return (l1 << 24) | (l2 << 16) | (l3 <<  8) | l4;
  693.         }
  694.     }

  695.     /** Extract a String from a record.
  696.      * @param record record to parse
  697.      * @param offset offset of the string within the record
  698.      * @param length maximal length of the string
  699.      * @return extracted string, with whitespace characters stripped
  700.      */
  701.     private String extractString(final byte[] record, final int offset, final int length) {
  702.         return new String(record, offset, length, StandardCharsets.US_ASCII).trim();
  703.     }

  704.     /** Local parser for header constants. */
  705.     private class ConstantsParser implements DataLoader {

  706.         /** Local constants map. */
  707.         private Map<String, Double> localConstants;

  708.        /** Get the local constants map.
  709.          * @return local constants map
  710.          */
  711.         public Map<String, Double> getConstants() {
  712.             return localConstants;
  713.         }

  714.         /** {@inheritDoc} */
  715.         public boolean stillAcceptsData() {
  716.             return localConstants == null;
  717.         }

  718.         /** {@inheritDoc} */
  719.         public void loadData(final InputStream input, final String name)
  720.             throws IOException, ParseException, OrekitException {

  721.             // read first header record
  722.             final byte[] first = readFirstRecord(input, name);

  723.             // the second record contains the values of the constants used for least-square filtering
  724.             final byte[] second = new byte[first.length];
  725.             if (!readInRecord(input, second, 0)) {
  726.                 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  727.             }

  728.             localConstants = parseConstants(first, second);

  729.         }

  730.     }

  731.     /** Local parser for Chebyshev polynomials. */
  732.     private class EphemerisParser implements DataLoader, TimeStampedGenerator<PosVelChebyshev> {

  733.         /** Set of Chebyshev polynomials read. */
  734.         private final SortedSet<PosVelChebyshev> entries;

  735.         /** Start of range we are interested in. */
  736.         private AbsoluteDate start;

  737.         /** End of range we are interested in. */
  738.         private AbsoluteDate end;

  739.         /** Simple constructor.
  740.          */
  741.         EphemerisParser() {
  742.             entries = new TreeSet<>(new ChronologicalComparator());
  743.         }

  744.         /** {@inheritDoc} */
  745.         public List<PosVelChebyshev> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {
  746.             try {

  747.                 // prepare reading
  748.                 entries.clear();
  749.                 if (existingDate == null) {
  750.                     // we want ephemeris data for the first time, set up an arbitrary first range
  751.                     start = date.shiftedBy(-FIFTY_DAYS);
  752.                     end   = date.shiftedBy(+FIFTY_DAYS);
  753.                 } else if (existingDate.compareTo(date) <= 0) {
  754.                     // we want to extend an existing range towards future dates
  755.                     start = existingDate;
  756.                     end   = date;
  757.                 } else {
  758.                     // we want to extend an existing range towards past dates
  759.                     start = date;
  760.                     end   = existingDate;
  761.                 }

  762.                 // get new entries in the specified data range
  763.                 if (!feed(this)) {
  764.                     throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
  765.                 }

  766.                 return new ArrayList<>(entries);

  767.             } catch (OrekitException oe) {
  768.                 throw new TimeStampedCacheException(oe);
  769.             }
  770.         }

  771.         /** {@inheritDoc} */
  772.         public boolean stillAcceptsData() {

  773.             // special case for Earth: we do not really load any ephemeris data
  774.             if (generateType == EphemerisType.EARTH) {
  775.                 return false;
  776.             }

  777.             // we have to look for data in all available ephemerides files as there may be
  778.             // data overlaps that result in incomplete data
  779.             if (entries.isEmpty()) {
  780.                 return true;
  781.             } else {
  782.                 // if the requested range is already filled, we do not need to look further
  783.                 return !(entries.first().getDate().compareTo(start) < 0 &&
  784.                          entries.last().getDate().compareTo(end)    > 0);
  785.             }

  786.         }

  787.         /** {@inheritDoc} */
  788.         public void loadData(final InputStream input, final String name)
  789.             throws IOException {

  790.             // read first header record
  791.             final byte[] first = readFirstRecord(input, name);

  792.             // the second record contains the values of the constants used for least-square filtering
  793.             final byte[] second = new byte[first.length];
  794.             if (!readInRecord(input, second, 0)) {
  795.                 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  796.             }

  797.             if (constants.get() == null) {
  798.                 constants.compareAndSet(null, parseConstants(first, second));
  799.             }

  800.             // check astronomical unit consistency
  801.             final double au = 1000 * extractDouble(first, HEADER_ASTRONOMICAL_UNIT_OFFSET);
  802.             if (au < 1.4e11 || au > 1.6e11) {
  803.                 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  804.             }
  805.             if (FastMath.abs(getLoadedAstronomicalUnit() - au) >= 10.0) {
  806.                 throw new OrekitException(OrekitMessages.INCONSISTENT_ASTRONOMICAL_UNIT_IN_FILES,
  807.                                           getLoadedAstronomicalUnit(), au);
  808.             }

  809.             // check Earth-Moon mass ratio consistency
  810.             final double emRat = extractDouble(first, HEADER_EM_RATIO_OFFSET);
  811.             if (emRat < 80 || emRat > 82) {
  812.                 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  813.             }
  814.             if (FastMath.abs(getLoadedEarthMoonMassRatio() - emRat) >= 1.0e-5) {
  815.                 throw new OrekitException(OrekitMessages.INCONSISTENT_EARTH_MOON_RATIO_IN_FILES,
  816.                                           getLoadedEarthMoonMassRatio(), emRat);
  817.             }

  818.             // parse first header record
  819.             parseFirstHeaderRecord(first, name);

  820.             if (startEpoch.compareTo(end) < 0 && finalEpoch.compareTo(start) > 0) {
  821.                 // this file contains data in the range we are looking for, read it
  822.                 final byte[] record = new byte[first.length];
  823.                 while (readInRecord(input, record, 0)) {
  824.                     final AbsoluteDate rangeStart = parseDataRecord(record);
  825.                     if (rangeStart.compareTo(end) > 0) {
  826.                         // we have already exceeded the range we were interested in,
  827.                         // we interrupt parsing here
  828.                         return;
  829.                     }
  830.                 }
  831.             }

  832.         }

  833.         /** Parse regular ephemeris record.
  834.          * @param record record to parse
  835.          * @return date of the last parsed chunk
  836.          */
  837.         private AbsoluteDate parseDataRecord(final byte[] record) {

  838.             // extract time range covered by the record
  839.             final AbsoluteDate rangeStart = extractDate(record, DATA_START_RANGE_OFFSET);
  840.             if (rangeStart.compareTo(startEpoch) < 0) {
  841.                 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE,
  842.                         rangeStart, startEpoch, finalEpoch, startEpoch.durationFrom(rangeStart));
  843.             }

  844.             final AbsoluteDate rangeEnd   = extractDate(record, DATE_END_RANGE_OFFSET);
  845.             if (rangeEnd.compareTo(finalEpoch) > 0) {
  846.                 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_AFTER,
  847.                         rangeEnd, startEpoch, finalEpoch, rangeEnd.durationFrom(finalEpoch));
  848.             }

  849.             if (rangeStart.compareTo(end) > 0 || rangeEnd.compareTo(start) < 0) {
  850.                 // we are not interested in this record, don't parse it
  851.                 return rangeEnd;
  852.             }

  853.             // loop over chunks inside the time range
  854.             AbsoluteDate chunkEnd = rangeStart;
  855.             final int nbChunks    = chunks;
  856.             final int nbCoeffs    = coeffs;
  857.             final int first       = firstIndex;
  858.             final double duration = chunksDuration;
  859.             for (int i = 0; i < nbChunks; ++i) {

  860.                 // set up chunk validity range
  861.                 final AbsoluteDate chunkStart = chunkEnd;
  862.                 chunkEnd = (i == nbChunks - 1) ? rangeEnd : rangeStart.shiftedBy((i + 1) * duration);

  863.                 // extract Chebyshev coefficients for the selected body
  864.                 // and convert them from kilometers to meters
  865.                 final double[] xCoeffs = new double[nbCoeffs];
  866.                 final double[] yCoeffs = new double[nbCoeffs];
  867.                 final double[] zCoeffs = new double[nbCoeffs];

  868.                 for (int k = 0; k < nbCoeffs; ++k) {
  869.                     // by now, only use the position components
  870.                     // if there are also velocity components contained in the file, ignore them
  871.                     final int index = first + components * i * nbCoeffs + k - 1;
  872.                     xCoeffs[k] = positionUnit * extractDouble(record, 8 * index);
  873.                     yCoeffs[k] = positionUnit * extractDouble(record, 8 * (index +  nbCoeffs));
  874.                     zCoeffs[k] = positionUnit * extractDouble(record, 8 * (index + 2 * nbCoeffs));
  875.                 }

  876.                 // build the position-velocity model for current chunk
  877.                 entries.add(new PosVelChebyshev(chunkStart, timeScale, duration, xCoeffs, yCoeffs, zCoeffs));

  878.             }

  879.             return rangeStart;

  880.         }

  881.     }

  882.     /** Raw position-velocity provider using ephemeris. */
  883.     private class EphemerisRawPVProvider implements RawPVProvider {


  884.         /** Retrieve correct Chebyshev polynomial for given date.
  885.          * @param date date
  886.          * @return PosVelChebyshev Chebyshev polynomial class for position-velocity-acceleration
  887.          */
  888.         private PosVelChebyshev getChebyshev(final AbsoluteDate date) {
  889.             PosVelChebyshev chebyshev;
  890.             try {
  891.                 chebyshev = ephemerides.getNeighbors(date).findFirst().get();
  892.             } catch (TimeStampedCacheException tce) {
  893.                 // we cannot bracket the date, check if the last available chunk covers the specified date
  894.                 chebyshev = ephemerides.getLatest();
  895.                 if (!chebyshev.inRange(date)) {
  896.                     // we were not able to recover from the error, the date is too far
  897.                     throw tce;
  898.                 }
  899.             }
  900.             return chebyshev;
  901.         }

  902.         /** {@inheritDoc} */
  903.         public PVCoordinates getRawPV(final AbsoluteDate date) {

  904.             // get raw PV from Chebyshev polynomials
  905.             final PosVelChebyshev chebyshev = getChebyshev(date);

  906.             // evaluate the Chebyshev polynomials
  907.             return chebyshev.getPositionVelocityAcceleration(date);

  908.         }

  909.         /** {@inheritDoc} */
  910.         public <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {

  911.             // get raw PV from Chebyshev polynomials
  912.             final PosVelChebyshev chebyshev = getChebyshev(date.toAbsoluteDate());

  913.             // evaluate the Chebyshev polynomials
  914.             return chebyshev.getPositionVelocityAcceleration(date);

  915.         }

  916.         /** {@inheritDoc} */
  917.         @Override
  918.         public Vector3D getRawPosition(final AbsoluteDate date) {
  919.             // get raw PV from Chebyshev polynomials
  920.             final PosVelChebyshev chebyshev = getChebyshev(date);

  921.             // evaluate the Chebyshev polynomials
  922.             return chebyshev.getPosition(date);
  923.         }

  924.         /** {@inheritDoc} */
  925.         @Override
  926.         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getRawPosition(final FieldAbsoluteDate<T> date) {
  927.             // get raw PV from Chebyshev polynomials
  928.             final PosVelChebyshev chebyshev = getChebyshev(date.toAbsoluteDate());

  929.             // evaluate the Chebyshev polynomials
  930.             return chebyshev.getPosition(date);
  931.         }

  932.     }

  933.     /** Raw position-velocity provider providing always zero. */
  934.     private static class ZeroRawPVProvider implements RawPVProvider {

  935.         /** {@inheritDoc} */
  936.         public PVCoordinates getRawPV(final AbsoluteDate date) {
  937.             return PVCoordinates.ZERO;
  938.         }

  939.         /** {@inheritDoc} */
  940.         public <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {
  941.             return FieldPVCoordinates.getZero(date.getField());
  942.         }

  943.     }

  944. }