JPLEphemeridesLoader.java
- /* Copyright 2002-2024 CS GROUP
- * Licensed to CS GROUP (CS) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * CS licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
- package org.orekit.bodies;
- import java.io.IOException;
- import java.io.InputStream;
- import java.nio.charset.StandardCharsets;
- import java.text.ParseException;
- import java.util.ArrayList;
- import java.util.HashMap;
- import java.util.List;
- import java.util.Map;
- import java.util.SortedSet;
- import java.util.TreeSet;
- import java.util.concurrent.atomic.AtomicReference;
- import org.hipparchus.CalculusFieldElement;
- import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
- import org.hipparchus.geometry.euclidean.threed.Vector3D;
- import org.hipparchus.util.FastMath;
- import org.orekit.annotation.DefaultDataContext;
- import org.orekit.data.AbstractSelfFeedingLoader;
- import org.orekit.data.DataContext;
- import org.orekit.data.DataLoader;
- import org.orekit.data.DataProvidersManager;
- import org.orekit.errors.OrekitException;
- import org.orekit.errors.OrekitInternalError;
- import org.orekit.errors.OrekitMessages;
- import org.orekit.errors.TimeStampedCacheException;
- import org.orekit.frames.Frame;
- import org.orekit.frames.Predefined;
- import org.orekit.time.AbsoluteDate;
- import org.orekit.time.ChronologicalComparator;
- import org.orekit.time.DateComponents;
- import org.orekit.time.FieldAbsoluteDate;
- import org.orekit.time.TimeComponents;
- import org.orekit.time.TimeScale;
- import org.orekit.time.TimeScales;
- import org.orekit.utils.Constants;
- import org.orekit.utils.FieldPVCoordinates;
- import org.orekit.utils.OrekitConfiguration;
- import org.orekit.utils.PVCoordinates;
- import org.orekit.utils.GenericTimeStampedCache;
- import org.orekit.utils.TimeStampedGenerator;
- import org.orekit.utils.units.UnitsConverter;
- /** Loader for JPL ephemerides binary files (DE 4xx) and similar formats (INPOP 06/08/10).
- * <p>JPL ephemerides binary files contain ephemerides for all solar system planets.</p>
- * <p>The JPL ephemerides binary files are recognized thanks to their base names,
- * which must match the pattern <code>[lu]nx[mp]####.ddd</code> (or
- * <code>[lu]nx[mp]####.ddd.gz</code> for gzip-compressed files) where # stands for a
- * digit character and where ddd is an ephemeris type (typically 405 or 406).</p>
- * <p>The loader supports files encoded in big-endian as well as in little-endian notation.
- * Usually, big-endian files are named <code>unx[mp]####.ddd</code>, while little-endian files
- * are named <code>lnx[mp]####.ddd</code>.</p>
- * <p>The IMCCE ephemerides binary files are recognized thanks to their base names,
- * which must match the pattern <code>inpop*.dat</code> (or
- * <code>inpop*.dat.gz</code> for gzip-compressed files) where * stands for any string.</p>
- * <p>The loader supports files encoded in big-endian as well as in little-endian notation.
- * Usually, big-endian files contain <code>bigendian</code> in their names, while little-endian files
- * contain <code>littleendian</code> in their names.</p>
- * <p>The loader supports files in TDB or TCB time scales.</p>
- * @author Luc Maisonobe
- */
- public class JPLEphemeridesLoader extends AbstractSelfFeedingLoader
- implements CelestialBodyLoader {
- /** Default supported files name pattern for JPL DE files. */
- public static final String DEFAULT_DE_SUPPORTED_NAMES = "^[lu]nx([mp](\\d\\d\\d\\d))+\\.(?:4\\d\\d)$";
- /** Default supported files name pattern for IMCCE INPOP files. */
- public static final String DEFAULT_INPOP_SUPPORTED_NAMES = "^inpop.*\\.dat$";
- /** 50 days in seconds. */
- private static final double FIFTY_DAYS = 50 * Constants.JULIAN_DAY;
- /** DE number used by INPOP files. */
- private static final int INPOP_DE_NUMBER = 100;
- /** Maximal number of constants in headers. */
- private static final int CONSTANTS_MAX_NUMBER = 400;
- /** Offset of the ephemeris type in first header record. */
- private static final int HEADER_EPHEMERIS_TYPE_OFFSET = 2840;
- /** Offset of the record size (for INPOP files) in first header record. */
- private static final int HEADER_RECORD_SIZE_OFFSET = 2856;
- /** Offset of the start epoch in first header record. */
- private static final int HEADER_START_EPOCH_OFFSET = 2652;
- /** Offset of the end epoch in first header record. */
- private static final int HEADER_END_EPOCH_OFFSET = 2660;
- /** Offset of the astronomical unit in first header record. */
- private static final int HEADER_ASTRONOMICAL_UNIT_OFFSET = 2680;
- /** Offset of the Earth-Moon mass ratio in first header record. */
- private static final int HEADER_EM_RATIO_OFFSET = 2688;
- /** Offset of Chebishev coefficients indices in first header record. */
- private static final int HEADER_CHEBISHEV_INDICES_OFFSET = 2696;
- /** Offset of libration coefficients indices in first header record. */
- private static final int HEADER_LIBRATION_INDICES_OFFSET = 2844;
- /** Offset of chunks duration in first header record. */
- private static final int HEADER_CHUNK_DURATION_OFFSET = 2668;
- /** Offset of the constants names in first header record. */
- private static final int HEADER_CONSTANTS_NAMES_OFFSET = 252;
- /** Offset of the constants values in second header record. */
- private static final int HEADER_CONSTANTS_VALUES_OFFSET = 0;
- /** Offset of the range start in the data records. */
- private static final int DATA_START_RANGE_OFFSET = 0;
- /** Offset of the range end in the data records. */
- private static final int DATE_END_RANGE_OFFSET = 8;
- /** The constant name for the astronomical unit. */
- private static final String CONSTANT_AU = "AU";
- /** The constant name for the earth-moon mass ratio. */
- private static final String CONSTANT_EMRAT = "EMRAT";
- /** List of supported ephemerides types. */
- public enum EphemerisType {
- /** Constant for solar system barycenter. */
- SOLAR_SYSTEM_BARYCENTER,
- /** Constant for the Sun. */
- SUN,
- /** Constant for Mercury. */
- MERCURY,
- /** Constant for Venus. */
- VENUS,
- /** Constant for the Earth-Moon barycenter. */
- EARTH_MOON,
- /** Constant for the Earth. */
- EARTH,
- /** Constant for the Moon. */
- MOON,
- /** Constant for Mars. */
- MARS,
- /** Constant for Jupiter. */
- JUPITER,
- /** Constant for Saturn. */
- SATURN,
- /** Constant for Uranus. */
- URANUS,
- /** Constant for Neptune. */
- NEPTUNE,
- /** Constant for Pluto. */
- PLUTO
- }
- /** Interface for raw position-velocity retrieval. */
- public interface RawPVProvider {
- /** Get the position-velocity at date.
- * @param date date at which the position-velocity is desired
- * @return position-velocity at the specified date
- */
- PVCoordinates getRawPV(AbsoluteDate date);
- /** Get the position at date.
- * @param date date at which the position is desired
- * @return position at the specified date
- */
- default Vector3D getRawPosition(final AbsoluteDate date) {
- return getRawPV(date).getPosition();
- }
- /** Get the position-velocity at date.
- * @param date date at which the position-velocity is desired
- * @param <T> type of the field elements
- * @return position-velocity at the specified date
- */
- <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(FieldAbsoluteDate<T> date);
- /** Get the position at date.
- * @param date date at which the position is desired
- * @param <T> type of the field elements
- * @return position at the specified date
- */
- default <T extends CalculusFieldElement<T>> FieldVector3D<T> getRawPosition(final FieldAbsoluteDate<T> date) {
- return getRawPV(date).getPosition();
- }
- }
- /** Ephemeris for selected body. */
- private final GenericTimeStampedCache<PosVelChebyshev> ephemerides;
- /** Constants defined in the file. */
- private final AtomicReference<Map<String, Double>> constants;
- /** Ephemeris type to generate. */
- private final EphemerisType generateType;
- /** Ephemeris type to load. */
- private final EphemerisType loadType;
- /** Time scales to use when loading data. */
- private final TimeScales timeScales;
- /** The GCRF implementation. */
- private final Frame gcrf;
- /** Current file start epoch. */
- private AbsoluteDate startEpoch;
- /** Current file final epoch. */
- private AbsoluteDate finalEpoch;
- /** Chunks duration (in seconds). */
- private double maxChunksDuration;
- /** Current file chunks duration (in seconds). */
- private double chunksDuration;
- /** Index of the first data for selected body. */
- private int firstIndex;
- /** Number of coefficients for selected body. */
- private int coeffs;
- /** Number of chunks for the selected body. */
- private int chunks;
- /** Number of components contained in the file. */
- private int components;
- /** Unit of the position coordinates (as a multiple of meters). */
- private double positionUnit;
- /** Time scale of the date coordinates. */
- private TimeScale timeScale;
- /** Indicator for binary file endianness. */
- private boolean bigEndian;
- /** Create a loader for JPL ephemerides binary files. This constructor uses the {@link
- * DataContext#getDefault() default data context}.
- *
- * @param supportedNames regular expression for supported files names
- * @param generateType ephemeris type to generate
- * @see #JPLEphemeridesLoader(String, EphemerisType, DataProvidersManager, TimeScales,
- * Frame)
- */
- @DefaultDataContext
- public JPLEphemeridesLoader(final String supportedNames, final EphemerisType generateType) {
- this(supportedNames, generateType,
- DataContext.getDefault().getDataProvidersManager(),
- DataContext.getDefault().getTimeScales(),
- DataContext.getDefault().getFrames().getGCRF());
- }
- /** Create a loader for JPL ephemerides binary files.
- * @param supportedNames regular expression for supported files names
- * @param generateType ephemeris type to generate
- * @param dataProvidersManager provides access to the ephemeris files.
- * @param timeScales used to access the TCB and TDB time scales while loading data.
- * @param gcrf Earth centered frame aligned with ICRF.
- * @since 10.1
- */
- public JPLEphemeridesLoader(final String supportedNames,
- final EphemerisType generateType,
- final DataProvidersManager dataProvidersManager,
- final TimeScales timeScales,
- final Frame gcrf) {
- super(supportedNames, dataProvidersManager);
- this.timeScales = timeScales;
- this.gcrf = gcrf;
- constants = new AtomicReference<>();
- this.generateType = generateType;
- if (generateType == EphemerisType.SOLAR_SYSTEM_BARYCENTER) {
- loadType = EphemerisType.EARTH_MOON;
- } else if (generateType == EphemerisType.EARTH_MOON) {
- loadType = EphemerisType.MOON;
- } else {
- loadType = generateType;
- }
- ephemerides = new GenericTimeStampedCache<>(
- 2, OrekitConfiguration.getCacheSlotsNumber(),
- Double.POSITIVE_INFINITY, FIFTY_DAYS,
- new EphemerisParser());
- maxChunksDuration = Double.NaN;
- chunksDuration = Double.NaN;
- }
- /** Load celestial body.
- * @param name name of the celestial body
- * @return loaded celestial body
- */
- public CelestialBody loadCelestialBody(final String name) {
- final double gm = getLoadedGravitationalCoefficient(generateType);
- final IAUPole iauPole = PredefinedIAUPoles
- .getIAUPole(generateType, timeScales);
- final double scale;
- final Frame definingFrameAlignedWithICRF;
- final RawPVProvider rawPVProvider;
- String inertialFrameName = null;
- String bodyOrientedFrameName = null;
- switch (generateType) {
- case SOLAR_SYSTEM_BARYCENTER : {
- scale = -1.0;
- final JPLEphemeridesLoader parentLoader = new JPLEphemeridesLoader(
- getSupportedNames(),
- EphemerisType.EARTH_MOON,
- getDataProvidersManager(),
- timeScales,
- gcrf);
- final CelestialBody parentBody =
- parentLoader.loadCelestialBody(CelestialBodyFactory.EARTH_MOON);
- definingFrameAlignedWithICRF = parentBody.getInertiallyOrientedFrame();
- rawPVProvider = new EphemerisRawPVProvider();
- inertialFrameName = Predefined.ICRF.getName();
- bodyOrientedFrameName = null;
- break;
- }
- case EARTH_MOON :
- scale = 1.0 / (1.0 + getLoadedEarthMoonMassRatio());
- definingFrameAlignedWithICRF = gcrf;
- rawPVProvider = new EphemerisRawPVProvider();
- break;
- case EARTH :
- scale = 1.0;
- definingFrameAlignedWithICRF = gcrf;
- rawPVProvider = new ZeroRawPVProvider();
- break;
- case MOON :
- scale = 1.0;
- definingFrameAlignedWithICRF = gcrf;
- rawPVProvider = new EphemerisRawPVProvider();
- break;
- default : {
- scale = 1.0;
- final JPLEphemeridesLoader parentLoader = new JPLEphemeridesLoader(
- getSupportedNames(),
- EphemerisType.SOLAR_SYSTEM_BARYCENTER,
- getDataProvidersManager(),
- timeScales,
- gcrf);
- final CelestialBody parentBody =
- parentLoader.loadCelestialBody(CelestialBodyFactory.SOLAR_SYSTEM_BARYCENTER);
- definingFrameAlignedWithICRF = parentBody.getInertiallyOrientedFrame();
- rawPVProvider = new EphemerisRawPVProvider();
- }
- }
- // build the celestial body
- return new JPLCelestialBody(name, getSupportedNames(), generateType, rawPVProvider,
- gm, scale, iauPole, definingFrameAlignedWithICRF,
- inertialFrameName, bodyOrientedFrameName);
- }
- /** Get astronomical unit.
- * @return astronomical unit in meters
- */
- public double getLoadedAstronomicalUnit() {
- return UnitsConverter.KILOMETRES_TO_METRES.convert(getLoadedConstant(CONSTANT_AU));
- }
- /** Get Earth/Moon mass ratio.
- * @return Earth/Moon mass ratio
- */
- public double getLoadedEarthMoonMassRatio() {
- return getLoadedConstant(CONSTANT_EMRAT);
- }
- /** Get the gravitational coefficient of a body.
- * @param body body for which the gravitational coefficient is requested
- * @return gravitational coefficient in m³/s²
- */
- public double getLoadedGravitationalCoefficient(final EphemerisType body) {
- // coefficient in au³/day²
- final double rawGM;
- switch (body) {
- case SOLAR_SYSTEM_BARYCENTER :
- return getLoadedGravitationalCoefficient(EphemerisType.SUN) +
- getLoadedGravitationalCoefficient(EphemerisType.MERCURY) +
- getLoadedGravitationalCoefficient(EphemerisType.VENUS) +
- getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) +
- getLoadedGravitationalCoefficient(EphemerisType.MARS) +
- getLoadedGravitationalCoefficient(EphemerisType.JUPITER) +
- getLoadedGravitationalCoefficient(EphemerisType.SATURN) +
- getLoadedGravitationalCoefficient(EphemerisType.URANUS) +
- getLoadedGravitationalCoefficient(EphemerisType.NEPTUNE) +
- getLoadedGravitationalCoefficient(EphemerisType.PLUTO);
- case SUN :
- rawGM = getLoadedConstant("GMS", "GM_Sun");
- break;
- case MERCURY :
- rawGM = getLoadedConstant("GM1", "GM_Mer");
- break;
- case VENUS :
- rawGM = getLoadedConstant("GM2", "GM_Ven");
- break;
- case EARTH_MOON :
- rawGM = getLoadedConstant("GMB", "GM_EMB");
- break;
- case EARTH :
- return getLoadedEarthMoonMassRatio() *
- getLoadedGravitationalCoefficient(EphemerisType.MOON);
- case MOON :
- return getLoadedGravitationalCoefficient(EphemerisType.EARTH_MOON) /
- (1.0 + getLoadedEarthMoonMassRatio());
- case MARS :
- rawGM = getLoadedConstant("GM4", "GM_Mar");
- break;
- case JUPITER :
- rawGM = getLoadedConstant("GM5", "GM_Jup");
- break;
- case SATURN :
- rawGM = getLoadedConstant("GM6", "GM_Sat");
- break;
- case URANUS :
- rawGM = getLoadedConstant("GM7", "GM_Ura");
- break;
- case NEPTUNE :
- rawGM = getLoadedConstant("GM8", "GM_Nep");
- break;
- case PLUTO :
- rawGM = getLoadedConstant("GM9", "GM_Plu");
- break;
- default :
- throw new OrekitInternalError(null);
- }
- final double au = getLoadedAstronomicalUnit();
- return rawGM * au * au * au / (Constants.JULIAN_DAY * Constants.JULIAN_DAY);
- }
- /** Get a constant defined in the ephemerides headers.
- * <p>Note that since constants are defined in the JPL headers
- * files, they are available as soon as one file is available, even
- * if it doesn't match the desired central date. This is because the
- * header must be parsed before the dates can be checked.</p>
- * <p>
- * There are alternate names for constants since for example JPL names are
- * different from INPOP names (Sun gravity: GMS or GM_Sun, Mars gravity:
- * GM4 or GM_Mar...).
- * </p>
- * @param names alternate names of the constant
- * @return value of the constant of NaN if the constant is not defined
- */
- public double getLoadedConstant(final String... names) {
- // lazy loading of constants
- Map<String, Double> map = constants.get();
- if (map == null) {
- final ConstantsParser parser = new ConstantsParser();
- if (!feed(parser)) {
- throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
- }
- map = parser.getConstants();
- constants.compareAndSet(null, map);
- }
- for (final String name : names) {
- if (map.containsKey(name)) {
- return map.get(name);
- }
- }
- return Double.NaN;
- }
- /** Get the maximal chunks duration.
- * @return chunks maximal duration in seconds
- */
- public double getMaxChunksDuration() {
- return maxChunksDuration;
- }
- /** Parse the first header record.
- * @param record first header record
- * @param name name of the file (or zip entry)
- */
- private void parseFirstHeaderRecord(final byte[] record, final String name) {
- // get the ephemerides type
- final int deNum = extractInt(record, HEADER_EPHEMERIS_TYPE_OFFSET);
- // as default, 3 polynomial coefficients for the Cartesian coordinates
- // (x, y, z) are contained in the file, positions are in kilometers
- // and times are in TDB
- components = 3;
- positionUnit = UnitsConverter.KILOMETRES_TO_METRES.convert(1.0);
- timeScale = timeScales.getTDB();
- if (deNum == INPOP_DE_NUMBER) {
- // an INPOP file may contain 6 components (including coefficients for the velocity vector)
- final double format = getLoadedConstant("FORMAT");
- if (!Double.isNaN(format) && (int) FastMath.IEEEremainder(format, 10) != 1) {
- components = 6;
- }
- // INPOP files may have their polynomials expressed in AU
- final double unite = getLoadedConstant("UNITE");
- if (!Double.isNaN(unite) && (int) unite == 0) {
- positionUnit = getLoadedAstronomicalUnit();
- }
- // INPOP files may have their times expressed in TCB
- final double timesc = getLoadedConstant("TIMESC");
- if (!Double.isNaN(timesc) && (int) timesc == 1) {
- timeScale = timeScales.getTCB();
- }
- }
- // extract covered date range
- startEpoch = extractDate(record, HEADER_START_EPOCH_OFFSET);
- finalEpoch = extractDate(record, HEADER_END_EPOCH_OFFSET);
- boolean ok = finalEpoch.compareTo(startEpoch) > 0;
- // indices of the Chebyshev coefficients for each ephemeris
- for (int i = 0; i < 12; ++i) {
- final int row1 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 12 * i);
- final int row2 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 4 + 12 * i);
- final int row3 = extractInt(record, HEADER_CHEBISHEV_INDICES_OFFSET + 8 + 12 * i);
- ok = ok && row1 >= 0 && row2 >= 0 && row3 >= 0;
- if (i == 0 && loadType == EphemerisType.MERCURY ||
- i == 1 && loadType == EphemerisType.VENUS ||
- i == 2 && loadType == EphemerisType.EARTH_MOON ||
- i == 3 && loadType == EphemerisType.MARS ||
- i == 4 && loadType == EphemerisType.JUPITER ||
- i == 5 && loadType == EphemerisType.SATURN ||
- i == 6 && loadType == EphemerisType.URANUS ||
- i == 7 && loadType == EphemerisType.NEPTUNE ||
- i == 8 && loadType == EphemerisType.PLUTO ||
- i == 9 && loadType == EphemerisType.MOON ||
- i == 10 && loadType == EphemerisType.SUN) {
- firstIndex = row1;
- coeffs = row2;
- chunks = row3;
- }
- }
- // compute chunks duration
- final double timeSpan = extractDouble(record, HEADER_CHUNK_DURATION_OFFSET);
- ok = ok && timeSpan > 0 && timeSpan < 100;
- chunksDuration = Constants.JULIAN_DAY * (timeSpan / chunks);
- if (Double.isNaN(maxChunksDuration)) {
- maxChunksDuration = chunksDuration;
- } else {
- maxChunksDuration = FastMath.max(maxChunksDuration, chunksDuration);
- }
- // sanity checks
- if (!ok) {
- throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
- }
- }
- /** Read first header record.
- * @param input input stream
- * @param name name of the file (or zip entry)
- * @return record record where to put bytes
- * @exception IOException if a read error occurs
- */
- private byte[] readFirstRecord(final InputStream input, final String name)
- throws IOException {
- // read first part of record, up to the ephemeris type
- final byte[] firstPart = new byte[HEADER_RECORD_SIZE_OFFSET + 4];
- if (!readInRecord(input, firstPart, 0)) {
- throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
- }
- // detect the endian format
- detectEndianess(firstPart);
- // get the ephemerides type
- final int deNum = extractInt(firstPart, HEADER_EPHEMERIS_TYPE_OFFSET);
- // the record size for this file
- final int recordSize;
- if (deNum == INPOP_DE_NUMBER) {
- // INPOP files have an extended DE format, which includes also the record size
- recordSize = extractInt(firstPart, HEADER_RECORD_SIZE_OFFSET) << 3;
- } else {
- // compute the record size for original JPL files
- recordSize = computeRecordSize(firstPart, name);
- }
- if (recordSize <= 0) {
- throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
- }
- // build a record with the proper size and finish read of the first complete record
- final int start = firstPart.length;
- final byte[] record = new byte[recordSize];
- System.arraycopy(firstPart, 0, record, 0, firstPart.length);
- if (!readInRecord(input, record, start)) {
- throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
- }
- return record;
- }
- /** Parse constants from first two header records.
- * @param first first header record
- * @param second second header record
- * @return map of parsed constants
- */
- private Map<String, Double> parseConstants(final byte[] first, final byte[] second) {
- final Map<String, Double> map = new HashMap<>();
- for (int i = 0; i < CONSTANTS_MAX_NUMBER; ++i) {
- // Note: for extracting the strings from the binary file, it makes no difference
- // if the file is stored in big-endian or little-endian notation
- final String constantName = extractString(first, HEADER_CONSTANTS_NAMES_OFFSET + i * 6, 6);
- if (constantName.length() == 0) {
- // no more constants to read
- break;
- }
- final double constantValue = extractDouble(second, HEADER_CONSTANTS_VALUES_OFFSET + 8 * i);
- map.put(constantName, constantValue);
- }
- // INPOP files do not have constants for AU and EMRAT, thus extract them from
- // the header record and create a constant for them to be consistent with JPL files
- if (!map.containsKey(CONSTANT_AU)) {
- map.put(CONSTANT_AU, extractDouble(first, HEADER_ASTRONOMICAL_UNIT_OFFSET));
- }
- if (!map.containsKey(CONSTANT_EMRAT)) {
- map.put(CONSTANT_EMRAT, extractDouble(first, HEADER_EM_RATIO_OFFSET));
- }
- return map;
- }
- /** Read bytes into the current record array.
- * @param input input stream
- * @param record record where to put bytes
- * @param start start index where to put bytes
- * @return true if record has been filled up
- * @exception IOException if a read error occurs
- */
- private boolean readInRecord(final InputStream input, final byte[] record, final int start)
- throws IOException {
- int index = start;
- while (index != record.length) {
- final int n = input.read(record, index, record.length - index);
- if (n < 0) {
- return false;
- }
- index += n;
- }
- return true;
- }
- /** Detect whether the JPL ephemerides file is stored in big-endian or
- * little-endian notation.
- * @param record the array containing the binary JPL header
- */
- private void detectEndianess(final byte[] record) {
- // default to big-endian
- bigEndian = true;
- // first try to read the DE number in big-endian format
- // the number is stored as unsigned int, so we have to convert it properly
- final long deNum = extractInt(record, HEADER_EPHEMERIS_TYPE_OFFSET) & 0xffffffffL;
- // simple heuristic: if the read value is larger than half the range of an integer
- // assume the file is in little-endian format
- if (deNum > (1 << 15)) {
- bigEndian = false;
- }
- }
- /** Calculate the record size of a JPL ephemerides file.
- * @param record the byte array containing the header record
- * @param name the name of the data file
- * @return the record size for this file
- */
- private int computeRecordSize(final byte[] record, final String name) {
- int recordSize = 0;
- boolean ok = true;
- // JPL files always have 3 position components
- final int nComp = 3;
- // iterate over the coefficient ptr array and sum up the record size
- // the coeffPtr array has the dimensions [12][nComp]
- for (int j = 0; j < 12; j++) {
- final int nCompCur = (j == 11) ? 2 : nComp;
- // Note: the array element coeffPtr[j][0] is not needed for the calculation
- final int idx = HEADER_CHEBISHEV_INDICES_OFFSET + j * nComp * 4;
- final int coeffPtr1 = extractInt(record, idx + 4);
- final int coeffPtr2 = extractInt(record, idx + 8);
- // sanity checks
- ok = ok && (coeffPtr1 >= 0 || coeffPtr2 >= 0);
- recordSize += coeffPtr1 * coeffPtr2 * nCompCur;
- }
- // the libration ptr array has the dimension [3]
- // Note: the array element libratPtr[0] is not needed for the calculation
- final int libratPtr1 = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 4);
- final int libratPtr2 = extractInt(record, HEADER_LIBRATION_INDICES_OFFSET + 8);
- // sanity checks
- ok = ok && (libratPtr1 >= 0 || libratPtr2 >= 0);
- recordSize += libratPtr1 * libratPtr2 * nComp + 2;
- recordSize <<= 3;
- if (!ok || recordSize <= 0) {
- throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
- }
- return recordSize;
- }
- /** Extract a date from a record.
- * @param record record to parse
- * @param offset offset of the double within the record
- * @return extracted date
- */
- private AbsoluteDate extractDate(final byte[] record, final int offset) {
- final double t = extractDouble(record, offset);
- int jDay = (int) FastMath.floor(t);
- double seconds = (t + 0.5 - jDay) * Constants.JULIAN_DAY;
- if (seconds >= Constants.JULIAN_DAY) {
- ++jDay;
- seconds -= Constants.JULIAN_DAY;
- }
- return new AbsoluteDate(new DateComponents(DateComponents.JULIAN_EPOCH, jDay),
- new TimeComponents(seconds), timeScale);
- }
- /** Extract a double from a record.
- * <p>Double numbers are stored according to IEEE 754 standard, with
- * most significant byte first.</p>
- * @param record record to parse
- * @param offset offset of the double within the record
- * @return extracted double
- */
- private double extractDouble(final byte[] record, final int offset) {
- final long l8 = ((long) record[offset + 0]) & 0xffl;
- final long l7 = ((long) record[offset + 1]) & 0xffl;
- final long l6 = ((long) record[offset + 2]) & 0xffl;
- final long l5 = ((long) record[offset + 3]) & 0xffl;
- final long l4 = ((long) record[offset + 4]) & 0xffl;
- final long l3 = ((long) record[offset + 5]) & 0xffl;
- final long l2 = ((long) record[offset + 6]) & 0xffl;
- final long l1 = ((long) record[offset + 7]) & 0xffl;
- final long l;
- if (bigEndian) {
- l = (l8 << 56) | (l7 << 48) | (l6 << 40) | (l5 << 32) |
- (l4 << 24) | (l3 << 16) | (l2 << 8) | l1;
- } else {
- l = (l1 << 56) | (l2 << 48) | (l3 << 40) | (l4 << 32) |
- (l5 << 24) | (l6 << 16) | (l7 << 8) | l8;
- }
- return Double.longBitsToDouble(l);
- }
- /** Extract an int from a record.
- * @param record record to parse
- * @param offset offset of the double within the record
- * @return extracted int
- */
- private int extractInt(final byte[] record, final int offset) {
- final int l4 = ((int) record[offset + 0]) & 0xff;
- final int l3 = ((int) record[offset + 1]) & 0xff;
- final int l2 = ((int) record[offset + 2]) & 0xff;
- final int l1 = ((int) record[offset + 3]) & 0xff;
- if (bigEndian) {
- return (l4 << 24) | (l3 << 16) | (l2 << 8) | l1;
- } else {
- return (l1 << 24) | (l2 << 16) | (l3 << 8) | l4;
- }
- }
- /** Extract a String from a record.
- * @param record record to parse
- * @param offset offset of the string within the record
- * @param length maximal length of the string
- * @return extracted string, with whitespace characters stripped
- */
- private String extractString(final byte[] record, final int offset, final int length) {
- return new String(record, offset, length, StandardCharsets.US_ASCII).trim();
- }
- /** Local parser for header constants. */
- private class ConstantsParser implements DataLoader {
- /** Local constants map. */
- private Map<String, Double> localConstants;
- /** Get the local constants map.
- * @return local constants map
- */
- public Map<String, Double> getConstants() {
- return localConstants;
- }
- /** {@inheritDoc} */
- public boolean stillAcceptsData() {
- return localConstants == null;
- }
- /** {@inheritDoc} */
- public void loadData(final InputStream input, final String name)
- throws IOException, ParseException, OrekitException {
- // read first header record
- final byte[] first = readFirstRecord(input, name);
- // the second record contains the values of the constants used for least-square filtering
- final byte[] second = new byte[first.length];
- if (!readInRecord(input, second, 0)) {
- throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
- }
- localConstants = parseConstants(first, second);
- }
- }
- /** Local parser for Chebyshev polynomials. */
- private class EphemerisParser implements DataLoader, TimeStampedGenerator<PosVelChebyshev> {
- /** Set of Chebyshev polynomials read. */
- private final SortedSet<PosVelChebyshev> entries;
- /** Start of range we are interested in. */
- private AbsoluteDate start;
- /** End of range we are interested in. */
- private AbsoluteDate end;
- /** Simple constructor.
- */
- EphemerisParser() {
- entries = new TreeSet<>(new ChronologicalComparator());
- }
- /** {@inheritDoc} */
- public List<PosVelChebyshev> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {
- try {
- // prepare reading
- entries.clear();
- if (existingDate == null) {
- // we want ephemeris data for the first time, set up an arbitrary first range
- start = date.shiftedBy(-FIFTY_DAYS);
- end = date.shiftedBy(+FIFTY_DAYS);
- } else if (existingDate.compareTo(date) <= 0) {
- // we want to extend an existing range towards future dates
- start = existingDate;
- end = date;
- } else {
- // we want to extend an existing range towards past dates
- start = date;
- end = existingDate;
- }
- // get new entries in the specified data range
- if (!feed(this)) {
- throw new OrekitException(OrekitMessages.NO_JPL_EPHEMERIDES_BINARY_FILES_FOUND);
- }
- return new ArrayList<>(entries);
- } catch (OrekitException oe) {
- throw new TimeStampedCacheException(oe);
- }
- }
- /** {@inheritDoc} */
- public boolean stillAcceptsData() {
- // special case for Earth: we do not really load any ephemeris data
- if (generateType == EphemerisType.EARTH) {
- return false;
- }
- // we have to look for data in all available ephemerides files as there may be
- // data overlaps that result in incomplete data
- if (entries.isEmpty()) {
- return true;
- } else {
- // if the requested range is already filled, we do not need to look further
- return !(entries.first().getDate().compareTo(start) < 0 &&
- entries.last().getDate().compareTo(end) > 0);
- }
- }
- /** {@inheritDoc} */
- public void loadData(final InputStream input, final String name)
- throws IOException {
- // read first header record
- final byte[] first = readFirstRecord(input, name);
- // the second record contains the values of the constants used for least-square filtering
- final byte[] second = new byte[first.length];
- if (!readInRecord(input, second, 0)) {
- throw new OrekitException(OrekitMessages.UNABLE_TO_READ_JPL_HEADER, name);
- }
- if (constants.get() == null) {
- constants.compareAndSet(null, parseConstants(first, second));
- }
- // check astronomical unit consistency
- final double au = 1000 * extractDouble(first, HEADER_ASTRONOMICAL_UNIT_OFFSET);
- if (au < 1.4e11 || au > 1.6e11) {
- throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
- }
- if (FastMath.abs(getLoadedAstronomicalUnit() - au) >= 10.0) {
- throw new OrekitException(OrekitMessages.INCONSISTENT_ASTRONOMICAL_UNIT_IN_FILES,
- getLoadedAstronomicalUnit(), au);
- }
- // check Earth-Moon mass ratio consistency
- final double emRat = extractDouble(first, HEADER_EM_RATIO_OFFSET);
- if (emRat < 80 || emRat > 82) {
- throw new OrekitException(OrekitMessages.NOT_A_JPL_EPHEMERIDES_BINARY_FILE, name);
- }
- if (FastMath.abs(getLoadedEarthMoonMassRatio() - emRat) >= 1.0e-5) {
- throw new OrekitException(OrekitMessages.INCONSISTENT_EARTH_MOON_RATIO_IN_FILES,
- getLoadedEarthMoonMassRatio(), emRat);
- }
- // parse first header record
- parseFirstHeaderRecord(first, name);
- if (startEpoch.compareTo(end) < 0 && finalEpoch.compareTo(start) > 0) {
- // this file contains data in the range we are looking for, read it
- final byte[] record = new byte[first.length];
- while (readInRecord(input, record, 0)) {
- final AbsoluteDate rangeStart = parseDataRecord(record);
- if (rangeStart.compareTo(end) > 0) {
- // we have already exceeded the range we were interested in,
- // we interrupt parsing here
- return;
- }
- }
- }
- }
- /** Parse regular ephemeris record.
- * @param record record to parse
- * @return date of the last parsed chunk
- */
- private AbsoluteDate parseDataRecord(final byte[] record) {
- // extract time range covered by the record
- final AbsoluteDate rangeStart = extractDate(record, DATA_START_RANGE_OFFSET);
- if (rangeStart.compareTo(startEpoch) < 0) {
- throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE,
- rangeStart, startEpoch, finalEpoch, startEpoch.durationFrom(rangeStart));
- }
- final AbsoluteDate rangeEnd = extractDate(record, DATE_END_RANGE_OFFSET);
- if (rangeEnd.compareTo(finalEpoch) > 0) {
- throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_AFTER,
- rangeEnd, startEpoch, finalEpoch, rangeEnd.durationFrom(finalEpoch));
- }
- if (rangeStart.compareTo(end) > 0 || rangeEnd.compareTo(start) < 0) {
- // we are not interested in this record, don't parse it
- return rangeEnd;
- }
- // loop over chunks inside the time range
- AbsoluteDate chunkEnd = rangeStart;
- final int nbChunks = chunks;
- final int nbCoeffs = coeffs;
- final int first = firstIndex;
- final double duration = chunksDuration;
- for (int i = 0; i < nbChunks; ++i) {
- // set up chunk validity range
- final AbsoluteDate chunkStart = chunkEnd;
- chunkEnd = (i == nbChunks - 1) ? rangeEnd : rangeStart.shiftedBy((i + 1) * duration);
- // extract Chebyshev coefficients for the selected body
- // and convert them from kilometers to meters
- final double[] xCoeffs = new double[nbCoeffs];
- final double[] yCoeffs = new double[nbCoeffs];
- final double[] zCoeffs = new double[nbCoeffs];
- for (int k = 0; k < nbCoeffs; ++k) {
- // by now, only use the position components
- // if there are also velocity components contained in the file, ignore them
- final int index = first + components * i * nbCoeffs + k - 1;
- xCoeffs[k] = positionUnit * extractDouble(record, 8 * index);
- yCoeffs[k] = positionUnit * extractDouble(record, 8 * (index + nbCoeffs));
- zCoeffs[k] = positionUnit * extractDouble(record, 8 * (index + 2 * nbCoeffs));
- }
- // build the position-velocity model for current chunk
- entries.add(new PosVelChebyshev(chunkStart, timeScale, duration, xCoeffs, yCoeffs, zCoeffs));
- }
- return rangeStart;
- }
- }
- /** Raw position-velocity provider using ephemeris. */
- private class EphemerisRawPVProvider implements RawPVProvider {
- /** Retrieve correct Chebyshev polynomial for given date.
- * @param date date
- * @return PosVelChebyshev Chebyshev polynomial class for position-velocity-acceleration
- */
- private PosVelChebyshev getChebyshev(final AbsoluteDate date) {
- PosVelChebyshev chebyshev;
- try {
- chebyshev = ephemerides.getNeighbors(date).findFirst().get();
- } catch (TimeStampedCacheException tce) {
- // we cannot bracket the date, check if the last available chunk covers the specified date
- chebyshev = ephemerides.getLatest();
- if (!chebyshev.inRange(date)) {
- // we were not able to recover from the error, the date is too far
- throw tce;
- }
- }
- return chebyshev;
- }
- /** {@inheritDoc} */
- public PVCoordinates getRawPV(final AbsoluteDate date) {
- // get raw PV from Chebyshev polynomials
- final PosVelChebyshev chebyshev = getChebyshev(date);
- // evaluate the Chebyshev polynomials
- return chebyshev.getPositionVelocityAcceleration(date);
- }
- /** {@inheritDoc} */
- public <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {
- // get raw PV from Chebyshev polynomials
- final PosVelChebyshev chebyshev = getChebyshev(date.toAbsoluteDate());
- // evaluate the Chebyshev polynomials
- return chebyshev.getPositionVelocityAcceleration(date);
- }
- /** {@inheritDoc} */
- @Override
- public Vector3D getRawPosition(final AbsoluteDate date) {
- // get raw PV from Chebyshev polynomials
- final PosVelChebyshev chebyshev = getChebyshev(date);
- // evaluate the Chebyshev polynomials
- return chebyshev.getPosition(date);
- }
- /** {@inheritDoc} */
- @Override
- public <T extends CalculusFieldElement<T>> FieldVector3D<T> getRawPosition(final FieldAbsoluteDate<T> date) {
- // get raw PV from Chebyshev polynomials
- final PosVelChebyshev chebyshev = getChebyshev(date.toAbsoluteDate());
- // evaluate the Chebyshev polynomials
- return chebyshev.getPosition(date);
- }
- }
- /** Raw position-velocity provider providing always zero. */
- private static class ZeroRawPVProvider implements RawPVProvider {
- /** {@inheritDoc} */
- public PVCoordinates getRawPV(final AbsoluteDate date) {
- return PVCoordinates.ZERO;
- }
- /** {@inheritDoc} */
- public <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> getRawPV(final FieldAbsoluteDate<T> date) {
- return FieldPVCoordinates.getZero(date.getField());
- }
- }
- }