Added SPICE Ephemeris as an intrinsic time type

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 

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	
-
-   */
-}

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_ */

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;

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 */

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

test/TestSpice.c

@@ -3,18 +3,76 @@
 #define _POSIX_C_SOURCE 200112L
 
 #include <stdio.h>
+#include <math.h>
+
 #include <das2/core.h>
 #include <das2/spice.h>
 
+#include <SpiceUsr.h>
+
+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)
 {
 
 	/* Exit on errors, log info messages and above */
 	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;
 }