Skip to content

Commit

Permalink
Time tests passing
Browse files Browse the repository at this point in the history
  • Loading branch information
@gcasadesus
gcasadesus committed Dec 29, 2024
1 parent cd8f3f0 commit d75be00
Show file tree
Hide file tree
Showing 9 changed files with 143 additions and 79 deletions.
6 changes: 3 additions & 3 deletions include/lupnt/core/constants.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,9 +71,9 @@ namespace lupnt {
static constexpr double MJD_J2000_TT = 51544.5; // [days]

// Vallado page 194
static constexpr double JD_COORDINATE_TAI = 2443144.5; // [days]
static constexpr double JD_COORDINATE_TT_TCG_TCB
= JD_COORDINATE_TAI + TT_TAI_OFFSET / SECS_DAY; // [days]
static constexpr double MJD_COORDINATE_TAI = 2443144.5 - JD_MJD_OFFSET; // [days]
static constexpr double MJD_COORDINATE_TT_TCG_TCB
= MJD_COORDINATE_TAI + TT_TAI_OFFSET / SECS_DAY; // [days]

static constexpr double L_B = 1.550505e-8;
static constexpr double L_G = 6.969290134e-10;
Expand Down
4 changes: 4 additions & 0 deletions include/lupnt/physics/time_converter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,9 @@

namespace lupnt {

extern const std::map<Time, std::string> time2string;
extern const std::map<std::string, Time> string2time;

Real ConvertTime(Real t, Time from, Time to);
VecX ConvertTime(VecX t, Time from, Time to);

Expand Down Expand Up @@ -42,6 +45,7 @@ namespace lupnt {
Real EarthRotationAngle(Real t_ut1);
Real Gregorian2MJD(int year, int month, int day, int hour = 0, int min = 0, Real sec = 0);
Real Gregorian2Time(int year, int month, int day, int hour = 0, int min = 0, Real sec = 0);
Real Gregorian2Time(const std::string& date);

Real GreenwichMeanSiderealTime(Real mjd_ut1);
Real GreenwichApparentSiderealTime(Real mjd_ut1);
Expand Down
6 changes: 0 additions & 6 deletions source/cpp/physics/spice_interface.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,12 +31,6 @@ namespace lupnt {
namespace spice {
bool spice_loaded = false;

const std::map<Time, std::string> time2string = {
{Time::UT1, "UT1"}, {Time::UTC, "UTC"}, {Time::TAI, "TAI"}, {Time::TDB, "TDB"},
{Time::TT, "TT"}, {Time::TCG, "TCG"}, {Time::TCB, "TCB"}, {Time::GPS, "GPS"},
{Time::JD_TT, "JDTDT"}, {Time::JD_TDB, "JDTDB"},
};

/**
* @brief load the Spice kernels
*
Expand Down
52 changes: 36 additions & 16 deletions source/cpp/physics/time_converter.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,18 @@ namespace lupnt {
// GPS Epoch 1980-01-06T00:00:00 UTC
// J2000 Epoch 2000-01-01T12:00:00 TT

const std::map<Time, std::string> time2string = {
{Time::UT1, "UT1"}, {Time::UTC, "UTC"}, {Time::TAI, "TAI"}, {Time::TDB, "TDB"},
{Time::TT, "TT"}, {Time::TCG, "TCG"}, {Time::TCB, "TCB"}, {Time::GPS, "GPS"},
{Time::JD_TT, "JDTDT"}, {Time::JD_TDB, "JDTDB"},
};

const std::map<std::string, Time> string2time = {
{"UT1", Time::UT1}, {"UTC", Time::UTC}, {"TAI", Time::TAI}, {"TDB", Time::TDB},
{"TT", Time::TT}, {"TCG", Time::TCG}, {"TCB", Time::TCB}, {"GPS", Time::GPS},
{"JDTDT", Time::JD_TT}, {"JDTDB", Time::JD_TDB},
};

/// @brief Convert time from one time system to another
/// @param t Time in the original time system
/// @param from Original time system
Expand Down Expand Up @@ -139,14 +151,14 @@ namespace lupnt {
Real TT2TAI(Real t_tt) { return t_tt - TT_TAI_OFFSET; }

Real TT2TCG(Real t_tt) {
Real jd_tt = Time2JD(t_tt);
Real tt_tcg = -L_G / (1.0 - L_G) * (jd_tt - JD_COORDINATE_TT_TCG_TCB) * SECS_DAY;
Real mjd_tt = Time2MJD(t_tt);
Real tt_tcg = -L_G / (1.0 - L_G) * (mjd_tt - MJD_COORDINATE_TT_TCG_TCB) * SECS_DAY;
return t_tt - tt_tcg;
}

Real TCG2TT(Real t_tcg) {
Real jd_tcg = Time2JD(t_tcg);
Real tt_tcg = -L_G * (jd_tcg - JD_COORDINATE_TT_TCG_TCB) * SECS_DAY;
Real mjd_tcg = Time2MJD(t_tcg);
Real tt_tcg = -L_G * (mjd_tcg - MJD_COORDINATE_TT_TCG_TCB) * SECS_DAY;
return t_tcg + tt_tcg;
}

Expand All @@ -168,14 +180,12 @@ namespace lupnt {
/// @param t_tt
/// @return
/// @ref
/// https://gssc.esa.int/navipedia/index.php/Transformations_between_Time_Systems#TDT_-_TDB.2C_TCB
/// @note Accurate to about 30 microseconds
/// Astrodynamics Convention and Modeling Reference for Lunar, Cislunar, and Libration Point
/// Orbits
Real TT2TDB(Real t_tt) {
double k = 1.657e-3;
double eb = 1.671e-2;
Real mean_anom = 6.239996 + 1.99096871e-7 * t_tt;
Real ecc_anom = mean_anom + eb * sin(mean_anom);
Real t_tdb = t_tt + k * sin(ecc_anom);
Real days_j2000_tt = t_tt / SECS_DAY;
Real me = (357.53 + 0.9856003 * days_j2000_tt) * RAD;
Real t_tdb = t_tt + 0.001658 * sin(me) + 0.000014 * sin(2 * me);
return t_tdb;
}

Expand All @@ -191,15 +201,18 @@ namespace lupnt {

Real TCB2TDB(Real t_tcb) {
const double tdb0 = -6.55e-5;
Real jd_tcb = Time2JD(t_tcb);
Real t_tdb = t_tcb - (1.55051976772e-8 * (jd_tcb - JD_COORDINATE_TT_TCG_TCB) * SECS_DAY + tdb0);
Real mjd_tcb = Time2MJD(t_tcb);
Real t_tdb
= t_tcb - (1.55051976772e-8 * (mjd_tcb - MJD_COORDINATE_TT_TCG_TCB) * SECS_DAY + tdb0);
return t_tdb;
}

Real TT2TCB(Real t_tt) {
Real t_tai = TT2TAI(t_tt);
Real jd_tai = Time2JD(t_tai);
Real t_tcb = t_tt + 1.5505197677e-8 * (jd_tai - JD_COORDINATE_TAI) * SECS_DAY;
Real t_tdb = TT2TDB(t_tt);
const double L_B = 1.550519768e-8;
Real offset = t_tai + TT2TAI(0) + (MJD_J2000_TT - MJD_COORDINATE_TAI) * SECS_DAY;
Real t_tcb = t_tdb + L_B * offset;
return t_tcb;
}

Expand Down Expand Up @@ -255,6 +268,13 @@ namespace lupnt {
return MJD2Time(mjd);
}

Real Gregorian2Time(const std::string& date) {
int year, month, day, hour, min;
double sec;
sscanf(date.c_str(), "%d-%d-%dT%d:%d:%lf", &year, &month, &day, &hour, &min, &sec);
return Gregorian2Time(year, month, day, hour, min, sec);
}

/// @brief Greenwich Mean Sidereal Time
/// @param mjd_ut1 UT1 (Modified Julian Date)
/// @return GMST [rad]
Expand Down Expand Up @@ -292,7 +312,7 @@ namespace lupnt {
sec = round(sec, precision);
ss << year << "-";
ss << std::setw(2) << std::setfill('0') << month << "-";
ss << std::setw(2) << std::setfill('0') << day << " ";
ss << std::setw(2) << std::setfill('0') << day << "T";
ss << std::setw(2) << std::setfill('0') << hour << ":";
ss << std::setw(2) << std::setfill('0') << min << ":";
ss << std::setw(2) << std::setfill('0') << floor(sec);
Expand Down
9 changes: 6 additions & 3 deletions source/python/bindings/py_constants.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,14 +39,17 @@ void init_constants(py::module& m) {
m.attr("DAYS_CENTURY") = py::float_(DAYS_CENTURY);
m.attr("DAYS_SEC") = py::float_(DAYS_SEC);

m.attr("JD_MJD_OFFSET") = py::float_(JD_MJD_OFFSET);
m.attr("TT_TAI_OFFSET") = py::float_(TT_TAI_OFFSET);
m.attr("A1_TAI_OFFSET") = py::float_(A1_TAI_OFFSET);

m.attr("JD_CCSDS_TAI") = py::float_(JD_CCSDS_TAI);
m.attr("JD_J2000_TT") = py::float_(JD_J2000_TT);
m.attr("MJD_CCSDS_TAI") = py::float_(MJD_CCSDS_TAI);
m.attr("MJD_J2000_TT") = py::float_(MJD_J2000_TT);

m.attr("JD_MJD_OFFSET") = py::float_(JD_MJD_OFFSET);
m.attr("TT_TAI_OFFSET") = py::float_(TT_TAI_OFFSET);
m.attr("A1_TAI_OFFSET") = py::float_(A1_TAI_OFFSET);
m.attr("MJD_COORDINATE_TAI") = py::float_(MJD_COORDINATE_TAI);
m.attr("MJD_COORDINATE_TT_TCG_TCB") = py::float_(MJD_COORDINATE_TT_TCG_TCB);

m.attr("L_B") = py::float_(L_B);
m.attr("L_G") = py::float_(L_G);
Expand Down
8 changes: 4 additions & 4 deletions source/python/bindings/py_spice_interface.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,14 @@ namespace py = pybind11;

class SpiceInterface {}; // dummy class

void init_spice_interface(py::module &m) {
void init_spice_interface(py::module& m) {
auto m_spice = m.def_submodule("spice", "This is A.");

m_spice.def("load_spice_kernel", &lupnt::spice::LoadSpiceKernel)
.def("extract_pck_coeffs", &lupnt::spice::ExtractPckCoeffs)
.def(
"get_frame_conversion_mat",
[](double t_tai, lupnt::Frame from, lupnt::Frame to) -> lupnt::Mat6d {
[](double t_tai, const std::string& from, const std::string& to) -> lupnt::Mat6d {
return lupnt::spice::GetFrameConversionMat(t_tai, from, to).cast<double>();
},
py::arg("t_tai"), py::arg("from"), py::arg("to"))
Expand Down Expand Up @@ -51,8 +51,8 @@ void init_spice_interface(py::module &m) {
py::arg("t_tai"), py::arg("center"), py::arg("target"))
.def(
"get_body_pos_spice",
[](double t_tai, lupnt::NaifId obs, lupnt::NaifId target, lupnt::Frame ref_frame,
std::string ab_orrection) -> lupnt::Vec3d {
[](double t_tai, lupnt::NaifId obs, lupnt::NaifId target, const std::string& ref_frame,
const std::string& ab_orrection) -> lupnt::Vec3d {
return lupnt::spice::GetBodyPosSpice(t_tai, obs, target, ref_frame, ab_orrection)
.cast<double>();
},
Expand Down
24 changes: 19 additions & 5 deletions test/cpp/physics/test_time_converter.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,9 @@ TEST_CASE("physics.ConvertTime") {
// All-to-all time conversions
for (auto time_in : time_list) {
Real t_tai = Gregorian2Time(2022, 12, 18, 18, 0, 0.0);
Real t_in = ConvertTime(t_tai, Time::TAI, time_in);
Real t_diff_in = ConvertTime(t_tai, Time::TAI, time_in);
for (auto time_out : time_list) {
Real t_out = ConvertTime(t_in, time_in, time_out);
Real t_out = ConvertTime(t_diff_in, time_in, time_out);
Real t_final = ConvertTime(t_out, time_out, Time::TAI);
RequireNear(t_tai, t_final, epsilon);
}
Expand All @@ -50,10 +50,24 @@ TEST_CASE("physics.ConvertTime") {
// Orekit
auto file = OpenTestDataFile("time_converter_orekit.txt");
std::string line;
int precision = 5;
Real t_utc;
while (std::getline(file, line)) {
std::string time, date;
std::string time, date_in;
double mjd_in, t_diff_in;
std::istringstream iss(line);
iss >> time >> date;
std::cout << time << " " << date << std::endl;
iss >> time >> date_in >> mjd_in >> t_diff_in;
if (time == "UTC") t_utc = Gregorian2Time(date_in);
Real t = ConvertTime(t_utc, Time::UTC, string2time.at(time));
Real t_j2000 = ConvertTime(0, Time::TT, string2time.at(time));
Real t_diff = t - t_j2000;
Real mjd = Time2MJD(t);
std::string date = Time2GregorianString(t, precision);
RequireNear(mjd, mjd_in, epsilon);
RequireNear(t_diff, t_diff_in, epsilon);
int len = date_in.size();
if (date.size() < len) len = date.size();
len--;
REQUIRE_THAT(date.substr(0, len), Equals(date_in.substr(0, len)));
}
}
57 changes: 50 additions & 7 deletions test/data/orekit.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,17 @@
"cells": [
{
"cell_type": "code",
"execution_count": null,
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"OpenJDK 64-Bit Server VM warning: Attempt to protect stack guard pages failed.\n",
"OpenJDK 64-Bit Server VM warning: Attempt to deallocate stack guard pages failed.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
Expand All @@ -31,9 +39,37 @@
},
{
"cell_type": "code",
"execution_count": 14,
"execution_count": 113,
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"UTC 2000-05-21T00:00:00.000 51685.0 12139264.184\n",
"TAI 2000-05-21T00:00:32.000 51685.000370370224 12139264.184\n",
"GPS 2000-05-21T00:00:13.000 51685.00015046308 12139264.184\n",
"TT 2000-05-21T00:01:04.184 51685.00074287038 12139264.184\n",
"TDB 2000-05-21T00:01:04.18513759455961 51685.00074288389 12139264.185210254\n",
"TCG 2000-05-21T00:01:04.69829349107975 51685.00074882293 12139264.192460205\n",
"TCB 2000-05-21T00:01:15.62708087917957 51685.0008753133 12139264.373431945\n",
"UTC 2025-05-21T00:00:00.000 60816.0 801057664.184\n",
"TAI 2025-05-21T00:00:37.000 60816.000428240746 801057669.184\n",
"GPS 2025-05-21T00:00:18.000 60816.000208333135 801057669.184\n",
"TT 2025-05-21T00:01:09.184 60816.00080074044 801057669.184\n",
"TDB 2025-05-21T00:01:09.18514761505224 60816.00080075394 801057669.1852202\n",
"TCG 2025-05-21T00:01:10.2481136167295 60816.000813056715 801057669.7422804\n",
"TCB 2025-05-21T00:01:32.8594267225875 60816.00107476208 801057681.6057777\n",
"UTC 2050-05-21T00:00:00.000 69947.0 1589976064.184\n",
"TAI 2050-05-21T00:00:37.000 69947.00042824075 1589976069.184\n",
"GPS 2050-05-21T00:00:18.000 69947.00020833313 1589976069.184\n",
"TT 2050-05-21T00:01:09.184 69947.00080074044 1589976069.184\n",
"TDB 2050-05-21T00:01:09.18515755794682 69947.00080075394 1589976069.1852303\n",
"TCG 2050-05-21T00:01:10.79793373889461 69947.00081942044 1589976070.2921004\n",
"TCB 2050-05-21T00:01:45.09177241087139 69947.00121633988 1589976093.8381236\n"
]
}
],
"source": [
"from org.orekit.time import TimeScalesFactory, AbsoluteDate\n",
"\n",
Expand All @@ -46,14 +82,21 @@
"tcb = TimeScalesFactory.getTCB()\n",
"\n",
"times = [utc, tai, gps, tt, tdb, tcg, tcb]\n",
"prec = 8\n",
"\n",
"def round_str(txt, prec):\n",
" return str(round(int(txt[:prec+1]),-1))[:prec]\n",
"\n",
"with open(\"time_converter_orekit.txt\", \"w\") as f:\n",
" for year in [1900, 1950, 2000, 2050, 2100]:\n",
" for year in [2000, 2025, 2050]:\n",
" date = AbsoluteDate(year, 5, 21, 0, 0, 0.0, utc)\n",
" for time in times:\n",
" f.write(time.toString().ljust(3))\n",
" f.write(\" \")\n",
" f.write(f\"{date.toString(time)}\\n\")"
" line = time.toString().ljust(3) + \" \"\n",
" line += date.toString(time).ljust(35) + \" \"\n",
" line += str(date.getMJD(time)).ljust(20) + \" \"\n",
" line += str(date.offsetFrom(AbsoluteDate.J2000_EPOCH, time)).ljust(10) + \"\\n\"\n",
" f.write(line)\n",
" print(line, end=\"\")"
]
}
],
Expand Down
56 changes: 21 additions & 35 deletions test/data/time_converter_orekit.txt
View file
Original file line number Diff line number Diff line change
@@ -1,35 +1,21 @@
UTC 1900-05-21T00:00:00.000
TAI 1900-05-21T00:00:00.000
GPS 1900-05-20T23:59:41.000
TT 1900-05-21T00:00:32.184
TDB 1900-05-21T00:00:32.1851177114488
TCG 1900-05-21T00:00:30.49895276545084
TCB 1900-05-20T23:59:54.69637786926563
UTC 1950-05-21T00:00:00.000
TAI 1950-05-21T00:00:00.000
GPS 1950-05-20T23:59:41.000
TT 1950-05-21T00:00:32.184
TDB 1950-05-21T00:00:32.18513798232179
TCG 1950-05-21T00:00:31.59859300978106
TCB 1950-05-21T00:00:19.16106963091725
UTC 2000-05-21T00:00:00.000
TAI 2000-05-21T00:00:32.000
GPS 2000-05-21T00:00:13.000
TT 2000-05-21T00:01:04.184
TDB 2000-05-21T00:01:04.18513759455961
TCG 2000-05-21T00:01:04.69829349107975
TCB 2000-05-21T00:01:15.62708087917957
UTC 2050-05-21T00:00:00.000
TAI 2050-05-21T00:00:37.000
GPS 2050-05-21T00:00:18.000
TT 2050-05-21T00:01:09.184
TDB 2050-05-21T00:01:09.18515755794682
TCG 2050-05-21T00:01:10.79793373889461
TCB 2050-05-21T00:01:45.09177241087139
UTC 2100-05-21T00:00:00.000
TAI 2100-05-21T00:00:37.000
GPS 2100-05-21T00:00:18.000
TT 2100-05-21T00:01:09.184
TDB 2100-05-21T00:01:09.18517720455521
TCG 2100-05-21T00:01:11.89757398322482
TCB 2100-05-21T00:02:09.55646354825841
UTC 2000-05-21T00:00:00.000 51685.0 12139264.184
TAI 2000-05-21T00:00:32.000 51685.000370370224 12139264.184
GPS 2000-05-21T00:00:13.000 51685.00015046308 12139264.184
TT 2000-05-21T00:01:04.184 51685.00074287038 12139264.184
TDB 2000-05-21T00:01:04.18513759455961 51685.00074288389 12139264.185210254
TCG 2000-05-21T00:01:04.69829349107975 51685.00074882293 12139264.192460205
TCB 2000-05-21T00:01:15.62708087917957 51685.0008753133 12139264.373431945
UTC 2025-05-21T00:00:00.000 60816.0 801057664.184
TAI 2025-05-21T00:00:37.000 60816.000428240746 801057669.184
GPS 2025-05-21T00:00:18.000 60816.000208333135 801057669.184
TT 2025-05-21T00:01:09.184 60816.00080074044 801057669.184
TDB 2025-05-21T00:01:09.18514761505224 60816.00080075394 801057669.1852202
TCG 2025-05-21T00:01:10.2481136167295 60816.000813056715 801057669.7422804
TCB 2025-05-21T00:01:32.8594267225875 60816.00107476208 801057681.6057777
UTC 2050-05-21T00:00:00.000 69947.0 1589976064.184
TAI 2050-05-21T00:00:37.000 69947.00042824075 1589976069.184
GPS 2050-05-21T00:00:18.000 69947.00020833313 1589976069.184
TT 2050-05-21T00:01:09.184 69947.00080074044 1589976069.184
TDB 2050-05-21T00:01:09.18515755794682 69947.00080075394 1589976069.1852303
TCG 2050-05-21T00:01:10.79793373889461 69947.00081942044 1589976070.2921004
TCB 2050-05-21T00:01:45.09177241087139 69947.00121633988 1589976093.8381236

0 comments on commit d75be00

Please sign in to comment.