diff --git a/README.md b/README.md index 1aa4afe..20d2fbc 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,7 @@ # das2C - version 3.0 (in work) -*This version will support parsing das3 streams without namespaces, das3 docs will - come latter on a subsequent release* +**Note: The main branch is under development. For a stable API use the 2.3 release.** Das servers typically provide data relavent to space plasma and magnetospheric physics research. To retrieve data, an HTTP GET request is posted to a das diff --git a/das2/spice.c b/das2/spice.c index c854c82..c2876c5 100644 --- a/das2/spice.c +++ b/das2/spice.c @@ -97,53 +97,3 @@ int das_print_spice_error(const char* sProgName) return DASERR_SPICE; } - -DasErrCode das_spice_dm2et(double* pEt, const das_datum* pDatum) -{ - /* TODO: Add in hook to pull from variable if needed for sequences, etc. ! */ - - return das_error(DASERR_NOTIMP, "Time to et conversion not yet implemented"); - - /* Not yet implemented ... - das_val_type vt = DasAry_valType(pCalc->pSrc); - das_units units = DasAry_units(pCalc->pSrc); - const ubyte* pData = DasAry_getAt(pCalc->pSrc, vt, aShape); - - if(pCalc->request.uFlags & XFORM_LOC){ - - / * Time usually comes in one of three forms: - 1. double, units - 2. long, units - 3. das_time, implied UTC - - Note that spice times don't have the same precision as TT2K over the - TT2K time range. There's nothing to be done about this other then not - calculating ephemeris vectors for every single point in a 20 MHz waveform! - * / - double et; - switch(vt){ - case vtDouble: - - - } - - - - - / * Extract the time values * / - das_datum dm; - DasVar_get(pTime, [0,0,0,0]); - - / * need a fast converter from TT2K to ephemeris time * / - double et = das_tt2K_to_ephem(*(int64_t*)(&dm)); - } - -#ifdef DEBUG - / * Now check it * / - das_time dt; - dt_from_tt2k() - -#endif - - */ -} \ No newline at end of file diff --git a/das2/spice.h b/das2/spice.h index e9dcf35..0a72174 100644 --- a/das2/spice.h +++ b/das2/spice.h @@ -49,11 +49,6 @@ int das_send_spice_err(int nDasVer, const char* sErrType); */ int das_print_spice_error(const char* sProgName); -/** Convert any das time datum to a spice ephemeris time - */ - -DasErrCode das_spice_dm2et(double* pEt, const das_datum* pDatum); - /** @} */ #endif /* _das2_spice_H_ */ diff --git a/das2/units.c b/das2/units.c index 5168422..cd53294 100644 --- a/das2/units.c +++ b/das2/units.c @@ -64,6 +64,7 @@ const char* UNIT_T1970 = "t1970"; const char* UNIT_NS1970 = "ns1970"; const char* UNIT_UTC = "UTC"; const char* UNIT_TT2000 = "TT2000"; +const char* UNIT_ET2000 = "ET2000"; /* Other common units */ const char* UNIT_SECONDS = "s"; @@ -180,24 +181,25 @@ bool units_init(const char* sProgName) g_lUnits[4] = UNIT_NS1970; g_lUnits[5] = UNIT_UTC; g_lUnits[6] = UNIT_TT2000; - g_lUnits[7] = UNIT_MILLISECONDS; - g_lUnits[8] = UNIT_MICROSECONDS; - g_lUnits[9] = UNIT_NANOSECONDS; - g_lUnits[10] = UNIT_SECONDS; - g_lUnits[11] = UNIT_HOURS; - g_lUnits[12] = UNIT_DAYS; - g_lUnits[13] = UNIT_HERTZ; - g_lUnits[14] = UNIT_KILO_HERTZ; - g_lUnits[15] = UNIT_MEGA_HERTZ; - g_lUnits[16] = UNIT_E_SPECDENS; - g_lUnits[17] = UNIT_B_SPECDENS; - g_lUnits[18] = UNIT_NT; - g_lUnits[19] = UNIT_NUMBER_DENS; - g_lUnits[20] = UNIT_DB; - g_lUnits[21] = UNIT_KM; - g_lUnits[22] = UNIT_EV; - g_lUnits[23] = UNIT_DIMENSIONLESS; - g_lUnits[24] = NULL; + g_lUnits[7] = UNIT_ET2000; + g_lUnits[8] = UNIT_MILLISECONDS; + g_lUnits[9] = UNIT_MICROSECONDS; + g_lUnits[10] = UNIT_NANOSECONDS; + g_lUnits[11] = UNIT_SECONDS; + g_lUnits[12] = UNIT_HOURS; + g_lUnits[13] = UNIT_DAYS; + g_lUnits[14] = UNIT_HERTZ; + g_lUnits[15] = UNIT_KILO_HERTZ; + g_lUnits[16] = UNIT_MEGA_HERTZ; + g_lUnits[17] = UNIT_E_SPECDENS; + g_lUnits[18] = UNIT_B_SPECDENS; + g_lUnits[19] = UNIT_NT; + g_lUnits[20] = UNIT_NUMBER_DENS; + g_lUnits[21] = UNIT_DB; + g_lUnits[22] = UNIT_KM; + g_lUnits[23] = UNIT_EV; + g_lUnits[24] = UNIT_DIMENSIONLESS; + g_lUnits[25] = NULL; return true; } @@ -246,6 +248,7 @@ das_units Units_interval(das_units unit){ if(unit == UNIT_NS1970) return UNIT_NANOSECONDS; if(unit == UNIT_UTC) return UNIT_SECONDS; if(unit == UNIT_TT2000) return UNIT_NANOSECONDS; + if(unit == UNIT_ET2000) return UNIT_SECONDS; return unit; } @@ -255,7 +258,7 @@ bool Units_isInterval(das_units unit){ return ( (unit == UNIT_US2000) || (unit == UNIT_MJ1958) || (unit == UNIT_T2000) || (unit == UNIT_T1970) || (unit == UNIT_NS1970) || (unit == UNIT_UTC) || - (unit == UNIT_TT2000) + (unit == UNIT_TT2000) || (unit == UNIT_ET2000) ); } @@ -265,7 +268,7 @@ bool Units_haveCalRep(das_units unit){ return (unit == UNIT_US2000) || (unit == UNIT_MJ1958) || (unit == UNIT_T2000) || (unit == UNIT_T1970) || (unit == UNIT_NS1970) || (unit == UNIT_UTC) || - (unit == UNIT_TT2000) ; + (unit == UNIT_TT2000) || (unit == UNIT_ET2000); } /* ************************************************************************* */ @@ -995,7 +998,7 @@ char* Units_toLabel(das_units unit, char* sBuf, int nLen) if(unit == NULL){ memset(sBuf, 0, nLen); return sBuf; } if((unit == UNIT_US2000) || (unit == UNIT_MJ1958) || (unit == UNIT_T2000) || (unit == UNIT_T1970) || (unit == UNIT_NS1970) || (unit == UNIT_UTC) - || (unit == UNIT_TT2000) ){ + || (unit == UNIT_TT2000)|| (unit == UNIT_ET2000) ){ snprintf(sBuf, nLen - 1, "UTC"); return sBuf; } @@ -1505,6 +1508,41 @@ bool Units_canConvert(das_units from, das_units to ) return true; } +/* Thanks SPICE leapseconds.tls kernel! -cwp */ +#define K 1.657e-3 +#define EB 1.671e-2 +#define M0 6.239996 +#define M1 1.99096871e-7 + +double Units_et2k_to_tt2k(double rET) +{ + double M = M0 + M1*rET; + double E = M + EB * sin(M); + double rTT = rET - K * sin( E ); + return rTT*1e9; +} + +/* TODO: + * Calculations below depend on sin(TT*1e-9) and sin(ET) being almost + * identical, which they are, however this conversion could be made to + * exactly match SPICE to 16 decimal places. To do so, change the TT + * to ET constants ever so slightly. + */ + +double Units_tt2k_to_et2k(double rTT) +{ + double M = M0 + M1*(rTT*1e-9); /* rTT ~ rET Assumption not full proof near J000 epoch */ + double E = M + EB * sin(M); + double rET = rTT*1e-9 + K * sin( E ); + return rET; +} + +#undef K +#undef EB +#undef M0 +#undef M1 + + /* Special conversion helpers, needs an upgrade when we switch to datums */ double _Units_convertToUS2000( double value, das_units fromUnits ) { @@ -1516,6 +1554,10 @@ double _Units_convertToUS2000( double value, das_units fromUnits ) if(fromUnits == UNIT_T1970) return (value - 946684800 ) * 1e6; if(fromUnits == UNIT_NS1970) return (value - 9.466848e+17 ) * 1e-3; if(fromUnits == UNIT_TT2000) return das_tt2K_to_us2K(value); + if(fromUnits == UNIT_ET2000){ + value = Units_et2k_to_tt2k(value); + return das_tt2K_to_us2K(value); + } das_error(DASERR_UNITS, "unsupported conversion to US2000 from %s\n", Units_toStr(fromUnits) ); @@ -1532,6 +1574,10 @@ double _Units_convertFromUS2000( double value, das_units toUnits ) { if(toUnits == UNIT_T1970) return value / 1e6 + 946684800; if(toUnits == UNIT_NS1970) return (value + 9.466848e+14) * 1e3; if(toUnits == UNIT_TT2000) return das_us2K_to_tt2K(value); + if(toUnits == UNIT_ET2000){ + value = das_us2K_to_tt2K(value); + return Units_tt2k_to_et2k(value); + } das_error(DASERR_UNITS, "unsupported conversion from US2000 to %s\n", Units_toStr(toUnits) @@ -1554,6 +1600,13 @@ double Units_convertTo(das_units to, double rFrom, das_units from) double rUs2k = 0.0, rTo = 0.0; if(Units_haveCalRep(to) && Units_haveCalRep(from)){ + /* TT2000 & ET2000 have no time loss across leap seconds, don't + go out to a lossy time when converting between the two */ + if((to == UNIT_ET2000)&&(from == UNIT_TT2000)) + return Units_tt2k_to_et2k(rFrom); + if((to == UNIT_TT2000)&&(from == UNIT_ET2000)) + return Units_et2k_to_tt2k(rFrom); + rUs2k = _Units_convertToUS2000( rFrom, from ); rTo = _Units_convertFromUS2000( rUs2k, to ); @@ -1736,7 +1789,6 @@ static int days[2][14] = { #define LEAP(y) ((y) % 4 ? 0 : ((y) % 100 ? 1 : ((y) % 400 ? 0 : 1))) - DasErrCode Units_convertToDt(das_time* pDt, double value, das_units epoch_units) { dt_null(pDt); @@ -1745,8 +1797,15 @@ DasErrCode Units_convertToDt(das_time* pDt, double value, das_units epoch_units) #define DAY_1958_TO_1970 (12*365 + 3) /* 12 years and 3 leap days */ double rUnix; + /* If this is SPICE ephmeris time we'll need to convert to TT2k first */ + if(epoch_units == UNIT_ET2000){ + value = Units_et2k_to_tt2k(value); + epoch_units = UNIT_TT2000; + } + /* If input is TT2K, use a dedicated converter that allows seconds field > 59 */ if(epoch_units == UNIT_TT2000){ + long long ntt2k = (long long) value; double yr, mo, dy, hr, mi, sc, ms, us, ns; @@ -1834,7 +1893,7 @@ DasErrCode Units_convertToDt(das_time* pDt, double value, das_units epoch_units) double Units_convertFromDt(das_units epoch_units, const das_time* pDt) { /* If input is TT2K, use a dedicated converter that allows seconds field > 59 */ - if(epoch_units == UNIT_TT2000){ + if((epoch_units == UNIT_TT2000)||(epoch_units == UNIT_ET2000)){ double sc = (int)pDt->second; double ms = (int)( (pDt->second - sc)*1e3 ); double us = (int)( (pDt->second - sc - ms*1e-3)*1e6 ); @@ -1846,7 +1905,14 @@ double Units_convertFromDt(das_units epoch_units, const das_time* pDt) double hr = pDt->hour; double mn = pDt->minute; long long ntt2k = das_utc_to_tt2K(yr, mt, dy, hr, mn, sc, ms, us, ns); - return (double)ntt2k; + + double rTT2k = (double)ntt2k; + + /* Out convert from TT2K if ET2000 is the actual desire */ + if(epoch_units == UNIT_ET2000) + return Units_tt2k_to_et2k(rTT2k); + else + return rTT2k; } double mj1958 = 0.0; diff --git a/das2/units.h b/das2/units.h index a237800..c655742 100644 --- a/das2/units.h +++ b/das2/units.h @@ -48,6 +48,11 @@ extern const char* UNIT_UTC; /** Time strings on the Gregorian Calendar */ last know leapsecond in the data time. */ extern const char* UNIT_TT2000; +/** SPICE ephem time, constant seconds for Solar Sys Center in J2000 epoch + * This is the same UTC epoch as TT2000, but does not move at 1 second per second + * from the Terrestrial point of view. */ +extern const char* UNIT_ET2000; + /* Other common units */ extern const char* UNIT_SECONDS; /** Units of SI Seconds */ extern const char* UNIT_HOURS; /** Units of 3600 seconds */ diff --git a/test/FakeLeapSeconds.txt b/test/FakeLeapSeconds.txt index cf777c0..936b9f5 100644 --- a/test/FakeLeapSeconds.txt +++ b/test/FakeLeapSeconds.txt @@ -1,7 +1,7 @@ ; FAKE! Leap Seconds Table - in the style used by CDF, for testing ONLY! ; ; Do NOT use this file in production code! Contains artificial -; leap socond at Jan. 1st 2021. +; leap second at Jan. 1st 2021. ; ; Comment lines starts with ";" at column 1. ; Year Month Day Leap Seconds Drift diff --git a/test/TestSpice.c b/test/TestSpice.c index 3e2008c..9501f20 100644 --- a/test/TestSpice.c +++ b/test/TestSpice.c @@ -3,9 +3,44 @@ #define _POSIX_C_SOURCE 200112L #include +#include + #include #include +#include + +void prnConvert(double rET){ + char sBuf[32] = {'\0'}; + et2utc_c(rET, "ISOC", 12, 31, sBuf); + das_time dt; + dt_parsetime(sBuf, &dt); + + int64_t nTT = dt_to_tt2k(&dt); + printf("ET %20.9f is %s UTC, %17ld TT (TT-ET is %f)\n", rET, sBuf, nTT, nTT*1e-9 - rET); + + /* TT and ET have the same epoch, so once you get to TT it's just: + ET = TT*1e-9 + K sin( E ) + K = 1.657e-3 + E = M + EB sin M + M = M0 + M1*ET + EB = 1.671e-2 + M0 = 6.239996 + M1 = 1.99096871e-7 + from there. + */ + + const double K = 1.657e-3 ; + const double EB = 1.671e-2; + const double M0 = 6.239996; + const double M1 = 1.99096871e-7; + double M = M0 + M1*(nTT*1e-9); + double E = M + EB * sin(M); + double rMyET = nTT*1e-9 + K * sin( E ); + + printf("myET %20.9f (delta Mine - Spice) %.12f ms\n", rMyET, (rMyET - rET)*1e6); +} + int main(int argc, char** argv) { @@ -13,8 +48,31 @@ int main(int argc, char** argv) das_init(argv[0], DASERR_DIS_EXIT, 0, DASLOG_INFO, NULL); das_spice_err_setup(); /* Make sure spice errors don't go to stdout */ - + printf("INFO: Can redirect spice errors from stdout, good\n"); - + + /* Load the leap seconds kernel */ + furnsh_c("test/leapseconds.tls"); + + /* Print ET 0 */ + prnConvert(-86400*366*20); + prnConvert(-86400*366); + prnConvert(-86400*274.5); + prnConvert(-86400*183); + prnConvert(-86400*91.5); + prnConvert(-86400*10.0); + prnConvert(-86400*1.0); + + prnConvert(86400*0.0); + + prnConvert(86400*1.0); + prnConvert(86400*10.0); + prnConvert(86400*91.5); + prnConvert(86400*183); /* Mean eccentricity should be opposite here */ + prnConvert(86400*274.5); + prnConvert(86400*366); + + prnConvert(86400*366*20); + return 0; } diff --git a/test/TestTT2000.c b/test/TestTT2000.c index ba952be..98d279a 100644 --- a/test/TestTT2000.c +++ b/test/TestTT2000.c @@ -15,6 +15,7 @@ #define _POSIX_C_SOURCE 200112L #include +#include #include @@ -96,6 +97,71 @@ struct map_t { double ms, us, ns; }; +/* Test conversions between TT2000 and ET2000. The difference should + be on the order of the cumulative sesonal general relativity effect + on time at Earth (~0.001658 seconds) +*/ + +bool testETtoTT(double rEt2k, int nTest) +{ + double rTt2k = Units_convertTo(UNIT_TT2000, rEt2k, UNIT_ET2000); + double rDiff = (rEt2k*1e3) - (rTt2k*1e-6); + double rExpectMsDiff = 1.658; + if(fabs(rDiff) > (rExpectMsDiff*1.1)){ + printf( + "ERROR: Test %d failed, magnitude of the difference between TT2000 " + "and SPICE ephemeris time is %.3f ms expected %.3f ms or less.\n", + nTest, rDiff, rExpectMsDiff + ); + return false; + } + /* For the times that have been selected there should be a measurable + difference. Two times a year they match exactly if the feeder + arrays change, change this code too */ + /* printf("Test %d: TT2000, ET2000 difference: %.3g milliseconds\n", nTest, rDiff); */ + if(fabs(rDiff) < 0.001){ + printf( + "ERROR: Test %d failed, expected some difference between TT2000 " + "and ET2000, found %.3f ms. This should only happen twice a year. " + "Did the unittest input arrays change?\n", nTest, fabs(rDiff) + ); + } + + /* For good measure test a trip out to us2000 and back */ + double rUs2k = Units_convertTo(UNIT_US2000, rEt2k, UNIT_ET2000); + double rEtRev = Units_convertTo(UNIT_ET2000, rUs2k, UNIT_US2000); + + das_time dt; + Units_convertToDt(&dt, rEt2k, UNIT_ET2000); + char sBuf[32]; + dt_isoc(sBuf, 31, &dt, 9); + if(fabs(rEt2k - rEtRev) > 10.0){ + printf( + "ERROR: Test %d, UTC %s round trip from ET2000 to us2000 diff: %.3f µs > 10 µs\n", + nTest, sBuf, fabs(rEt2k - rEtRev) + ); + return false; + } + + return true; +} + +bool testTTtoET(double rTt2k, int nTest) +{ + double rEt2k = Units_convertTo(UNIT_ET2000, rTt2k, UNIT_TT2000); + double rDiff = (rTt2k*1e-6) - (rEt2k*1e3) ; + double rExpectMsDiff = 1.658; + if(fabs(rDiff) > (rExpectMsDiff*1.1)){ + printf( + "ERROR: Test %d failed, magnitude of the difference between TT2000 " + "and SPICE ephemeris time is %.3f ms expected %.3f ms or less.\n", + nTest, rDiff, rExpectMsDiff + ); + return false; + } + return true; +} + int main(int argc, char** argv) { int i; /* always need one of these */ @@ -299,7 +365,35 @@ int main(int argc, char** argv) printf("ERROR: Test 8 failed, time calculations not altered by external table\n"); return 13; } - + + /* Test conversions between the two leap-second aware time scales */ + int nTest = 9; + double aSolarTime[] = { + 1262304000.0, /* sometime around July 2040 */ + 946728000.0, /* sometime around July 2030 */ + 646704000.0, /* sometime around July 2020 */ + 315619200.0, /* sometime around January 1st, 2010 */ + 157852800.0, /* sometime around January 1st, 2005 */ + 0.0, /* At the J2000 epoch */ + -632448000.0 /* sometime around July 1980 */ + }; + double aEarthTime[] = { + 1.262304e+18, /* sometime around January 2040 */ + 9.46728e+17, /* sometime around July 2030 */ + 6.46704e+17, /* sometime around July 2020 */ + 3.156192e+17, /* sometime around January 1st, 2010 */ + 1.578528e+17, /* sometime around January 1st, 2005 */ + 0.0, /* At the J2000 epoch */ + -6.31152e+17 /* sometime around July 1980 */ + }; + for(int i = 0; i < 7; ++i){ + if(!testETtoTT(aSolarTime[i], nTest)) + return nTest; + if(!testTTtoET(aEarthTime[i], nTest)) + return nTest; + ++nTest; + } + printf("INFO: All TT2000 tests passed\n"); return 0; } diff --git a/utilities/das3_spice.c b/utilities/das3_spice.c index cf5e402..631c1fc 100644 --- a/utilities/das3_spice.c +++ b/utilities/das3_spice.c @@ -1059,6 +1059,14 @@ DasErrCode onDataSet(DasStream* pSdIn, int iPktId, DasDs* pDsIn, void* pUser) /* ************************************************************************ */ /* Data output */ +DasErrCode _dm2et(double* pOutput, const das_datum* pInput) +{ + if(pInput->vt == vtTime) + return Units_convertFromDt(UNIT_ET2000, (const das_time*)pInput); + + return Units_convertTo(UNIT_ET2000, das_datum_toDbl(pInput), pInput->units); +} + /* We want the input dataset here because we need to know how many provided values we're going to get */ DasErrCode _writeLocation(DasDs* pDsIn, XCalc* pCalc) @@ -1090,7 +1098,8 @@ DasErrCode _writeLocation(DasDs* pDsIn, XCalc* pCalc) for(; !iter.done; das_uniq_iter_next(&iter)){ DasVar_get(pCalc->pTime, iter.index, &dm); - das_spice_dm2et(&rEt, &dm); + if(!_dm2et(&rEt, &dm)) + return das_error(PERR, "Conversion of datum to ephemeris time failed"); spkezp_c( pReq->nBodyId, rEt, pReq->aOutFrame, "NONE", pReq->nOutCenter, aRecOut, &rLt @@ -1187,7 +1196,10 @@ DasErrCode _writeRotation(DasDs* pDsIn, XCalc* pCalc) for(; !iter.done; das_uniq_iter_next(&iter)){ DasVar_get(pCalc->pTime, iter.index, &dm); - das_spice_dm2et(&rEt, &dm); + + if(!_dm2et(&rEt, &dm)) + return das_error(PERR, "Conversion of datum to ephemeris time failed"); + pxform_c(pReq->aInFrame, pReq->aOutFrame, rEt, mRot); /* Get rot matrix */ DasVar_get(pVarIn, iter.index, &dm);