JPLEphemeridesLoader.java
/* Copyright 2002-2023 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());
}
}
}