The current celestia code supports only JPL Big Endian ephemerides
DE406/DE405/DE200.
The big endian ephemeris DE406 comes in chunk of 300 years. If you
want to concatenate several pieces, using the software provided by JPL
(which does not handle endianess), you need a big endian platform. The
little endian DE406 comes in one piece of 6000 years, so you don't
even need to concatenate this one.
The patch proposed adds the following:
- support for both endianess with autodetection of endianess
- support for all JPL DE ephemeris (from DE200) with autodetection of
record size.
- support for JPL DE compatible ephemeris such as the INPOP ephemeris
from Paris observatory (only the ephemeris files in the legacy JPL DE
compatible format and with timescale TDB, not the INPOP native format,
or the TCB files). http://www.imcce.fr/inpop/download10e.php
The method to autodetect ephemeris type, record size and endianess is
based on the CALCEPH software.
It would probably be better if possible to wrap around the CALCEPH
library (so that it would also support INPOP native format) but this
would be more work.
In addition to the ephemeris number and the start and end date, the
record size and endianess are also printed to the log (which you get
by pressing the tilde key).
If the ephemeris is an INPOP ephemeris, then the ephemeris number is
not printed to the log: it just says "INPOP" (I have been too lazy to
find out where this is stored).
It has been tested with:
- DE406 big endian
- DE406 little endian
- DE405 big endian
- DE423 little endian
- INPOP10A little endian
However, I have not performed any accurate regression testing. I have
just checked that the planets orbits do not look strange and that the
log is reporting the correct endianess, record size and ephemeris
name.
Only tested on Linux/x86_64, would appreciate if it could be tested
further and eventually be merged.
Code: Select all
Index: src/celephem/jpleph.cpp
===================================================================
--- src/celephem/jpleph.cpp (revision 5224)
+++ src/celephem/jpleph.cpp (working copy)
@@ -21,9 +21,9 @@
using namespace Eigen;
using namespace std;
-static const unsigned int DE200RecordSize = 826;
-static const unsigned int DE405RecordSize = 1018;
-static const unsigned int DE406RecordSize = 728;
+//static const unsigned int DE200RecordSize = 826;
+//static const unsigned int DE405RecordSize = 1018;
+//static const unsigned int DE406RecordSize = 728;
static const unsigned int NConstants = 400;
static const unsigned int ConstantNameLength = 6;
@@ -33,22 +33,22 @@
static const int LabelSize = 84;
-// Read a big-endian 32-bit unsigned integer
-static int32 readUint(istream& in)
+// Read a big-endian or little endian 32-bit unsigned integer
+static int32 readUint(istream& in, bool swap)
{
int32 ret;
in.read((char*) &ret, sizeof(int32));
- BE_TO_CPU_INT32(ret, ret);
+ if (swap) ret= bswap_32(ret);
return (uint32) ret;
}
-// Read a big-endian 64-bit IEEE double--if the native double format isn't
+// Read a big-endian or little endian 64-bit IEEE double--if the native double format isn't
// IEEE 754, there will be troubles.
-static double readDouble(istream& in)
+static double readDouble(istream& in, bool swap)
{
double d;
in.read((char*) &d, sizeof(double));
- BE_TO_CPU_DOUBLE(d, d);
+ if (swap) d = bswap_double(d);
return d;
}
@@ -86,7 +86,16 @@
return endDate;
}
+unsigned int JPLEphemeris::getRecordSize() const
+{
+ return recordSize;
+}
+bool JPLEphemeris::getByteSwap() const
+{
+ return swapBytes;
+}
+
// Return the position of an object relative to the solar system barycenter
// or the Earth (in the case of the Moon) at a specified TDB Julian date tjd.
// If tjd is outside the span covered by the ephemeris it is clamped to a
@@ -171,26 +180,87 @@
JPLEphemeris* JPLEphemeris::load(istream& in)
{
+
JPLEphemeris* eph = NULL;
+ eph = new JPLEphemeris();
+ if (eph == NULL)
+ return NULL;
+ // Figure out ephemeris type and endianess
+
// Skip past three header labels
in.ignore(LabelSize * 3);
if (!in.good())
return NULL;
-
// Skip past the constant names
in.ignore(NConstants * ConstantNameLength);
if (!in.good())
return NULL;
+ if (!in.good())
+ return NULL;
+ // Skip past the start time, end time, and time interval
+ in.ignore(3 * 8);
+ if (!in.good())
+ return NULL;
+ // Skip past number of constants with valid values
+ in.ignore(4);
+ if (!in.good())
+ return NULL;
+ // Skip past AU and Earth-Moon ratio
+ in.ignore(2 * 8);
+ if (!in.good())
+ return NULL;
+ // Skip past coefficient information for each item in the ephemeris
+ in.ignore(4 * 3 * JPLEph_NItems);
+ if (!in.good())
+ return NULL;
- eph = new JPLEphemeris();
- if (eph == NULL)
+ // read DE number
+ uint32 deNum = readUint(in,false);
+ uint32 deNum2 = (uint32) bswap_32(deNum);
+
+ if (deNum==100) {
+ // INPOP ephemeris with same endianess as CPU
+ eph->swapBytes=false;
+ eph->DENum=deNum;
+ }
+ else if (deNum2==100) {
+ // INPOP ephemeris with different endianess
+ eph->swapBytes=true;
+ eph->DENum=deNum2;
+ }
+ else if ((deNum>(1U<<15)) and (deNum2>=200)){
+ // DE ephemeris with different endianess
+ eph->swapBytes=true;
+ eph->DENum=deNum2;
+ }
+ else if ((deNum<=(1U<<15)) and (deNum>=200)){
+ // DE ephemeris with same endianess as CPU
+ eph->swapBytes=false;
+ eph->DENum=deNum;
+ }
+ else {
+ delete eph;
return NULL;
+ }
+
+ // Rewind input file
+ in.seekg(0);
+ // Skip past three header labels
+ in.ignore(LabelSize * 3);
+ if (!in.good())
+ return NULL;
+
+ // Skip past the constant names
+ in.ignore(NConstants * ConstantNameLength);
+ if (!in.good())
+ return NULL;
+
// Read the start time, end time, and time interval
- eph->startDate = readDouble(in);
- eph->endDate = readDouble(in);
- eph->daysPerInterval = readDouble(in);
+ eph->startDate = readDouble(in,eph->swapBytes);
+ eph->endDate = readDouble(in,eph->swapBytes);
+ eph->daysPerInterval = readDouble(in,eph->swapBytes);
if (!in.good())
{
delete eph;
@@ -198,54 +268,52 @@
}
// Number of constants with valid values; not useful for us
- (void) readUint(in);
+ (void) readUint(in,eph->swapBytes);
- eph->au = readDouble(in); // kilometers per astronomical unit
- eph->earthMoonMassRatio = readDouble(in);
+ eph->au = readDouble(in,eph->swapBytes); // kilometers per astronomical unit
+ eph->earthMoonMassRatio = readDouble(in,eph->swapBytes);
// Read the coefficient information for each item in the ephemeris
unsigned int i;
+ eph->recordSize=0;
for (i = 0; i < JPLEph_NItems; i++)
{
- eph->coeffInfo[i].offset = readUint(in) - 3;
- eph->coeffInfo[i].nCoeffs = readUint(in);
- eph->coeffInfo[i].nGranules = readUint(in);
+ eph->coeffInfo[i].offset = readUint(in,eph->swapBytes) - 3;
+ eph->coeffInfo[i].nCoeffs = readUint(in,eph->swapBytes);
+ eph->coeffInfo[i].nGranules = readUint(in,eph->swapBytes);
+ eph->recordSize+=eph->coeffInfo[i].nCoeffs*eph->coeffInfo[i].nGranules*((i==11)?2:3); // last item is the nutation ephemeris (only 2 components)
}
if (!in.good())
{
delete eph;
return NULL;
}
+ // Skip DE number
+ in.ignore(4);
- eph->DENum = readUint(in);
+ eph->librationCoeffInfo.offset = readUint(in,eph->swapBytes);
+ eph->librationCoeffInfo.nCoeffs = readUint(in,eph->swapBytes);
+ eph->librationCoeffInfo.nGranules = readUint(in,eph->swapBytes);
+ eph->recordSize+=eph->librationCoeffInfo.nCoeffs*eph->librationCoeffInfo.nGranules*3;
+ eph->recordSize+=2; // record start and end time
- switch (eph->DENum)
+ if (!in.good())
{
- case 200:
- eph->recordSize = DE200RecordSize;
- break;
- case 405:
- eph->recordSize = DE405RecordSize;
- break;
- case 406:
- eph->recordSize = DE406RecordSize;
- break;
- default:
delete eph;
return NULL;
}
- eph->librationCoeffInfo.offset = readUint(in);
- eph->librationCoeffInfo.nCoeffs = readUint(in);
- eph->librationCoeffInfo.nGranules = readUint(in);
- if (!in.good())
- {
- delete eph;
- return NULL;
+ // if INPOP ephemeris, read record size
+ if (deNum==100) {
+ eph->recordSize=readUint(in,eph->swapBytes);
+ // Skip past the rest of the record
+ in.ignore(eph->recordSize * 8 - 2860);
}
+ else {
+ // Skip past the rest of the record
+ in.ignore(eph->recordSize * 8 - 2856);
+ }
- // Skip past the rest of the record
- in.ignore(eph->recordSize * 8 - 2856);
// The next record contains constant values (which we don't need)
in.ignore(eph->recordSize * 8);
if (!in.good())
@@ -259,14 +327,14 @@
eph->records.resize(nRecords);
for (i = 0; i < nRecords; i++)
{
- eph->records[i].t0 = readDouble(in);
- eph->records[i].t1 = readDouble(in);
+ eph->records[i].t0 = readDouble(in,eph->swapBytes);
+ eph->records[i].t1 = readDouble(in,eph->swapBytes);
// Allocate coefficient array for this record; the first two
// 'coefficients' are actually the start and end time (t0 and t1)
eph->records[i].coeffs = new double[eph->recordSize - 2];
for (unsigned int j = 0; j < eph->recordSize - 2; j++)
- eph->records[i].coeffs[j] = readDouble(in);
+ eph->records[i].coeffs[j] = readDouble(in,eph->swapBytes);
// Make sure that we read this record successfully
if (!in.good())
Index: src/celephem/customorbit.cpp
===================================================================
--- src/celephem/customorbit.cpp (revision 5224)
+++ src/celephem/customorbit.cpp (working copy)
@@ -3160,12 +3160,20 @@
if (in.good())
jpleph = JPLEphemeris::load(in);
if (jpleph != NULL)
- {
- clog << "Loaded DE" << jpleph->getDENumber() <<
- " ephemeris. Valid from JD" <<
+ {
+ if (jpleph->getDENumber()!=100) {
+ clog << "Loaded DE" << jpleph->getDENumber();
+ }
+ else {
+ clog << "Loaded INPOP";
+ }
+ clog <<" ephemeris. Valid from JD" <<
setprecision(8) <<
jpleph->getStartDate() << " to JD" <<
jpleph->getEndDate() << '\n';
+ clog << "Ephemeris record size: "<< jpleph->getRecordSize()<< " doubles";
+ if (jpleph->getByteSwap()) clog << " with non-native endianess";
+ clog << "." << endl;
}
}
Index: src/celephem/jpleph.h
===================================================================
--- src/celephem/jpleph.h (revision 5224)
+++ src/celephem/jpleph.h (working copy)
@@ -59,6 +59,7 @@
class JPLEphemeris
{
private:
+
JPLEphemeris();
public:
@@ -71,6 +72,8 @@
unsigned int getDENumber() const;
double getStartDate() const;
double getEndDate() const;
+ bool getByteSwap() const;
+ unsigned int getRecordSize() const;
private:
JPLEphCoeffInfo coeffInfo[JPLEph_NItems];
@@ -85,7 +88,7 @@
unsigned int DENum; // ephemeris version
unsigned int recordSize; // number of doubles per record
-
+ bool swapBytes;
std::vector<JPLEphRecord> records;
};
