JPLEphemeridesLoader.java

  1. /* Copyright 2002-2018 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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;

  19. import java.io.IOException;
  20. import java.io.InputStream;
  21. import java.io.UnsupportedEncodingException;
  22. import java.text.ParseException;
  23. import java.util.ArrayList;
  24. import java.util.Comparator;
  25. import java.util.HashMap;
  26. import java.util.List;
  27. import java.util.Map;
  28. import java.util.SortedSet;
  29. import java.util.TreeSet;
  30. import java.util.concurrent.atomic.AtomicReference;

  32. import org.hipparchus.RealFieldElement;
  33. import org.hipparchus.util.FastMath;
  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.FramesFactory;
  42. import org.orekit.time.AbsoluteDate;
  43. import org.orekit.time.DateComponents;
  44. import org.orekit.time.FieldAbsoluteDate;
  45. import org.orekit.time.TimeComponents;
  46. import org.orekit.time.TimeScale;
  47. import org.orekit.time.TimeScalesFactory;
  48. import org.orekit.utils.Constants;
  49. import org.orekit.utils.FieldPVCoordinates;
  50. import org.orekit.utils.OrekitConfiguration;
  51. import org.orekit.utils.PVCoordinates;
  52. import org.orekit.utils.GenericTimeStampedCache;
  53. 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 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)$";

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

  81.     /** 50 days in seconds. */
  82.     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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  135.     /** List of supported ephemerides types. */
  136.     public enum EphemerisType {

  138.         /** Constant for solar system barycenter. */
  139.         SOLAR_SYSTEM_BARYCENTER,

  141.         /** Constant for the Sun. */
  142.         SUN,

  144.         /** Constant for Mercury. */
  145.         MERCURY,

  147.         /** Constant for Venus. */
  148.         VENUS,

  150.         /** Constant for the Earth-Moon barycenter. */
  151.         EARTH_MOON,

  153.         /** Constant for the Earth. */
  154.         EARTH,

  156.         /** Constant for the Moon. */
  157.         MOON,

  159.         /** Constant for Mars. */
  160.         MARS,

  162.         /** Constant for Jupiter. */
  163.         JUPITER,

  165.         /** Constant for Saturn. */
  166.         SATURN,

  168.         /** Constant for Uranus. */
  169.         URANUS,

  171.         /** Constant for Neptune. */
  172.         NEPTUNE,

  174.         /** Constant for Pluto. */
  175.         PLUTO

  177.     }

  179.     /** Interface for raw position-velocity retrieval. */
  180.     public interface RawPVProvider {

  182.         /** Get the position-velocity at date.
  183.          * @param date date at which the position-velocity is desired
  184.          * @return position-velocity at the specified date
  185.          * @exception OrekitException if the date is not available to the loader
  186.          */
  187.         PVCoordinates getRawPV(AbsoluteDate date) throws OrekitException;

  189.         /** Get the position-velocity at date.
  190.          * @param date date at which the position-velocity is desired
  191.          * @param <T> type of the field elements
  192.          * @return position-velocity at the specified date
  193.          * @exception OrekitException if the date is not available to the loader
  194.          */
  195.         <T extends RealFieldElement<T>> FieldPVCoordinates<T> getRawPV(FieldAbsoluteDate<T> date) throws OrekitException;

  197.     }

  199.     /** Regular expression for supported files names. */
  200.     private final String supportedNames;

  202.     /** Ephemeris for selected body. */
  203.     private final GenericTimeStampedCache<PosVelChebyshev> ephemerides;

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

  208.     /** Ephemeris type to generate. */
  209.     private final EphemerisType generateType;

  211.     /** Ephemeris type to load. */
  212.     private final EphemerisType loadType;

  214.     /** Current file start epoch. */
  215.     private AbsoluteDate startEpoch;

  217.     /** Current file final epoch. */
  218.     private AbsoluteDate finalEpoch;

  220.     /** Chunks duration (in seconds). */
  221.     private double maxChunksDuration;

  223.     /** Current file chunks duration (in seconds). */
  224.     private double chunksDuration;

  226.     /** Index of the first data for selected body. */
  227.     private int firstIndex;

  229.     /** Number of coefficients for selected body. */
  230.     private int coeffs;

  232.     /** Number of chunks for the selected body. */
  233.     private int chunks;

  235.     /** Number of components contained in the file. */
  236.     private int components;

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

  241.     /** Time scale of the date coordinates. */
  242.     private TimeScale timeScale;

  244.     /** Indicator for binary file endianness. */
  245.     private boolean bigEndian;

  247.     /** Create a loader for JPL ephemerides binary files.
  248.      * @param supportedNames regular expression for supported files names
  249.      * @param generateType ephemeris type to generate
  250.      * @exception OrekitException if the header constants cannot be read
  251.      */
  252.     public JPLEphemeridesLoader(final String supportedNames, final EphemerisType generateType)
  253.         throws OrekitException {

  255.         this.supportedNames = supportedNames;
  256.         constants = new AtomicReference<Map<String, Double>>();

  258.         this.generateType  = generateType;
  259.         if (generateType == EphemerisType.SOLAR_SYSTEM_BARYCENTER) {
  260.             loadType = EphemerisType.EARTH_MOON;
  261.         } else if (generateType == EphemerisType.EARTH_MOON) {
  262.             loadType = EphemerisType.MOON;
  263.         } else {
  264.             loadType = generateType;
  265.         }

  267.         ephemerides = new GenericTimeStampedCache<PosVelChebyshev>(2, OrekitConfiguration.getCacheSlotsNumber(),
  268.                                                                    Double.POSITIVE_INFINITY, FIFTY_DAYS,
  269.                                                                    new EphemerisParser());
  270.         maxChunksDuration = Double.NaN;
  271.         chunksDuration    = Double.NaN;

  273.     }

  275.     /** Load celestial body.
  276.      * @param name name of the celestial body
  277.      * @return loaded celestial body
  278.      * @throws OrekitException if the body cannot be loaded
  279.      */
  280.     public CelestialBody loadCelestialBody(final String name) throws OrekitException {

  282.         final double gm       = getLoadedGravitationalCoefficient(generateType);
  283.         final IAUPole iauPole = PredefinedIAUPoles.getIAUPole(generateType);
  284.         final double scale;
  285.         final Frame definingFrameAlignedWithICRF;
  286.         final RawPVProvider rawPVProvider;
  287.         switch (generateType) {
  288.             case SOLAR_SYSTEM_BARYCENTER : {
  289.                 scale = -1.0;
  290.                 final JPLEphemeridesLoader parentLoader =
  291.                         new JPLEphemeridesLoader(supportedNames, EphemerisType.EARTH_MOON);
  292.                 final CelestialBody parentBody =
  293.                         parentLoader.loadCelestialBody(CelestialBodyFactory.EARTH_MOON);
  294.                 definingFrameAlignedWithICRF = parentBody.getInertiallyOrientedFrame();
  295.                 rawPVProvider = new EphemerisRawPVProvider();
  296.                 break;
  297.             }
  298.             case EARTH_MOON :
  299.                 scale         = 1.0 / (1.0 + getLoadedEarthMoonMassRatio());
  300.                 definingFrameAlignedWithICRF =  FramesFactory.getGCRF();
  301.                 rawPVProvider = new EphemerisRawPVProvider();
  302.                 break;
  303.             case EARTH :
  304.                 scale         = 1.0;
  305.                 definingFrameAlignedWithICRF = FramesFactory.getGCRF();
  306.                 rawPVProvider = new ZeroRawPVProvider();
  307.                 break;
  308.             case MOON :
  309.                 scale         =  1.0;
  310.                 definingFrameAlignedWithICRF =  FramesFactory.getGCRF();
  311.                 rawPVProvider = new EphemerisRawPVProvider();
  312.                 break;
  313.             default : {
  314.                 scale = 1.0;
  315.                 final JPLEphemeridesLoader parentLoader =
  316.                         new JPLEphemeridesLoader(supportedNames, EphemerisType.SOLAR_SYSTEM_BARYCENTER);
  317.                 final CelestialBody parentBody =
  318.                         parentLoader.loadCelestialBody(CelestialBodyFactory.SOLAR_SYSTEM_BARYCENTER);
  319.                 definingFrameAlignedWithICRF = parentBody.getInertiallyOrientedFrame();
  320.                 rawPVProvider = new EphemerisRawPVProvider();
  321.             }
  322.         }

  324.         // build the celestial body
  325.         return new JPLCelestialBody(name, supportedNames, generateType, rawPVProvider,
  326.                                     gm, scale, iauPole, definingFrameAlignedWithICRF);

  328.     }

  330.     /** Get astronomical unit.
  331.      * @return astronomical unit in meters
  332.      * @exception OrekitException if constants cannot be loaded
  333.      */
  334.     public double getLoadedAstronomicalUnit() throws OrekitException {
  335.         return 1000.0 * getLoadedConstant(CONSTANT_AU);
  336.     }

  338.     /** Get Earth/Moon mass ratio.
  339.      * @return Earth/Moon mass ratio
  340.      * @exception OrekitException if constants cannot be loaded
  341.      */
  342.     public double getLoadedEarthMoonMassRatio() throws OrekitException {
  343.         return getLoadedConstant(CONSTANT_EMRAT);
  344.     }

  346.     /** Get the gravitational coefficient of a body.
  347.      * @param body body for which the gravitational coefficient is requested
  348.      * @return gravitational coefficient in m³/s²
  349.      * @exception OrekitException if constants cannot be loaded
  350.      */
  351.     public double getLoadedGravitationalCoefficient(final EphemerisType body)
  352.         throws OrekitException {

  354.         // coefficient in au³/day²
  355.         final double rawGM;
  356.         switch (body) {
  357.             case SOLAR_SYSTEM_BARYCENTER :
  358.                 return getLoadedGravitationalCoefficient(EphemerisType.SUN)        +
  359.                         getLoadedGravitationalCoefficient(EphemerisType.MERCURY)    +
  360.                         getLoadedGravitationalCoefficient(EphemerisType.VENUS)      +
  361.                         getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) +
  362.                         getLoadedGravitationalCoefficient(EphemerisType.MARS)       +
  363.                         getLoadedGravitationalCoefficient(EphemerisType.JUPITER)    +
  364.                         getLoadedGravitationalCoefficient(EphemerisType.SATURN)     +
  365.                         getLoadedGravitationalCoefficient(EphemerisType.URANUS)     +
  366.                         getLoadedGravitationalCoefficient(EphemerisType.NEPTUNE)    +
  367.                         getLoadedGravitationalCoefficient(EphemerisType.PLUTO);
  368.             case SUN :
  369.                 rawGM = getLoadedConstant("GMS", "GM_Sun");
  370.                 break;
  371.             case MERCURY :
  372.                 rawGM = getLoadedConstant("GM1", "GM_Mer");
  373.                 break;
  374.             case VENUS :
  375.                 rawGM = getLoadedConstant("GM2", "GM_Ven");
  376.                 break;
  377.             case EARTH_MOON :
  378.                 rawGM = getLoadedConstant("GMB", "GM_EMB");
  379.                 break;
  380.             case EARTH :
  381.                 return getLoadedEarthMoonMassRatio() *
  382.                         getLoadedGravitationalCoefficient(EphemerisType.MOON);
  383.             case MOON :
  384.                 return getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) /
  385.                         (1.0 + getLoadedEarthMoonMassRatio());
  386.             case MARS :
  387.                 rawGM = getLoadedConstant("GM4", "GM_Mar");
  388.                 break;
  389.             case JUPITER :
  390.                 rawGM = getLoadedConstant("GM5", "GM_Jup");
  391.                 break;
  392.             case SATURN :
  393.                 rawGM = getLoadedConstant("GM6", "GM_Sat");
  394.                 break;
  395.             case URANUS :
  396.                 rawGM = getLoadedConstant("GM7", "GM_Ura");
  397.                 break;
  398.             case NEPTUNE :
  399.                 rawGM = getLoadedConstant("GM8", "GM_Nep");
  400.                 break;
  401.             case PLUTO :
  402.                 rawGM = getLoadedConstant("GM9", "GM_Plu");
  403.                 break;
  404.             default :
  405.                 throw new OrekitInternalError(null);
  406.         }

  408.         final double au    = getLoadedAstronomicalUnit();
  409.         return rawGM * au * au * au / (Constants.JULIAN_DAY * Constants.JULIAN_DAY);

  411.     }

  413.     /** Get a constant defined in the ephemerides headers.
  414.      * <p>Note that since constants are defined in the JPL headers
  415.      * files, they are available as soon as one file is available, even
  416.      * if it doesn't match the desired central date. This is because the
  417.      * header must be parsed before the dates can be checked.</p>
  418.      * <p>
  419.      * There are alternate names for constants since for example JPL names are
  420.      * different from INPOP names (Sun gravity: GMS or GM_Sun, Mars gravity:
  421.      * GM4 or GM_Mar...).
  422.      * </p>
  423.      * @param names alternate names of the constant
  424.      * @return value of the constant of NaN if the constant is not defined
  425.      * @exception OrekitException if constants cannot be loaded
  426.      */
  427.     public double getLoadedConstant(final String... names) throws OrekitException {

  429.         // lazy loading of constants
  430.         Map<String, Double> map = constants.get();
  431.         if (map == null) {
  432.             final ConstantsParser parser = new ConstantsParser();
  433.             if (!DataProvidersManager.getInstance().feed(supportedNames, parser)) {
  434.                 throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
  435.             }
  436.             map = parser.getConstants();
  437.             constants.compareAndSet(null, map);
  438.         }

  440.         for (final String name : names) {
  441.             if (map.containsKey(name)) {
  442.                 return map.get(name).doubleValue();
  443.             }
  444.         }

  446.         return Double.NaN;

  448.     }

  450.     /** Get the maximal chunks duration.
  451.      * @return chunks maximal duration in seconds
  452.      */
  453.     public double getMaxChunksDuration() {
  454.         return maxChunksDuration;
  455.     }

  457.     /** Parse the first header record.
  458.      * @param record first header record
  459.      * @param name name of the file (or zip entry)
  460.      * @exception OrekitException if the header is not a JPL ephemerides binary file header
  461.      */
  462.     private void parseFirstHeaderRecord(final byte[] record, final String name)
  463.         throws OrekitException {

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

  468.         // as default, 3 polynomial coefficients for the Cartesian coordinates
  469.         // (x, y, z) are contained in the file, positions are in kilometers
  470.         // and times are in TDB
  471.         components   = 3;
  472.         positionUnit = 1000.0;
  473.         timeScale    = TimeScalesFactory.getTDB();

  475.         if (deNum == INPOP_DE_NUMBER) {
  476.             // an INPOP file may contain 6 components (including coefficients for the velocity vector)
  477.             final double format = getLoadedConstant("FORMAT");
  478.             if (!Double.isNaN(format) && (int) FastMath.IEEEremainder(format, 10) != 1) {
  479.                 components = 6;
  480.             }

  482.             // INPOP files may have their polynomials expressed in AU
  483.             final double unite = getLoadedConstant("UNITE");
  484.             if (!Double.isNaN(unite) && (int) unite == 0) {
  485.                 positionUnit = getLoadedAstronomicalUnit();
  486.             }

  488.             // INPOP files may have their times expressed in TCB
  489.             final double timesc = getLoadedConstant("TIMESC");
  490.             if (!Double.isNaN(timesc) && (int) timesc == 1) {
  491.                 timeScale = TimeScalesFactory.getTCB();
  492.             }

  494.         }

  496.         // extract covered date range
  497.         startEpoch = extractDate(record, HEADER_START_EPOCH_OFFSET);
  498.         finalEpoch = extractDate(record, HEADER_END_EPOCH_OFFSET);
  499.         boolean ok = finalEpoch.compareTo(startEpoch) > 0;

  501.         // indices of the Chebyshev coefficients for each ephemeris
  502.         for (int i = 0; i < 12; ++i) {
  503.             final int row1 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET     + 12 * i);
  504.             final int row2 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 4 + 12 * i);
  505.             final int row3 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 8 + 12 * i);
  506.             ok = ok && (row1 >= 0) && (row2 >= 0) && (row3 >= 0);
  507.             if (((i ==  0) && (loadType == EphemerisType.MERCURY))    ||
  508.                     ((i ==  1) && (loadType == EphemerisType.VENUS))      ||
  509.                     ((i ==  2) && (loadType == EphemerisType.EARTH_MOON)) ||
  510.                     ((i ==  3) && (loadType == EphemerisType.MARS))       ||
  511.                     ((i ==  4) && (loadType == EphemerisType.JUPITER))    ||
  512.                     ((i ==  5) && (loadType == EphemerisType.SATURN))     ||
  513.                     ((i ==  6) && (loadType == EphemerisType.URANUS))     ||
  514.                     ((i ==  7) && (loadType == EphemerisType.NEPTUNE))    ||
  515.                     ((i ==  8) && (loadType == EphemerisType.PLUTO))      ||
  516.                     ((i ==  9) && (loadType == EphemerisType.MOON))       ||
  517.                     ((i == 10) && (loadType == EphemerisType.SUN))) {
  518.                 firstIndex = row1;
  519.                 coeffs     = row2;
  520.                 chunks     = row3;
  521.             }
  522.         }

  524.         // compute chunks duration
  525.         final double timeSpan = extractDouble(record, HEADER_CHUNK_DURATION_OFFSET);
  526.         ok = ok && (timeSpan > 0) && (timeSpan < 100);
  527.         chunksDuration = Constants.JULIAN_DAY * (timeSpan / chunks);
  528.         if (Double.isNaN(maxChunksDuration)) {
  529.             maxChunksDuration = chunksDuration;
  530.         } else {
  531.             maxChunksDuration = FastMath.max(maxChunksDuration, chunksDuration);
  532.         }

  534.         // sanity checks
  535.         if (!ok) {
  536.             throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  537.         }

  539.     }

  541.     /** Read first header record.
  542.      * @param input input stream
  543.      * @param name name of the file (or zip entry)
  544.      * @return record record where to put bytes
  545.      * @exception OrekitException if the stream does not contain a JPL ephemeris
  546.      * @exception IOException if a read error occurs
  547.      */
  548.     private byte[] readFirstRecord(final InputStream input, final String name)
  549.         throws OrekitException, IOException {

  551.         // read first part of record, up to the ephemeris type
  552.         final byte[] firstPart = new byte[HEADER_RECORD_SIZE_OFFSET + 4];
  553.         if (!readInRecord(input, firstPart, 0)) {
  554.             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  555.         }

  557.         // detect the endian format
  558.         detectEndianess(firstPart);

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

  563.         // the record size for this file
  564.         int recordSize = 0;

  566.         if (deNum == INPOP_DE_NUMBER) {
  567.             // INPOP files have an extended DE format, which includes also the record size
  568.             recordSize = extractInt(firstPart, HEADER_RECORD_SIZE_OFFSET) << 3;
  569.         } else {
  570.             // compute the record size for original JPL files
  571.             recordSize = computeRecordSize(firstPart, name);
  572.         }

  574.         if (recordSize <= 0) {
  575.             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  576.         }

  578.         // build a record with the proper size and finish read of the first complete record
  579.         final int start = firstPart.length;
  580.         final byte[] record = new byte[recordSize];
  581.         System.arraycopy(firstPart, 0, record, 0, firstPart.length);
  582.         if (!readInRecord(input, record, start)) {
  583.             throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  584.         }

  586.         return record;

  588.     }

  590.     /** Parse constants from first two header records.
  591.      * @param first first header record
  592.      * @param second second header record
  593.      * @param name name of the file (or zip entry)
  594.      * @return map of parsed constants
  595.      */
  596.     private Map<String, Double> parseConstants(final byte[] first, final byte[] second, final String name) {

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

  600.         for (int i = 0; i < CONSTANTS_MAX_NUMBER; ++i) {
  601.             // Note: for extracting the strings from the binary file, it makes no difference
  602.             //       if the file is stored in big-endian or little-endian notation
  603.             final String constantName = extractString(first, HEADER_CONSTANTS_NAMES_OFFSET + i * 6, 6);
  604.             if (constantName.length() == 0) {
  605.                 // no more constants to read
  606.                 break;
  607.             }
  608.             final double constantValue = extractDouble(second, HEADER_CONSTANTS_VALUES_OFFSET + 8 * i);
  609.             map.put(constantName, constantValue);
  610.         }

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

  618.         if (!map.containsKey(CONSTANT_EMRAT)) {
  619.             map.put(CONSTANT_EMRAT, extractDouble(first, HEADER_EM_RATIO_OFFSET));
  620.         }

  622.         return map;

  624.     }

  626.     /** Read bytes into the current record array.
  627.      * @param input input stream
  628.      * @param record record where to put bytes
  629.      * @param start start index where to put bytes
  630.      * @return true if record has been filled up
  631.      * @exception IOException if a read error occurs
  632.      */
  633.     private boolean readInRecord(final InputStream input, final byte[] record, final int start)
  634.         throws IOException {
  635.         int index = start;
  636.         while (index != record.length) {
  637.             final int n = input.read(record, index, record.length - index);
  638.             if (n < 0) {
  639.                 return false;
  640.             }
  641.             index += n;
  642.         }
  643.         return true;
  644.     }

  646.     /** Detect whether the JPL ephemerides file is stored in big-endian or
  647.      * little-endian notation.
  648.      * @param record the array containing the binary JPL header
  649.      */
  650.     private void detectEndianess(final byte[] record) {

  652.         // default to big-endian
  653.         bigEndian = true;

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

  659.         // simple heuristic: if the read value is larger than half the range of an integer
  660.         //                   assume the file is in little-endian format
  661.         if (deNum > (1 << 15)) {
  662.             bigEndian = false;
  663.         }

  665.     }

  667.     /** Calculate the record size of a JPL ephemerides file.
  668.      * @param record the byte array containing the header record
  669.      * @param name the name of the data file
  670.      * @return the record size for this file
  671.      * @throws OrekitException if the file contains unexpected data
  672.      */
  673.     private int computeRecordSize(final byte[] record, final String name)
  674.         throws OrekitException {

  676.         int recordSize = 0;
  677.         boolean ok = true;
  678.         // JPL files always have 3 position components
  679.         final int nComp = 3;

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

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

  691.             // sanity checks
  692.             ok = ok && (coeffPtr1 >= 0 || coeffPtr2 >= 0);

  694.             recordSize += coeffPtr1 * coeffPtr2 * nCompCur;
  695.         }

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

  702.         // sanity checks
  703.         ok = ok && (libratPtr1 >= 0 || libratPtr2 >= 0);

  705.         recordSize += libratPtr1 * libratPtr2 * nComp + 2;
  706.         recordSize <<= 3;

  708.         if (!ok || recordSize <= 0) {
  709.             throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  710.         }

  712.         return recordSize;

  714.     }

  716.     /** Extract a date from a record.
  717.      * @param record record to parse
  718.      * @param offset offset of the double within the record
  719.      * @return extracted date
  720.      */
  721.     private AbsoluteDate extractDate(final byte[] record, final int offset) {

  723.         final double t = extractDouble(record, offset);
  724.         int    jDay    = (int) FastMath.floor(t);
  725.         double seconds = (t + 0.5 - jDay) * Constants.JULIAN_DAY;
  726.         if (seconds >= Constants.JULIAN_DAY) {
  727.             ++jDay;
  728.             seconds -= Constants.JULIAN_DAY;
  729.         }
  730.         return new AbsoluteDate(new DateComponents(DateComponents.JULIAN_EPOCH, jDay),
  731.                                 new TimeComponents(seconds), timeScale);
  732.     }

  734.     /** Extract a double from a record.
  735.      * <p>Double numbers are stored according to IEEE 754 standard, with
  736.      * most significant byte first.</p>
  737.      * @param record record to parse
  738.      * @param offset offset of the double within the record
  739.      * @return extracted double
  740.      */
  741.     private double extractDouble(final byte[] record, final int offset) {
  742.         final long l8 = ((long) record[offset + 0]) & 0xffl;
  743.         final long l7 = ((long) record[offset + 1]) & 0xffl;
  744.         final long l6 = ((long) record[offset + 2]) & 0xffl;
  745.         final long l5 = ((long) record[offset + 3]) & 0xffl;
  746.         final long l4 = ((long) record[offset + 4]) & 0xffl;
  747.         final long l3 = ((long) record[offset + 5]) & 0xffl;
  748.         final long l2 = ((long) record[offset + 6]) & 0xffl;
  749.         final long l1 = ((long) record[offset + 7]) & 0xffl;
  750.         final long l;
  751.         if (bigEndian) {
  752.             l = (l8 << 56) | (l7 << 48) | (l6 << 40) | (l5 << 32) |
  753.                 (l4 << 24) | (l3 << 16) | (l2 <<  8) | l1;
  754.         } else {
  755.             l = (l1 << 56) | (l2 << 48) | (l3 << 40) | (l4 << 32) |
  756.                 (l5 << 24) | (l6 << 16) | (l7 <<  8) | l8;
  757.         }
  758.         return Double.longBitsToDouble(l);
  759.     }

  761.     /** Extract an int from a record.
  762.      * @param record record to parse
  763.      * @param offset offset of the double within the record
  764.      * @return extracted int
  765.      */
  766.     private int extractInt(final byte[] record, final int offset) {
  767.         final int l4 = ((int) record[offset + 0]) & 0xff;
  768.         final int l3 = ((int) record[offset + 1]) & 0xff;
  769.         final int l2 = ((int) record[offset + 2]) & 0xff;
  770.         final int l1 = ((int) record[offset + 3]) & 0xff;

  772.         if (bigEndian) {
  773.             return (l4 << 24) | (l3 << 16) | (l2 <<  8) | l1;
  774.         } else {
  775.             return (l1 << 24) | (l2 << 16) | (l3 <<  8) | l4;
  776.         }
  777.     }

  779.     /** Extract a String from a record.
  780.      * @param record record to parse
  781.      * @param offset offset of the string within the record
  782.      * @param length maximal length of the string
  783.      * @return extracted string, with whitespace characters stripped
  784.      */
  785.     private String extractString(final byte[] record, final int offset, final int length) {
  786.         try {
  787.             return new String(record, offset, length, "US-ASCII").trim();
  788.         } catch (UnsupportedEncodingException uee) {
  789.             throw new OrekitInternalError(uee);
  790.         }
  791.     }

  793.     /** Local parser for header constants. */
  794.     private class ConstantsParser implements DataLoader {

  796.         /** Local constants map. */
  797.         private Map<String, Double> localConstants;

  799.        /** Get the local constants map.
  800.          * @return local constants map
  801.          */
  802.         public Map<String, Double> getConstants() {
  803.             return localConstants;
  804.         }

  806.         /** {@inheritDoc} */
  807.         public boolean stillAcceptsData() {
  808.             return localConstants == null;
  809.         }

  811.         /** {@inheritDoc} */
  812.         public void loadData(final InputStream input, final String name)
  813.             throws IOException, ParseException, OrekitException {

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

  818.             // the second record contains the values of the constants used for least-square filtering
  819.             final byte[] second = new byte[first.length];
  820.             if (!readInRecord(input, second, 0)) {
  821.                 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  822.             }

  824.             localConstants = parseConstants(first, second, name);

  826.         }

  828.     }

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

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

  836.         /** Start of range we are interested in. */
  837.         private AbsoluteDate start;

  839.         /** End of range we are interested in. */
  840.         private AbsoluteDate end;

  842.         /** Simple constructor.
  843.          */
  844.         EphemerisParser() {
  845.             entries = new TreeSet<PosVelChebyshev>(new Comparator<PosVelChebyshev>() {
  846.                 public int compare(final PosVelChebyshev o1, final PosVelChebyshev o2) {
  847.                     return o1.getDate().compareTo(o2.getDate());
  848.                 }
  849.             });
  850.         }

  852.         /** {@inheritDoc} */
  853.         public List<PosVelChebyshev> generate(final AbsoluteDate existingDate, final AbsoluteDate date)
  854.             throws TimeStampedCacheException {
  855.             try {

  857.                 // prepare reading
  858.                 entries.clear();
  859.                 if (existingDate == null) {
  860.                     // we want ephemeris data for the first time, set up an arbitrary first range
  861.                     start = date.shiftedBy(-FIFTY_DAYS);
  862.                     end   = date.shiftedBy(+FIFTY_DAYS);
  863.                 } else if (existingDate.compareTo(date) <= 0) {
  864.                     // we want to extend an existing range towards future dates
  865.                     start = existingDate;
  866.                     end   = date;
  867.                 } else {
  868.                     // we want to extend an existing range towards past dates
  869.                     start = date;
  870.                     end   = existingDate;
  871.                 }

  873.                 // get new entries in the specified data range
  874.                 if (!DataProvidersManager.getInstance().feed(supportedNames, this)) {
  875.                     throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
  876.                 }

  878.                 return new ArrayList<PosVelChebyshev>(entries);

  880.             } catch (OrekitException oe) {
  881.                 throw new TimeStampedCacheException(oe);
  882.             }
  883.         }

  885.         /** {@inheritDoc} */
  886.         public boolean stillAcceptsData() {

  888.             // special case for Earth: we do not really load any ephemeris data
  889.             if (generateType == EphemerisType.EARTH) {
  890.                 return false;
  891.             }

  893.             // we have to look for data in all available ephemerides files as there may be
  894.             // data overlaps that result in incomplete data
  895.             if (entries.isEmpty()) {
  896.                 return true;
  897.             } else {
  898.                 // if the requested range is already filled, we do not need to look further
  899.                 return !(entries.first().getDate().compareTo(start) < 0 &&
  900.                          entries.last().getDate().compareTo(end)    > 0);
  901.             }

  903.         }

  905.         /** {@inheritDoc} */
  906.         public void loadData(final InputStream input, final String name)
  907.             throws OrekitException, IOException {

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

  912.             // the second record contains the values of the constants used for least-square filtering
  913.             final byte[] second = new byte[first.length];
  914.             if (!readInRecord(input, second, 0)) {
  915.                 throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
  916.             }

  918.             if (constants.get() == null) {
  919.                 constants.compareAndSet(null, parseConstants(first, second, name));
  920.             }

  922.             // check astronomical unit consistency
  923.             final double au = 1000 * extractDouble(first, HEADER_ASTRONOMICAL_UNIT_OFFSET);
  924.             if ((au < 1.4e11) || (au > 1.6e11)) {
  925.                 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  926.             }
  927.             if (FastMath.abs(getLoadedAstronomicalUnit() - au) >= 10.0) {
  928.                 throw new OrekitException(OrekitMessages.INCONSISTENT_ASTRONOMICAL_UNIT_IN_FILES,
  929.                                           getLoadedAstronomicalUnit(), au);
  930.             }

  932.             // check Earth-Moon mass ratio consistency
  933.             final double emRat = extractDouble(first, HEADER_EM_RATIO_OFFSET);
  934.             if ((emRat < 80) || (emRat > 82)) {
  935.                 throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
  936.             }
  937.             if (FastMath.abs(getLoadedEarthMoonMassRatio() - emRat) >= 1.0e-5) {
  938.                 throw new OrekitException(OrekitMessages.INCONSISTENT_EARTH_MOON_RATIO_IN_FILES,
  939.                                           getLoadedEarthMoonMassRatio(), emRat);
  940.             }

  942.             // parse first header record
  943.             parseFirstHeaderRecord(first, name);

  945.             if (startEpoch.compareTo(end) < 0 && finalEpoch.compareTo(start) > 0) {
  946.                 // this file contains data in the range we are looking for, read it
  947.                 final byte[] record = new byte[first.length];
  948.                 while (readInRecord(input, record, 0)) {
  949.                     final AbsoluteDate rangeStart = parseDataRecord(record);
  950.                     if (rangeStart.compareTo(end) > 0) {
  951.                         // we have already exceeded the range we were interested in,
  952.                         // we interrupt parsing here
  953.                         return;
  954.                     }
  955.                 }
  956.             }

  958.         }

  960.         /** Parse regular ephemeris record.
  961.          * @param record record to parse
  962.          * @return date of the last parsed chunk
  963.          * @exception OrekitException if the header is not a JPL ephemerides binary file header
  964.          */
  965.         private AbsoluteDate parseDataRecord(final byte[] record) throws OrekitException {

  967.             // extract time range covered by the record
  968.             final AbsoluteDate rangeStart = extractDate(record, DATA_START_RANGE_OFFSET);
  969.             if (rangeStart.compareTo(startEpoch) < 0) {
  970.                 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE, rangeStart, startEpoch, finalEpoch);
  971.             }

  973.             final AbsoluteDate rangeEnd   = extractDate(record, DATE_END_RANGE_OFFSET);
  974.             if (rangeEnd.compareTo(finalEpoch) > 0) {
  975.                 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE, rangeEnd, startEpoch, finalEpoch);
  976.             }

  978.             if (rangeStart.compareTo(end) > 0 || rangeEnd.compareTo(start) < 0) {
  979.                 // we are not interested in this record, don't parse it
  980.                 return rangeEnd;
  981.             }

  983.             // loop over chunks inside the time range
  984.             AbsoluteDate chunkEnd = rangeStart;
  985.             final int nbChunks    = chunks;
  986.             final int nbCoeffs    = coeffs;
  987.             final int first       = firstIndex;
  988.             final double duration = chunksDuration;
  989.             for (int i = 0; i < nbChunks; ++i) {

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

  995.                 // extract Chebyshev coefficients for the selected body
  996.                 // and convert them from kilometers to meters
  997.                 final double[] xCoeffs = new double[nbCoeffs];
  998.                 final double[] yCoeffs = new double[nbCoeffs];
  999.                 final double[] zCoeffs = new double[nbCoeffs];

  1001.                 for (int k = 0; k < nbCoeffs; ++k) {
  1002.                     // by now, only use the position components
  1003.                     // if there are also velocity components contained in the file, ignore them
  1004.                     final int index = first + components * i * nbCoeffs + k - 1;
  1005.                     xCoeffs[k] = positionUnit * extractDouble(record, 8 * index);
  1006.                     yCoeffs[k] = positionUnit * extractDouble(record, 8 * (index +  nbCoeffs));
  1007.                     zCoeffs[k] = positionUnit * extractDouble(record, 8 * (index + 2 * nbCoeffs));
  1008.                 }

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

  1013.             }

  1015.             return rangeStart;

  1017.         }

  1019.     }

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

  1024.         /** {@inheritDoc} */
  1025.         public PVCoordinates getRawPV(final AbsoluteDate date) throws TimeStampedCacheException {

  1027.             // get raw PV from Chebyshev polynomials
  1028.             PosVelChebyshev chebyshev;
  1029.             try {
  1030.                 chebyshev = ephemerides.getNeighbors(date).findFirst().get();
  1031.             } catch (TimeStampedCacheException tce) {
  1032.                 // we cannot bracket the date, check if the last available chunk covers the specified date
  1033.                 chebyshev = ephemerides.getLatest();
  1034.                 if (!chebyshev.inRange(date)) {
  1035.                     // we were not able to recover from the error, the date is too far
  1036.                     throw tce;
  1037.                 }
  1038.             }

  1040.             // evaluate the Chebyshev polynomials
  1041.             return chebyshev.getPositionVelocityAcceleration(date);

  1043.         }

  1045.         /** {@inheritDoc} */
  1046.         public <T extends RealFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date)
  1047.             throws TimeStampedCacheException {

  1049.             // get raw PV from Chebyshev polynomials
  1050.             PosVelChebyshev chebyshev;
  1051.             try {
  1052.                 chebyshev = ephemerides.getNeighbors(date.toAbsoluteDate()).findFirst().get();
  1053.             } catch (TimeStampedCacheException tce) {
  1054.                 // we cannot bracket the date, check if the last available chunk covers the specified date
  1055.                 chebyshev = ephemerides.getLatest();
  1056.                 if (!chebyshev.inRange(date.toAbsoluteDate())) {
  1057.                     // we were not able to recover from the error, the date is too far
  1058.                     throw tce;
  1059.                 }
  1060.             }

  1062.             // evaluate the Chebyshev polynomials
  1063.             return chebyshev.getPositionVelocityAcceleration(date);

  1065.         }

  1067.     }

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

  1072.         /** {@inheritDoc} */
  1073.         public PVCoordinates getRawPV(final AbsoluteDate date) {
  1074.             return PVCoordinates.ZERO;
  1075.         }

  1077.         /** {@inheritDoc} */
  1078.         public <T extends RealFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {
  1079.             return FieldPVCoordinates.getZero(date.getField());
  1080.         }

  1082.     }

  1084. }
Created with JaCoCo 0.7.9.201702052155