JPLEphemeridesLoader.java

  1. /* Copyright 2002-2020 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.bodies;

  18. 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.RealFieldElement;
  30. import org.hipparchus.util.FastMath;
  31. import org.orekit.annotation.DefaultDataContext;
  32. import org.orekit.data.AbstractSelfFeedingLoader;
  33. import org.orekit.data.DataContext;
  34. import org.orekit.data.DataLoader;
  35. import org.orekit.data.DataProvidersManager;
  36. import org.orekit.errors.OrekitException;
  37. import org.orekit.errors.OrekitInternalError;
  38. import org.orekit.errors.OrekitMessages;
  39. import org.orekit.errors.TimeStampedCacheException;
  40. import org.orekit.frames.Frame;
  41. import org.orekit.frames.Predefined;
  42. import org.orekit.time.AbsoluteDate;
  43. import org.orekit.time.ChronologicalComparator;
  44. import org.orekit.time.DateComponents;
  45. import org.orekit.time.FieldAbsoluteDate;
  46. import org.orekit.time.TimeComponents;
  47. import org.orekit.time.TimeScale;
  48. import org.orekit.time.TimeScales;
  49. import org.orekit.utils.Constants;
  50. import org.orekit.utils.FieldPVCoordinates;
  51. import org.orekit.utils.OrekitConfiguration;
  52. import org.orekit.utils.PVCoordinates;
  53. import org.orekit.utils.GenericTimeStampedCache;
  54. import org.orekit.utils.TimeStampedGenerator;

  55. /** Loader for JPL ephemerides binary files (DE 4xx) and similar formats (INPOP 06/08/10).
  56.  * <p>JPL ephemerides binary files contain ephemerides for all solar system planets.</p>
  57.  * <p>The JPL ephemerides binary files are recognized thanks to their base names,
  58.  * which must match the pattern <code>[lu]nx[mp]####.ddd</code> (or
  59.  * <code>[lu]nx[mp]####.ddd.gz</code> for gzip-compressed files) where # stands for a
  60.  * digit character and where ddd is an ephemeris type (typically 405 or 406).</p>
  61.  * <p>The loader supports files encoded in big-endian as well as in little-endian notation.
  62.  * Usually, big-endian files are named <code>unx[mp]####.ddd</code>, while little-endian files
  63.  * are named <code>lnx[mp]####.ddd</code>.</p>
  64.  * <p>The IMCCE ephemerides binary files are recognized thanks to their base names,
  65.  * which must match the pattern <code>inpop*.dat</code> (or
  66.  * <code>inpop*.dat.gz</code> for gzip-compressed files) where * stands for any string.</p>
  67.  * <p>The loader supports files encoded in big-endian as well as in little-endian notation.
  68.  * Usually, big-endian files contain <code>bigendian</code> in their names, while little-endian files
  69.  * contain <code>littleendian</code> in their names.</p>
  70.  * <p>The loader supports files in TDB or TCB time scales.</p>
  71.  * @author Luc Maisonobe
  72.  */
  73. public class JPLEphemeridesLoader extends AbstractSelfFeedingLoader
  74.         implements CelestialBodyLoader {

  75.     /** Default supported files name pattern for JPL DE files. */
  76.     public static final String DEFAULT_DE_SUPPORTED_NAMES = "^[lu]nx([mp](\\d\\d\\d\\d))+\\.(?:4\\d\\d)$";

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  115.     /** List of supported ephemerides types. */
  116.     public enum EphemerisType {

  117.         /** Constant for solar system barycenter. */
  118.         SOLAR_SYSTEM_BARYCENTER,

  119.         /** Constant for the Sun. */
  120.         SUN,

  121.         /** Constant for Mercury. */
  122.         MERCURY,

  123.         /** Constant for Venus. */
  124.         VENUS,

  125.         /** Constant for the Earth-Moon barycenter. */
  126.         EARTH_MOON,

  127.         /** Constant for the Earth. */
  128.         EARTH,

  129.         /** Constant for the Moon. */
  130.         MOON,

  131.         /** Constant for Mars. */
  132.         MARS,

  133.         /** Constant for Jupiter. */
  134.         JUPITER,

  135.         /** Constant for Saturn. */
  136.         SATURN,

  137.         /** Constant for Uranus. */
  138.         URANUS,

  139.         /** Constant for Neptune. */
  140.         NEPTUNE,

  141.         /** Constant for Pluto. */
  142.         PLUTO

  143.     }

  144.     /** Interface for raw position-velocity retrieval. */
  145.     public interface RawPVProvider {

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

  151.         /** Get the position-velocity at date.
  152.          * @param date date at which the position-velocity is desired
  153.          * @param <T> type of the field elements
  154.          * @return position-velocity at the specified date
  155.          */
  156.         <T extends RealFieldElement<T>> FieldPVCoordinates<T> getRawPV(FieldAbsoluteDate<T> date);

  157.     }

  158.     /** Ephemeris for selected body. */
  159.     private final GenericTimeStampedCache<PosVelChebyshev> ephemerides;

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

  162.     /** Ephemeris type to generate. */
  163.     private final EphemerisType generateType;

  164.     /** Ephemeris type to load. */
  165.     private final EphemerisType loadType;

  166.     /** Time scales to use when loading data. */
  167.     private final TimeScales timeScales;

  168.     /** The GCRF implementation. */
  169.     private final Frame gcrf;

  170.     /** Current file start epoch. */
  171.     private AbsoluteDate startEpoch;

  172.     /** Current file final epoch. */
  173.     private AbsoluteDate finalEpoch;

  174.     /** Chunks duration (in seconds). */
  175.     private double maxChunksDuration;

  176.     /** Current file chunks duration (in seconds). */
  177.     private double chunksDuration;

  178.     /** Index of the first data for selected body. */
  179.     private int firstIndex;

  180.     /** Number of coefficients for selected body. */
  181.     private int coeffs;

  182.     /** Number of chunks for the selected body. */
  183.     private int chunks;

  184.     /** Number of components contained in the file. */
  185.     private int components;

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

  188.     /** Time scale of the date coordinates. */
  189.     private TimeScale timeScale;

  190.     /** Indicator for binary file endianness. */
  191.     private boolean bigEndian;

  192.     /** Create a loader for JPL ephemerides binary files. This constructor uses the {@link
  193.      * DataContext#getDefault() default data context}.
  194.      *
  195.      * @param supportedNames regular expression for supported files names
  196.      * @param generateType ephemeris type to generate
  197.      * @see #JPLEphemeridesLoader(String, EphemerisType, DataProvidersManager, TimeScales,
  198.      * Frame)
  199.      */
  200.     @DefaultDataContext
  201.     public JPLEphemeridesLoader(final String supportedNames, final EphemerisType generateType) {
  202.         this(supportedNames, generateType,
  203.                 DataContext.getDefault().getDataProvidersManager(),
  204.                 DataContext.getDefault().getTimeScales(),
  205.                 DataContext.getDefault().getFrames().getGCRF());
  206.     }

  207.     /** Create a loader for JPL ephemerides binary files.
  208.      * @param supportedNames regular expression for supported files names
  209.      * @param generateType ephemeris type to generate
  210.      * @param dataProvidersManager provides access to the ephemeris files.
  211.      * @param timeScales used to access the TCB and TDB time scales while loading data.
  212.      * @param gcrf Earth centered frame aligned with ICRF.
  213.      * @since 10.1
  214.      */
  215.     public JPLEphemeridesLoader(final String supportedNames,
  216.                                 final EphemerisType generateType,
  217.                                 final DataProvidersManager dataProvidersManager,
  218.                                 final TimeScales timeScales,
  219.                                 final Frame gcrf) {
  220.         super(supportedNames, dataProvidersManager);

  221.         this.timeScales = timeScales;
  222.         this.gcrf = gcrf;
  223.         constants = new AtomicReference<>();

  224.         this.generateType  = generateType;
  225.         if (generateType == EphemerisType.SOLAR_SYSTEM_BARYCENTER) {
  226.             loadType = EphemerisType.EARTH_MOON;
  227.         } else if (generateType == EphemerisType.EARTH_MOON) {
  228.             loadType = EphemerisType.MOON;
  229.         } else {
  230.             loadType = generateType;
  231.         }

  232.         ephemerides = new GenericTimeStampedCache<>(
  233.                 2, OrekitConfiguration.getCacheSlotsNumber(),
  234.                 Double.POSITIVE_INFINITY, FIFTY_DAYS,
  235.                 new EphemerisParser());
  236.         maxChunksDuration = Double.NaN;
  237.         chunksDuration    = Double.NaN;

  238.     }

  239.     /** Load celestial body.
  240.      * @param name name of the celestial body
  241.      * @return loaded celestial body
  242.      */
  243.     public CelestialBody loadCelestialBody(final String name) {

  244.         final double gm       = getLoadedGravitationalCoefficient(generateType);
  245.         final IAUPole iauPole = PredefinedIAUPoles
  246.                 .getIAUPole(generateType, timeScales);
  247.         final double scale;
  248.         final Frame definingFrameAlignedWithICRF;
  249.         final RawPVProvider rawPVProvider;
  250.         String inertialFrameName = null;
  251.         String bodyOrientedFrameName = null;
  252.         switch (generateType) {
  253.             case SOLAR_SYSTEM_BARYCENTER : {
  254.                 scale = -1.0;
  255.                 final JPLEphemeridesLoader parentLoader = new JPLEphemeridesLoader(
  256.                         getSupportedNames(),
  257.                         EphemerisType.EARTH_MOON,
  258.                         getDataProvidersManager(),
  259.                         timeScales,
  260.                         gcrf);
  261.                 final CelestialBody parentBody =
  262.                         parentLoader.loadCelestialBody(CelestialBodyFactory.EARTH_MOON);
  263.                 definingFrameAlignedWithICRF = parentBody.getInertiallyOrientedFrame();
  264.                 rawPVProvider = new EphemerisRawPVProvider();
  265.                 inertialFrameName = Predefined.ICRF.getName();
  266.                 bodyOrientedFrameName = null;
  267.                 break;
  268.             }
  269.             case EARTH_MOON :
  270.                 scale         = 1.0 / (1.0 + getLoadedEarthMoonMassRatio());
  271.                 definingFrameAlignedWithICRF = gcrf;
  272.                 rawPVProvider = new EphemerisRawPVProvider();
  273.                 break;
  274.             case EARTH :
  275.                 scale         = 1.0;
  276.                 definingFrameAlignedWithICRF = gcrf;
  277.                 rawPVProvider = new ZeroRawPVProvider();
  278.                 break;
  279.             case MOON :
  280.                 scale         =  1.0;
  281.                 definingFrameAlignedWithICRF = gcrf;
  282.                 rawPVProvider = new EphemerisRawPVProvider();
  283.                 break;
  284.             default : {
  285.                 scale = 1.0;
  286.                 final JPLEphemeridesLoader parentLoader = new JPLEphemeridesLoader(
  287.                         getSupportedNames(),
  288.                         EphemerisType.SOLAR_SYSTEM_BARYCENTER,
  289.                         getDataProvidersManager(),
  290.                         timeScales,
  291.                         gcrf);
  292.                 final CelestialBody parentBody =
  293.                         parentLoader.loadCelestialBody(CelestialBodyFactory.SOLAR_SYSTEM_BARYCENTER);
  294.                 definingFrameAlignedWithICRF = parentBody.getInertiallyOrientedFrame();
  295.                 rawPVProvider = new EphemerisRawPVProvider();
  296.             }
  297.         }

  298.         // build the celestial body
  299.         return new JPLCelestialBody(name, getSupportedNames(), generateType, rawPVProvider,
  300.                                     gm, scale, iauPole, definingFrameAlignedWithICRF,
  301.                                     inertialFrameName, bodyOrientedFrameName);

  302.     }

  303.     /** Get astronomical unit.
  304.      * @return astronomical unit in meters
  305.      */
  306.     public double getLoadedAstronomicalUnit() {
  307.         return 1000.0 * getLoadedConstant(CONSTANT_AU);
  308.     }

  309.     /** Get Earth/Moon mass ratio.
  310.      * @return Earth/Moon mass ratio
  311.      */
  312.     public double getLoadedEarthMoonMassRatio() {
  313.         return getLoadedConstant(CONSTANT_EMRAT);
  314.     }

  315.     /** Get the gravitational coefficient of a body.
  316.      * @param body body for which the gravitational coefficient is requested
  317.      * @return gravitational coefficient in m³/s²
  318.      */
  319.     public double getLoadedGravitationalCoefficient(final EphemerisType body) {

  320.         // coefficient in au³/day²
  321.         final double rawGM;
  322.         switch (body) {
  323.             case SOLAR_SYSTEM_BARYCENTER :
  324.                 return getLoadedGravitationalCoefficient(EphemerisType.SUN)        +
  325.                         getLoadedGravitationalCoefficient(EphemerisType.MERCURY)    +
  326.                         getLoadedGravitationalCoefficient(EphemerisType.VENUS)      +
  327.                         getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) +
  328.                         getLoadedGravitationalCoefficient(EphemerisType.MARS)       +
  329.                         getLoadedGravitationalCoefficient(EphemerisType.JUPITER)    +
  330.                         getLoadedGravitationalCoefficient(EphemerisType.SATURN)     +
  331.                         getLoadedGravitationalCoefficient(EphemerisType.URANUS)     +
  332.                         getLoadedGravitationalCoefficient(EphemerisType.NEPTUNE)    +
  333.                         getLoadedGravitationalCoefficient(EphemerisType.PLUTO);
  334.             case SUN :
  335.                 rawGM = getLoadedConstant("GMS", "GM_Sun");
  336.                 break;
  337.             case MERCURY :
  338.                 rawGM = getLoadedConstant("GM1", "GM_Mer");
  339.                 break;
  340.             case VENUS :
  341.                 rawGM = getLoadedConstant("GM2", "GM_Ven");
  342.                 break;
  343.             case EARTH_MOON :
  344.                 rawGM = getLoadedConstant("GMB", "GM_EMB");
  345.                 break;
  346.             case EARTH :
  347.                 return getLoadedEarthMoonMassRatio() *
  348.                         getLoadedGravitationalCoefficient(EphemerisType.MOON);
  349.             case MOON :
  350.                 return getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) /
  351.                         (1.0 + getLoadedEarthMoonMassRatio());
  352.             case MARS :
  353.                 rawGM = getLoadedConstant("GM4", "GM_Mar");
  354.                 break;
  355.             case JUPITER :
  356.                 rawGM = getLoadedConstant("GM5", "GM_Jup");
  357.                 break;
  358.             case SATURN :
  359.                 rawGM = getLoadedConstant("GM6", "GM_Sat");
  360.                 break;
  361.             case URANUS :
  362.                 rawGM = getLoadedConstant("GM7", "GM_Ura");
  363.                 break;
  364.             case NEPTUNE :
  365.                 rawGM = getLoadedConstant("GM8", "GM_Nep");
  366.                 break;
  367.             case PLUTO :
  368.                 rawGM = getLoadedConstant("GM9", "GM_Plu");
  369.                 break;
  370.             default :
  371.                 throw new OrekitInternalError(null);
  372.         }

  373.         final double au    = getLoadedAstronomicalUnit();
  374.         return rawGM * au * au * au / (Constants.JULIAN_DAY * Constants.JULIAN_DAY);

  375.     }

  376.     /** Get a constant defined in the ephemerides headers.
  377.      * <p>Note that since constants are defined in the JPL headers
  378.      * files, they are available as soon as one file is available, even
  379.      * if it doesn't match the desired central date. This is because the
  380.      * header must be parsed before the dates can be checked.</p>
  381.      * <p>
  382.      * There are alternate names for constants since for example JPL names are
  383.      * different from INPOP names (Sun gravity: GMS or GM_Sun, Mars gravity:
  384.      * GM4 or GM_Mar...).
  385.      * </p>
  386.      * @param names alternate names of the constant
  387.      * @return value of the constant of NaN if the constant is not defined
  388.      */
  389.     public double getLoadedConstant(final String... names) {

  390.         // lazy loading of constants
  391.         Map<String, Double> map = constants.get();
  392.         if (map == null) {
  393.             final ConstantsParser parser = new ConstantsParser();
  394.             if (!feed(parser)) {
  395.                 throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
  396.             }
  397.             map = parser.getConstants();
  398.             constants.compareAndSet(null, map);
  399.         }

  400.         for (final String name : names) {
  401.             if (map.containsKey(name)) {
  402.                 return map.get(name);
  403.             }
  404.         }

  405.         return Double.NaN;

  406.     }

  407.     /** Get the maximal chunks duration.
  408.      * @return chunks maximal duration in seconds
  409.      */
  410.     public double getMaxChunksDuration() {
  411.         return maxChunksDuration;
  412.     }

  413.     /** Parse the first header record.
  414.      * @param record first header record
  415.      * @param name name of the file (or zip entry)
  416.      */
  417.     private void parseFirstHeaderRecord(final byte[] record, final String name) {

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

  420.         // as default, 3 polynomial coefficients for the Cartesian coordinates
  421.         // (x, y, z) are contained in the file, positions are in kilometers
  422.         // and times are in TDB
  423.         components   = 3;
  424.         positionUnit = 1000.0;
  425.         timeScale    = timeScales.getTDB();

  426.         if (deNum == INPOP_DE_NUMBER) {
  427.             // an INPOP file may contain 6 components (including coefficients for the velocity vector)
  428.             final double format = getLoadedConstant("FORMAT");
  429.             if (!Double.isNaN(format) && (int) FastMath.IEEEremainder(format, 10) != 1) {
  430.                 components = 6;
  431.             }

  432.             // INPOP files may have their polynomials expressed in AU
  433.             final double unite = getLoadedConstant("UNITE");
  434.             if (!Double.isNaN(unite) && (int) unite == 0) {
  435.                 positionUnit = getLoadedAstronomicalUnit();
  436.             }

  437.             // INPOP files may have their times expressed in TCB
  438.             final double timesc = getLoadedConstant("TIMESC");
  439.             if (!Double.isNaN(timesc) && (int) timesc == 1) {
  440.                 timeScale = timeScales.getTCB();
  441.             }

  442.         }

  443.         // extract covered date range
  444.         startEpoch = extractDate(record, HEADER_START_EPOCH_OFFSET);
  445.         finalEpoch = extractDate(record, HEADER_END_EPOCH_OFFSET);
  446.         boolean ok = finalEpoch.compareTo(startEpoch) > 0;

  447.         // indices of the Chebyshev coefficients for each ephemeris
  448.         for (int i = 0; i < 12; ++i) {
  449.             final int row1 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET     + 12 * i);
  450.             final int row2 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 4 + 12 * i);
  451.             final int row3 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 8 + 12 * i);
  452.             ok = ok && (row1 >= 0) && (row2 >= 0) && (row3 >= 0);
  453.             if (((i ==  0) && (loadType == EphemerisType.MERCURY))    ||
  454.                     ((i ==  1) && (loadType == EphemerisType.VENUS))      ||
  455.                     ((i ==  2) && (loadType == EphemerisType.EARTH_MOON)) ||
  456.                     ((i ==  3) && (loadType == EphemerisType.MARS))       ||
  457.                     ((i ==  4) && (loadType == EphemerisType.JUPITER))    ||
  458.                     ((i ==  5) && (loadType == EphemerisType.SATURN))     ||
  459.                     ((i ==  6) && (loadType == EphemerisType.URANUS))     ||
  460.                     ((i ==  7) && (loadType == EphemerisType.NEPTUNE))    ||
  461.                     ((i ==  8) && (loadType == EphemerisType.PLUTO))      ||
  462.                     ((i ==  9) && (loadType == EphemerisType.MOON))       ||
  463.                     ((i == 10) && (loadType == EphemerisType.SUN))) {
  464.                 firstIndex = row1;
  465.                 coeffs     = row2;
  466.                 chunks     = row3;
  467.             }
  468.         }

  469.         // compute chunks duration
  470.         final double timeSpan = extractDouble(record, HEADER_CHUNK_DURATION_OFFSET);
  471.         ok = ok && (timeSpan > 0) && (timeSpan < 100);
  472.         chunksDuration = Constants.JULIAN_DAY * (timeSpan / chunks);
  473.         if (Double.isNaN(maxChunksDuration)) {
  474.             maxChunksDuration = chunksDuration;
  475.         } else {
  476.             maxChunksDuration = FastMath.max(maxChunksDuration, chunksDuration);
  477.         }

  478.         // sanity checks
  479.         if (!ok) {
  480.             throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  481.         }

  482.     }

  483.     /** Read first header record.
  484.      * @param input input stream
  485.      * @param name name of the file (or zip entry)
  486.      * @return record record where to put bytes
  487.      * @exception IOException if a read error occurs
  488.      */
  489.     private byte[] readFirstRecord(final InputStream input, final String name)
  490.         throws IOException {

  491.         // read first part of record, up to the ephemeris type
  492.         final byte[] firstPart = new byte[HEADER_RECORD_SIZE_OFFSET + 4];
  493.         if (!readInRecord(input, firstPart, 0)) {
  494.             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  495.         }

  496.         // detect the endian format
  497.         detectEndianess(firstPart);

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

  500.         // the record size for this file
  501.         final int recordSize;

  502.         if (deNum == INPOP_DE_NUMBER) {
  503.             // INPOP files have an extended DE format, which includes also the record size
  504.             recordSize = extractInt(firstPart, HEADER_RECORD_SIZE_OFFSET) << 3;
  505.         } else {
  506.             // compute the record size for original JPL files
  507.             recordSize = computeRecordSize(firstPart, name);
  508.         }

  509.         if (recordSize <= 0) {
  510.             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  511.         }

  512.         // build a record with the proper size and finish read of the first complete record
  513.         final int start = firstPart.length;
  514.         final byte[] record = new byte[recordSize];
  515.         System.arraycopy(firstPart, 0, record, 0, firstPart.length);
  516.         if (!readInRecord(input, record, start)) {
  517.             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  518.         }

  519.         return record;

  520.     }

  521.     /** Parse constants from first two header records.
  522.      * @param first first header record
  523.      * @param second second header record
  524.      * @return map of parsed constants
  525.      */
  526.     private Map<String, Double> parseConstants(final byte[] first, final byte[] second) {

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

  528.         for (int i = 0; i < CONSTANTS_MAX_NUMBER; ++i) {
  529.             // Note: for extracting the strings from the binary file, it makes no difference
  530.             //       if the file is stored in big-endian or little-endian notation
  531.             final String constantName = extractString(first, HEADER_CONSTANTS_NAMES_OFFSET + i * 6, 6);
  532.             if (constantName.length() == 0) {
  533.                 // no more constants to read
  534.                 break;
  535.             }
  536.             final double constantValue = extractDouble(second, HEADER_CONSTANTS_VALUES_OFFSET + 8 * i);
  537.             map.put(constantName, constantValue);
  538.         }

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

  544.         if (!map.containsKey(CONSTANT_EMRAT)) {
  545.             map.put(CONSTANT_EMRAT, extractDouble(first, HEADER_EM_RATIO_OFFSET));
  546.         }

  547.         return map;

  548.     }

  549.     /** Read bytes into the current record array.
  550.      * @param input input stream
  551.      * @param record record where to put bytes
  552.      * @param start start index where to put bytes
  553.      * @return true if record has been filled up
  554.      * @exception IOException if a read error occurs
  555.      */
  556.     private boolean readInRecord(final InputStream input, final byte[] record, final int start)
  557.         throws IOException {
  558.         int index = start;
  559.         while (index != record.length) {
  560.             final int n = input.read(record, index, record.length - index);
  561.             if (n < 0) {
  562.                 return false;
  563.             }
  564.             index += n;
  565.         }
  566.         return true;
  567.     }

  568.     /** Detect whether the JPL ephemerides file is stored in big-endian or
  569.      * little-endian notation.
  570.      * @param record the array containing the binary JPL header
  571.      */
  572.     private void detectEndianess(final byte[] record) {

  573.         // default to big-endian
  574.         bigEndian = true;

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

  578.         // simple heuristic: if the read value is larger than half the range of an integer
  579.         //                   assume the file is in little-endian format
  580.         if (deNum > (1 << 15)) {
  581.             bigEndian = false;
  582.         }

  583.     }

  584.     /** Calculate the record size of a JPL ephemerides file.
  585.      * @param record the byte array containing the header record
  586.      * @param name the name of the data file
  587.      * @return the record size for this file
  588.      */
  589.     private int computeRecordSize(final byte[] record, final String name) {

  590.         int recordSize = 0;
  591.         boolean ok = true;
  592.         // JPL files always have 3 position components
  593.         final int nComp = 3;

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

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

  602.             // sanity checks
  603.             ok = ok && (coeffPtr1 >= 0 || coeffPtr2 >= 0);

  604.             recordSize += coeffPtr1 * coeffPtr2 * nCompCur;
  605.         }

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

  610.         // sanity checks
  611.         ok = ok && (libratPtr1 >= 0 || libratPtr2 >= 0);

  612.         recordSize += libratPtr1 * libratPtr2 * nComp + 2;
  613.         recordSize <<= 3;

  614.         if (!ok || recordSize <= 0) {
  615.             throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  616.         }

  617.         return recordSize;

  618.     }

  619.     /** Extract a date from a record.
  620.      * @param record record to parse
  621.      * @param offset offset of the double within the record
  622.      * @return extracted date
  623.      */
  624.     private AbsoluteDate extractDate(final byte[] record, final int offset) {

  625.         final double t = extractDouble(record, offset);
  626.         int    jDay    = (int) FastMath.floor(t);
  627.         double seconds = (t + 0.5 - jDay) * Constants.JULIAN_DAY;
  628.         if (seconds >= Constants.JULIAN_DAY) {
  629.             ++jDay;
  630.             seconds -= Constants.JULIAN_DAY;
  631.         }
  632.         return new AbsoluteDate(new DateComponents(DateComponents.JULIAN_EPOCH, jDay),
  633.                                 new TimeComponents(seconds), timeScale);
  634.     }

  635.     /** Extract a double from a record.
  636.      * <p>Double numbers are stored according to IEEE 754 standard, with
  637.      * most significant byte first.</p>
  638.      * @param record record to parse
  639.      * @param offset offset of the double within the record
  640.      * @return extracted double
  641.      */
  642.     private double extractDouble(final byte[] record, final int offset) {
  643.         final long l8 = ((long) record[offset + 0]) & 0xffl;
  644.         final long l7 = ((long) record[offset + 1]) & 0xffl;
  645.         final long l6 = ((long) record[offset + 2]) & 0xffl;
  646.         final long l5 = ((long) record[offset + 3]) & 0xffl;
  647.         final long l4 = ((long) record[offset + 4]) & 0xffl;
  648.         final long l3 = ((long) record[offset + 5]) & 0xffl;
  649.         final long l2 = ((long) record[offset + 6]) & 0xffl;
  650.         final long l1 = ((long) record[offset + 7]) & 0xffl;
  651.         final long l;
  652.         if (bigEndian) {
  653.             l = (l8 << 56) | (l7 << 48) | (l6 << 40) | (l5 << 32) |
  654.                 (l4 << 24) | (l3 << 16) | (l2 <<  8) | l1;
  655.         } else {
  656.             l = (l1 << 56) | (l2 << 48) | (l3 << 40) | (l4 << 32) |
  657.                 (l5 << 24) | (l6 << 16) | (l7 <<  8) | l8;
  658.         }
  659.         return Double.longBitsToDouble(l);
  660.     }

  661.     /** Extract an int from a record.
  662.      * @param record record to parse
  663.      * @param offset offset of the double within the record
  664.      * @return extracted int
  665.      */
  666.     private int extractInt(final byte[] record, final int offset) {
  667.         final int l4 = ((int) record[offset + 0]) & 0xff;
  668.         final int l3 = ((int) record[offset + 1]) & 0xff;
  669.         final int l2 = ((int) record[offset + 2]) & 0xff;
  670.         final int l1 = ((int) record[offset + 3]) & 0xff;

  671.         if (bigEndian) {
  672.             return (l4 << 24) | (l3 << 16) | (l2 <<  8) | l1;
  673.         } else {
  674.             return (l1 << 24) | (l2 << 16) | (l3 <<  8) | l4;
  675.         }
  676.     }

  677.     /** Extract a String from a record.
  678.      * @param record record to parse
  679.      * @param offset offset of the string within the record
  680.      * @param length maximal length of the string
  681.      * @return extracted string, with whitespace characters stripped
  682.      */
  683.     private String extractString(final byte[] record, final int offset, final int length) {
  684.         return new String(record, offset, length, StandardCharsets.US_ASCII).trim();
  685.     }

  686.     /** Local parser for header constants. */
  687.     private class ConstantsParser implements DataLoader {

  688.         /** Local constants map. */
  689.         private Map<String, Double> localConstants;

  690.        /** Get the local constants map.
  691.          * @return local constants map
  692.          */
  693.         public Map<String, Double> getConstants() {
  694.             return localConstants;
  695.         }

  696.         /** {@inheritDoc} */
  697.         public boolean stillAcceptsData() {
  698.             return localConstants == null;
  699.         }

  700.         /** {@inheritDoc} */
  701.         public void loadData(final InputStream input, final String name)
  702.             throws IOException, ParseException, OrekitException {

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

  705.             // the second record contains the values of the constants used for least-square filtering
  706.             final byte[] second = new byte[first.length];
  707.             if (!readInRecord(input, second, 0)) {
  708.                 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  709.             }

  710.             localConstants = parseConstants(first, second);

  711.         }

  712.     }

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

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

  717.         /** Start of range we are interested in. */
  718.         private AbsoluteDate start;

  719.         /** End of range we are interested in. */
  720.         private AbsoluteDate end;

  721.         /** Simple constructor.
  722.          */
  723.         EphemerisParser() {
  724.             entries = new TreeSet<>(new ChronologicalComparator());
  725.         }

  726.         /** {@inheritDoc} */
  727.         public List<PosVelChebyshev> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {
  728.             try {

  729.                 // prepare reading
  730.                 entries.clear();
  731.                 if (existingDate == null) {
  732.                     // we want ephemeris data for the first time, set up an arbitrary first range
  733.                     start = date.shiftedBy(-FIFTY_DAYS);
  734.                     end   = date.shiftedBy(+FIFTY_DAYS);
  735.                 } else if (existingDate.compareTo(date) <= 0) {
  736.                     // we want to extend an existing range towards future dates
  737.                     start = existingDate;
  738.                     end   = date;
  739.                 } else {
  740.                     // we want to extend an existing range towards past dates
  741.                     start = date;
  742.                     end   = existingDate;
  743.                 }

  744.                 // get new entries in the specified data range
  745.                 if (!feed(this)) {
  746.                     throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
  747.                 }

  748.                 return new ArrayList<>(entries);

  749.             } catch (OrekitException oe) {
  750.                 throw new TimeStampedCacheException(oe);
  751.             }
  752.         }

  753.         /** {@inheritDoc} */
  754.         public boolean stillAcceptsData() {

  755.             // special case for Earth: we do not really load any ephemeris data
  756.             if (generateType == EphemerisType.EARTH) {
  757.                 return false;
  758.             }

  759.             // we have to look for data in all available ephemerides files as there may be
  760.             // data overlaps that result in incomplete data
  761.             if (entries.isEmpty()) {
  762.                 return true;
  763.             } else {
  764.                 // if the requested range is already filled, we do not need to look further
  765.                 return !(entries.first().getDate().compareTo(start) < 0 &&
  766.                          entries.last().getDate().compareTo(end)    > 0);
  767.             }

  768.         }

  769.         /** {@inheritDoc} */
  770.         public void loadData(final InputStream input, final String name)
  771.             throws IOException {

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

  774.             // the second record contains the values of the constants used for least-square filtering
  775.             final byte[] second = new byte[first.length];
  776.             if (!readInRecord(input, second, 0)) {
  777.                 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  778.             }

  779.             if (constants.get() == null) {
  780.                 constants.compareAndSet(null, parseConstants(first, second));
  781.             }

  782.             // check astronomical unit consistency
  783.             final double au = 1000 * extractDouble(first, HEADER_ASTRONOMICAL_UNIT_OFFSET);
  784.             if ((au < 1.4e11) || (au > 1.6e11)) {
  785.                 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  786.             }
  787.             if (FastMath.abs(getLoadedAstronomicalUnit() - au) >= 10.0) {
  788.                 throw new OrekitException(OrekitMessages.INCONSISTENT_ASTRONOMICAL_UNIT_IN_FILES,
  789.                                           getLoadedAstronomicalUnit(), au);
  790.             }

  791.             // check Earth-Moon mass ratio consistency
  792.             final double emRat = extractDouble(first, HEADER_EM_RATIO_OFFSET);
  793.             if ((emRat < 80) || (emRat > 82)) {
  794.                 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  795.             }
  796.             if (FastMath.abs(getLoadedEarthMoonMassRatio() - emRat) >= 1.0e-5) {
  797.                 throw new OrekitException(OrekitMessages.INCONSISTENT_EARTH_MOON_RATIO_IN_FILES,
  798.                                           getLoadedEarthMoonMassRatio(), emRat);
  799.             }

  800.             // parse first header record
  801.             parseFirstHeaderRecord(first, name);

  802.             if (startEpoch.compareTo(end) < 0 && finalEpoch.compareTo(start) > 0) {
  803.                 // this file contains data in the range we are looking for, read it
  804.                 final byte[] record = new byte[first.length];
  805.                 while (readInRecord(input, record, 0)) {
  806.                     final AbsoluteDate rangeStart = parseDataRecord(record);
  807.                     if (rangeStart.compareTo(end) > 0) {
  808.                         // we have already exceeded the range we were interested in,
  809.                         // we interrupt parsing here
  810.                         return;
  811.                     }
  812.                 }
  813.             }

  814.         }

  815.         /** Parse regular ephemeris record.
  816.          * @param record record to parse
  817.          * @return date of the last parsed chunk
  818.          */
  819.         private AbsoluteDate parseDataRecord(final byte[] record) {

  820.             // extract time range covered by the record
  821.             final AbsoluteDate rangeStart = extractDate(record, DATA_START_RANGE_OFFSET);
  822.             if (rangeStart.compareTo(startEpoch) < 0) {
  823.                 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE, rangeStart, startEpoch, finalEpoch);
  824.             }

  825.             final AbsoluteDate rangeEnd   = extractDate(record, DATE_END_RANGE_OFFSET);
  826.             if (rangeEnd.compareTo(finalEpoch) > 0) {
  827.                 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE, rangeEnd, startEpoch, finalEpoch);
  828.             }

  829.             if (rangeStart.compareTo(end) > 0 || rangeEnd.compareTo(start) < 0) {
  830.                 // we are not interested in this record, don't parse it
  831.                 return rangeEnd;
  832.             }

  833.             // loop over chunks inside the time range
  834.             AbsoluteDate chunkEnd = rangeStart;
  835.             final int nbChunks    = chunks;
  836.             final int nbCoeffs    = coeffs;
  837.             final int first       = firstIndex;
  838.             final double duration = chunksDuration;
  839.             for (int i = 0; i < nbChunks; ++i) {

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

  843.                 // extract Chebyshev coefficients for the selected body
  844.                 // and convert them from kilometers to meters
  845.                 final double[] xCoeffs = new double[nbCoeffs];
  846.                 final double[] yCoeffs = new double[nbCoeffs];
  847.                 final double[] zCoeffs = new double[nbCoeffs];

  848.                 for (int k = 0; k < nbCoeffs; ++k) {
  849.                     // by now, only use the position components
  850.                     // if there are also velocity components contained in the file, ignore them
  851.                     final int index = first + components * i * nbCoeffs + k - 1;
  852.                     xCoeffs[k] = positionUnit * extractDouble(record, 8 * index);
  853.                     yCoeffs[k] = positionUnit * extractDouble(record, 8 * (index +  nbCoeffs));
  854.                     zCoeffs[k] = positionUnit * extractDouble(record, 8 * (index + 2 * nbCoeffs));
  855.                 }

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

  858.             }

  859.             return rangeStart;

  860.         }

  861.     }

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

  864.         /** {@inheritDoc} */
  865.         public PVCoordinates getRawPV(final AbsoluteDate date) {

  866.             // get raw PV from Chebyshev polynomials
  867.             PosVelChebyshev chebyshev;
  868.             try {
  869.                 chebyshev = ephemerides.getNeighbors(date).findFirst().get();
  870.             } catch (TimeStampedCacheException tce) {
  871.                 // we cannot bracket the date, check if the last available chunk covers the specified date
  872.                 chebyshev = ephemerides.getLatest();
  873.                 if (!chebyshev.inRange(date)) {
  874.                     // we were not able to recover from the error, the date is too far
  875.                     throw tce;
  876.                 }
  877.             }

  878.             // evaluate the Chebyshev polynomials
  879.             return chebyshev.getPositionVelocityAcceleration(date);

  880.         }

  881.         /** {@inheritDoc} */
  882.         public <T extends RealFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {

  883.             // get raw PV from Chebyshev polynomials
  884.             PosVelChebyshev chebyshev;
  885.             try {
  886.                 chebyshev = ephemerides.getNeighbors(date.toAbsoluteDate()).findFirst().get();
  887.             } catch (TimeStampedCacheException tce) {
  888.                 // we cannot bracket the date, check if the last available chunk covers the specified date
  889.                 chebyshev = ephemerides.getLatest();
  890.                 if (!chebyshev.inRange(date.toAbsoluteDate())) {
  891.                     // we were not able to recover from the error, the date is too far
  892.                     throw tce;
  893.                 }
  894.             }

  895.             // evaluate the Chebyshev polynomials
  896.             return chebyshev.getPositionVelocityAcceleration(date);

  897.         }

  898.     }

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

  901.         /** {@inheritDoc} */
  902.         public PVCoordinates getRawPV(final AbsoluteDate date) {
  903.             return PVCoordinates.ZERO;
  904.         }

  905.         /** {@inheritDoc} */
  906.         public <T extends RealFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {
  907.             return FieldPVCoordinates.getZero(date.getField());
  908.         }

  909.     }

  910. }