Index: trunk/psLib/src/astronomy/psTime.c =================================================================== --- trunk/psLib/src/astronomy/psTime.c (revision 4029) +++ trunk/psLib/src/astronomy/psTime.c (revision 4051) @@ -10,6 +10,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.61 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-05-25 20:26:55 $ + * @version $Revision: 1.62 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-05-31 21:47:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -24,4 +24,5 @@ #include "psTime.h" #include "psError.h" +#include "psLogMsg.h" #include "psMemory.h" #include "psAbort.h" @@ -69,9 +70,30 @@ static psMetadata *timeMetadata = NULL; +// Offset to convert terrestrial time(TT) to international atomic time(TAI) +#define TAI_TT_OFFSET_SECONDS 32 +#define TAI_TT_OFFSET_NANOSECONDS 184000000 + +// Offset from converting to MJD +#define MJD_EPOCH_OFFSET 40587.0 + +// Offset for converting to JD +#define JD_EPOCH_OFFSET 2440587.5 + +// Offset of year 0000 from epoch +#define YEAR_0000_SEC -62125920000.0 + +// Offset of year 9999 from epoch +#define YEAR_9999_SEC 253202544000.0 + /** Static function prototypes */ static char *cleanString(char *inString, int sLen); static char* getToken(char **inString, char *delimiter, psParseErrorType *status); -psF64 searchTables(psF64 index, psU64 column, psLookupStatusType *status, char *metadataTableNames[], psU32 nTables); - +static psF64 searchTables(psF64 index, psU64 column, char *metadataTableNames[], + psU32 nTables, psLookupStatusType* status); +static psTime* convertTimeTAIUTC(psTime* time); +static psTime* convertTimeUTCTAI(psTime* time); +static psTime* convertTimeTAITT(psTime* time); +static psTime* convertTimeTTTAI(psTime* time); +static psTime* convertTimeUTCUT1(psTime* time); /** Removes leading and trailing whitespace and # characters from a string. The cleaned string is a new null @@ -83,14 +105,13 @@ char *cleaned = NULL; - ptrB = inString; - /* Skip over leading # or whitespace */ + // Skip over leading # or whitespace while (isspace(*ptrB) || *ptrB=='#') { ptrB++; } - /* Skip over trailing whitespace, null terminators, and # characters */ - ptrE = inString + sLen; + // Skip over trailing whitespace, null terminators, and # characters + ptrE = inString + sLen - 1; while(isspace(*ptrE) || *ptrE=='\0' || *ptrE=='#') { ptrE--; @@ -113,5 +134,4 @@ int sLen = 0; - // Skip over leading whitespace while(isspace(**inString)) { @@ -121,4 +141,5 @@ // Length of token, not including delimiter sLen = strcspn(*inString, delimiter); + if(sLen) { @@ -127,5 +148,6 @@ // Move to end of token - (*inString) += sLen; + (*inString) += (sLen+1); + } else if(**inString!='\0' && sLen==0) { *status = PS_PARSE_ERROR_GENERAL; @@ -148,45 +170,58 @@ -/** Searches time tables in priority order and performs interpolation if input index value is within a table. - * If the index value is out of range, the status is set accordingly. */ -psF64 searchTables(psF64 index, psU64 column, psLookupStatusType *status, char *metadataTableNames[], psU32 nTables) -{ - psU32 i = 0; - char *tableName = NULL; - psF64 result = 0.0; - psLookupTable *table = NULL; - psMetadataItem *tableMetadataItem = NULL; - - +// Searches time tables in priority order and performs interpolation if input index value is within a table. +// If the index value is out of range, the status is set accordingly. +static psF64 searchTables(psF64 index, psU64 column, char *metadataTableNames[], + psU32 nTables, psLookupStatusType* status) +{ + char* tableName = NULL; + psF64 result = NAN; + psLookupTable* table = NULL; + psMetadataItem* tableMetadataItem = NULL; // Check time metadata. Function call reports errors. if(timeMetadata == NULL) { if(!p_psTimeInit(p_psGetConfigFileName())) - return 0.0; - } - - // Search each table in priority order: ser7, eopc,finals - for(i=0; idata.V; - PS_ASSERT_PTR_NON_NULL(table,0.0); - - // Attempt to interpolate table - *status = PS_LOOKUP_SUCCESS; - result = psLookupTableInterpolate(table, index, 2, status); - if(*status == PS_LOOKUP_SUCCESS) { - return result; - } else if(*status == PS_LOOKUP_ERROR) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED_NAME, - tableName); - return 0.0; + + // Check that table is not NULL + PS_ASSERT_PTR_NON_NULL(table,NAN); + + // Check if index within to/from range + if(index >= table->validFrom ) { + if(index <= table->validTo) { + // Attempt to interpolate table + result = psLookupTableInterpolate(table, index, column); + *status = PS_LOOKUP_SUCCESS; + if(!isnan(result)) { + break; + } + } else { + *status = PS_LOOKUP_PAST_BOTTOM; + } + } else { + *status = PS_LOOKUP_PAST_TOP; } } @@ -197,10 +232,14 @@ bool p_psTimeInit(const char *fileName) { + psS32 numLines = 0; bool foundTable = false; char *tableDir = NULL; char *tableNames = NULL; + char *tableFormats = NULL; char *namesPtr = NULL; + char *formatPtr = NULL; char *metadataNamesPtr = NULL; char *tableName = NULL; + char *tableFormat = NULL; char *fullTableName = NULL; psS32 i = 0; @@ -210,8 +249,9 @@ psVector *tablesFrom = NULL; psVector *tablesTo = NULL; + psVector *tablesIndex = NULL; psMetadataItem *metadataItem = NULL; psLookupTable *table = NULL; psParseErrorType status = PS_PARSE_SUCCESS; - char metadataTableNames[4][MAX_STRING_LENGTH] = {"ser7", "eopc", "finals", "tai"}; + char metadataTableNames[4][MAX_STRING_LENGTH] = {"daily", "eopc", "finals", "tai"}; // Read time config file @@ -251,4 +291,5 @@ psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.tables.to"); + psFree(tablesFrom); return false; } @@ -258,4 +299,22 @@ psFree(tablesFrom); psFree(tablesTo); + return false; + } + + // Get index columns for the tables + metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.index"); + if(metadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.tables.index"); + psFree(tablesFrom); + psFree(tablesTo); + return false; + } + tablesIndex = psVectorCopy(tablesIndex, metadataItem->data.V, PS_TYPE_S32); + if(tablesIndex->n != numTables) { + psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_BAD_VECTOR,tablesIndex->n,numTables); + psFree(tablesFrom); + psFree(tablesTo); + psFree(tablesIndex); return false; } @@ -268,4 +327,5 @@ psFree(tablesFrom); psFree(tablesTo); + psFree(tablesIndex); return false; } @@ -280,8 +340,24 @@ psFree(tablesFrom); psFree(tablesTo); + psFree(tablesIndex); psFree(tableDir); return false; } tableNames = psStringCopy(metadataItem->data.V); + + // Get table format strings + metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.format"); + if(metadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE,true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.tables.format"); + psFree(tablesFrom); + psFree(tablesTo); + psFree(tablesIndex); + psFree(tableDir); + psFree(tableNames); + return false; + } + tableFormats = psStringCopy(metadataItem->data.V); + formatPtr = tableFormats; // Read time tables @@ -289,5 +365,6 @@ while((tableName=getToken(&namesPtr, " ", &status)) != NULL) { - // Form path with table name, adding one to length for last '/' that may not occur in string in cong file + // Form path with table name, adding one to length for last '/' that may not occur + // in string in cong file fullTableName = (char*)psAlloc(strlen(tableDir)+strlen(tableName)+1+1); @@ -298,8 +375,15 @@ strcat(fullTableName, tableName); + // Get table format + tableFormat = getToken(&formatPtr,",",&status); + if(tableFormat == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.tables.format"); + } + // Create and read table if(i < numTables) { - table = psLookupTableAlloc(fullTableName, tablesFrom->data.F64[i], tablesTo->data.F64[i]); - table = psLookupTableRead(table); + table = psLookupTableAlloc(fullTableName, (const char*)tableFormat, tablesIndex->data.S32[i]); + numLines = psLookupTableRead(table); } else { psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_TABLE_COUNT, i+1, numTables); @@ -321,4 +405,5 @@ psFree(fullTableName); psFree(tableName); + psFree(tableFormat); psFree(table); i++; @@ -333,4 +418,6 @@ psFree(tablesFrom); psFree(tablesTo); + psFree(tablesIndex); + psFree(tableFormats); return true; @@ -351,7 +438,7 @@ psTime *outTime = NULL; - // Error checks - if(type!=PS_TIME_TAI && type!=PS_TIME_UTC) { + if(type!=PS_TIME_TAI && type!=PS_TIME_UTC && type!=PS_TIME_UT1 && + type!=PS_TIME_TT) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, @@ -360,6 +447,8 @@ } + // Allocate memory for structure outTime = (psTime*)psAlloc(sizeof(psTime)); + // Initialize members outTime->sec = 0; outTime->nsec = 0; @@ -375,8 +464,13 @@ psTime *time = NULL; - // Allocate psTime struct time = psTimeAlloc(type); + // Verify time structure allocated + if(time == NULL) { + return NULL; + } + + // Get the system time if (gettimeofday(&now, (struct timezone *)0) == -1) { psError(PS_ERR_OS_CALL_FAILED, true, @@ -391,5 +485,5 @@ // Add most leapseconds to UTC time to get TAI time if necessary if(type == PS_TIME_TAI) { - time->sec += psTimeGetTAIDelta(time); + time->sec += p_psTimeGetTAIDelta(time); } @@ -397,25 +491,239 @@ } +static psTime* convertTimeTAIUTC(psTime* time) +{ + psF64 deltaTAI = 0.0; + psS64 deltaSec = 0; + psU32 deltaNsec = 0; + psF64 deltaUTC = 0.0; + + // Determine delta to convert between UTC and TAI + deltaTAI = p_psTimeGetTAIDelta(time); + deltaSec = (psS64)(deltaTAI); + deltaNsec = (psU32)((deltaTAI - (psF64)deltaSec) * 1e9); + + // Determine seconds + time->sec -= deltaSec; + + // Check for underflow in nsec + if(deltaNsec > time->nsec) { + // Borrow second + time->nsec += 1e9; + time->sec--; + } + + // Determine nsec + time->nsec -= deltaNsec; + + // Check for overflow in nsec + if(time->nsec >= 1e9) { + time->nsec -= 1e9; + time->sec++; + } + + // Set new type + time->type = PS_TIME_UTC; + + // Check if leapsecond present in delta + deltaUTC = p_psTimeGetTAIDelta(time); + if(fabs(deltaTAI-deltaUTC) >= 1.0) { + time->sec++; + } + + return time; +} + +static psTime* convertTimeUTCTAI(psTime* time) +{ + psF64 delta = 0.0; + psS64 deltaSec = 0; + psU32 deltaNsec = 0; + + // Determine delta to convert between UTC and TAI + delta = p_psTimeGetTAIDelta(time); + + deltaSec = (psS64)(delta); + deltaNsec = (psU32)((delta - (psF64)deltaSec) * 1e9); + + // Determine seconds + time->sec += deltaSec; + + // Determine nsec + time->nsec += deltaNsec; + + // Check for overflow in nsec + if(time->nsec >= 1e9) { + time->nsec -= 1e9; + time->sec++; + } + + // Set new type + time->type = PS_TIME_TAI; + + return time; +} + +static psTime* convertTimeTAITT(psTime* time) +{ + // Add TT offset + time->sec += TAI_TT_OFFSET_SECONDS; + time->nsec += TAI_TT_OFFSET_NANOSECONDS; + + // Check for overflow in nsec + if(time->nsec >= 1e9) { + time->nsec -= 1e9; + time->sec++; + } + + // Set new type + time->type = PS_TIME_TT; + + return time; +} + +static psTime* convertTimeTTTAI(psTime* time) +{ + // Subtract TT offset + time->sec -= TAI_TT_OFFSET_SECONDS; + + // Check for nsec underflow + if(TAI_TT_OFFSET_NANOSECONDS > time->nsec) { + // Borrow second + time->sec--; + time->nsec += 1e9; + } + time->nsec -= TAI_TT_OFFSET_NANOSECONDS; + + // Check for overflow in nsec + if(time->nsec >= 1e9) { + time->nsec -= 1e9; + time->sec++; + } + + // Set new type + time->type = PS_TIME_TAI; + + return time; +} + +static psTime* convertTimeUTCUT1(psTime* time) +{ + psS64 ut1utc = 0; + + // Get UT1-UTC value + ut1utc = (psS64)(psTimeGetUT1Delta(time,PS_IERS_A) * 1e9); + + // Since UTC is within 0.9 sec of UT1 then nsec member is the member affected + if((ut1utc < 0) && (abs(ut1utc) > time->nsec)) { + // Borrow from sec + time->sec--; + if(time->leapsecond) { + time->leapsecond = false; + } else { + time->leapsecond = psTimeIsLeapSecond(time); + } + // Add to nsec + time->nsec += 1e9; + } + time->nsec += ut1utc; + + // Check for overflow in nsec + if(time->nsec >= 1e9) { + time->nsec -= 1e9; + time->sec++; + if(time->leapsecond) { + time->leapsecond = false; + time->sec--; + } else { + time->leapsecond = psTimeIsLeapSecond(time); + } + } + + // Set new type + time->type = PS_TIME_UT1; + + return time; +} + psTime* psTimeConvert(psTime *time, psTimeType type) { - double delta = 0.0; - - // Error checks PS_ASSERT_PTR_NON_NULL(time,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,1e9-1,NULL); - - if (time->type == type) { // time already right type. That was easy! + PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),time); + + // If the input type is UT1 then return time and generate error message + if(time->type == PS_TIME_UT1) { + psError(PS_ERR_BAD_PARAMETER_VALUE,true,"Cannot convert from UT1 time type"); return time; } - delta = psTimeGetTAIDelta(time); - - if (type == PS_TIME_UTC) { - time->sec = time->sec - delta; // User wants UTC time. Subtract leapseconds. - } else if(type == PS_TIME_TAI) { - time->sec = time->sec + delta; // User wants TAI time. Add leapseconds. + // If the time to convert to is the same as psTime the return time + if (time->type == type) { + return time; + } + + // Convert from TAI to UTC, TT, UT1 + if(time->type == PS_TIME_TAI) { + // Convert from TAI to UTC + if(type == PS_TIME_UTC) { + time = convertTimeTAIUTC(time); + // Convert from TAI to TT + } else if(type == PS_TIME_TT) { + time = convertTimeTAITT(time); + // Convert from TAI to UT1 + } else if(type == PS_TIME_UT1) { + // Convert to UTC first + time = convertTimeTAIUTC(time); + // Convert UTC to UT1 + time = convertTimeUTCUT1(time); + // Convert from TAI to unknown time type + } else { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, type); + } + // Convert from TT to TAI, UTC, UT1 + } else if(time->type == PS_TIME_TT) { + // Convert from TT to UTC + if(type == PS_TIME_UTC) { + // Convert to TAI time first + time = convertTimeTTTAI(time); + // Convert from TAI to UTC + time = convertTimeTAIUTC(time); + // Convert from TT to TAI + } else if(type == PS_TIME_TAI) { + time = convertTimeTTTAI(time); + // Convert from TT to UT1 + } else if(type == PS_TIME_UT1) { + // Convert to UTC first + // Convert to TAI time first + time = convertTimeTTTAI(time); + // Convert from TAI to UTC + time = convertTimeTAIUTC(time); + // Convert from UTC to UT1 + time = convertTimeUTCUT1(time); + // Convert from TT to unknown time type + } else { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, type); + } + // Convert from UTC to TAI, TT, UT1 + } else if(time->type == PS_TIME_UTC) { + // Convert UTC to TAI + if(type == PS_TIME_TAI) { + time = convertTimeUTCTAI(time); + // Convert UTC to TT + } else if(type == PS_TIME_TT) { + // Convert to TAI time first + time = convertTimeUTCTAI(time); + // Convert TAI to TT + time = convertTimeTAITT(time); + // Convert UTC to UT1 + } else if(type == PS_TIME_UT1) { + time = convertTimeUTCUT1(time); + // Convert UTC to unknown time type + } else { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, type); + } + // Convert unknown time type } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, type); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, time->type); } @@ -423,5 +731,5 @@ } -double psTimeToLMST(psTime *time, double longitude) +psF64 psTimeToLMST(psTime *time, psF64 longitude) { psF64 jdTdtDays = 0.0; @@ -433,5 +741,4 @@ psF64 t = 0.0; psF64 tu = 0.0; - psF64 ut1UtcDbl = 0.0; psF64 const1 = 24110.5493771; psF64 const2 = 8639877.3173760; @@ -440,34 +747,30 @@ psF64 const5 = -0.0000062; psF64 const6 = 0.0000013; - psTime *ut1UtcDelta = NULL; psTime *tdtTime = NULL; - psTime *taiTime = NULL; - psTime *utcTime = NULL; psTime *ut1Time = NULL; - // Error checks PS_ASSERT_PTR_NON_NULL(time,NAN); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,1e9-1,NAN); - - // Calculate TAI or UTC time based on type of time user passes - if(time->type == PS_TIME_TAI) { - taiTime = psMemIncrRefCounter(time); - utcTime = psTimeAlloc(PS_TIME_UTC); - utcTime->sec = taiTime->sec - psTimeGetTAIDelta(time); - utcTime->nsec = taiTime->nsec; - } else if(time->type == PS_TIME_UTC) { - utcTime = psMemIncrRefCounter(time); - taiTime = psTimeAlloc(PS_TIME_TAI); - taiTime->sec = utcTime->sec + psTimeGetTAIDelta(time); - taiTime->nsec = utcTime->nsec; - } - - // Convert Universal Time (UTC) to UT1 - ut1UtcDbl = psTimeGetUT1Delta(taiTime); - utcTime->type = PS_TIME_TAI; // Don't allow addition of leapseconds to ut1Time - ut1Time = psTimeMath(utcTime, ut1UtcDbl); - utcTime->type = PS_TIME_UTC; - + PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN); + + // Verify input time is not in UT1 seconds + if(time->type == PS_TIME_UT1) { + psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_TYPE_INCORRECT,time->type); + return NAN; + } + + // Determine time reference to TT + tdtTime = psTimeAlloc(time->type); + tdtTime->sec = time->sec; + tdtTime->nsec = time->nsec; + tdtTime->leapsecond = time->leapsecond; + tdtTime = psTimeConvert(tdtTime,PS_TIME_TT); + + // Determine time reference to UT1 + ut1Time = psTimeAlloc(time->type); + ut1Time->sec = time->sec; + ut1Time->nsec = time->nsec; + ut1Time->leapsecond = time->leapsecond; + ut1Time = psTimeConvert(ut1Time,PS_TIME_UT1); // Calculate UT1 as Julian Centuries since J2000.0 @@ -476,7 +779,4 @@ t = (jdUt1Days - 2451545.0)/36525.0; - // Calculate Terrestial Dynamical Time (TDT) - tdtTime = psTimeMath(taiTime, 32.184); - // Calculate TDT as Julian centuries since J2000.0 jdTdtDays = psTimeToJD(tdtTime); @@ -486,6 +786,8 @@ fracDays = fmod(mjdUt1Days, 1.0); - // Calculate Greenwich Mean Sidereal Time (GMST) in radians. Equation set up to minimize multiplications. - gmstRad = fracDays*TWOPI+(const1+const2*tu+t*(const3+t*(const4+t*(const5+const6*t))))*S2R; + // Calculate Greenwich Mean Sidereal Time (GMST) in radians. + // Equation set up to minimize multiplications. + gmstRad = fracDays*TWOPI + + (const1+const2*tu+t*(const3+t*(const4+t*(const5+const6*t))))*S2R; // Place GMST between 0 and 2*pi @@ -496,38 +798,43 @@ // Free temporary structs - psFree(ut1UtcDelta); psFree(ut1Time); psFree(tdtTime); - psFree(utcTime); - psFree(taiTime); return lmstRad; } -double psTimeGetUT1Delta(const psTime *time) -{ - psU32 nTables = 3; - psF64 mjd = 0.0; - psF64 result = 0.0; - psF64 dut2ut1 = 0.0; - psF64 t = 0.0; - psVector *dut = NULL; - psMetadataItem *tableMetadataItem = NULL; - psLookupStatusType status = PS_LOOKUP_SUCCESS; - char *metadataTableNames[3] = {"ser7Table", "eopcTable", "finalsTable"}; - +psF64 psTimeGetUT1Delta(const psTime *time, psTimeBulletin bulletin) +{ + psU32 nTables = 2; + psF64 mjd = 0.0; + psF64 result = 0.0; + psU64 tableColumn = 0; + psF64 dut2ut1 = 0.0; + psF64 t = 0.0; + psVector* dut = NULL; + psMetadataItem* tableMetadataItem = NULL; + psLookupStatusType status = PS_LOOKUP_SUCCESS; + char* metadataTableNames[2] = {"dailyTable", "finalsTable"}; // Error checks PS_ASSERT_PTR_NON_NULL(time,NAN); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,1e9-1,NAN); - - if(time->type != PS_TIME_TAI) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_INCORRECT, time->type); - return 0.0; + PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN); + + // Check for invalid bulletin specified + if((bulletin != PS_IERS_A) && (bulletin != PS_IERS_B)) { + psError(PS_ERR_BAD_PARAMETER_VALUE,true,"Invalid bulletin specified %d",bulletin); + return NAN; + } + + // Set lookup table column based on Bullentin + if(bulletin == PS_IERS_A) { + tableColumn = 3; + } else { + tableColumn = 6; } // Attempt to find value through table lookup and interpolation mjd = psTimeToMJD(time); - result = searchTables(mjd, 0, &status, metadataTableNames, nTables); + result = searchTables(mjd,tableColumn,metadataTableNames,nTables,&status); // Value could not be found through table lookup and interpolation @@ -535,10 +842,12 @@ // Date too early for tables. Get default time delta value from metadata, and issue warning. - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES, mjd, "UT1-UTC"); - + psLogMsg(__func__,PS_LOG_WARN,PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES,mjd,"UT1-UTC"); + + // Lookup value from time metadata loaded from psTime.config tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.dut"); if(tableMetadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.before.dut"); - return 0.0; + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.before.dut"); + return NAN; } result = tableMetadataItem->data.F64; @@ -552,14 +861,18 @@ t = 2000.0 + (MJD - 51544.03)/365.2422 */ - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES, mjd, "UT1-UTC"); - + // Generate warning of postdate information + psLogMsg(__func__,PS_LOG_WARN,PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES, mjd, "UT1-UTC"); + + // Lookup values to calculate prediction tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.dut"); if(tableMetadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.dut"); - return 0.0; + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.predict.dut"); + return NAN; } dut = (psVector*)tableMetadataItem->data.V; - PS_ASSERT_PTR_NON_NULL(dut,0.0); - + PS_ASSERT_PTR_NON_NULL(dut,NAN); + + // Calculate predication of future UT1-UTC t = 2000.0 + (mjd - 51544.03)/365.2422; dut2ut1 = 0.022*sin(TWOPI*t) - 0.012*cos(TWOPI*t) - 0.006*sin(4.0*M_PI*t) + 0.007*cos(4.0*M_PI*t); @@ -568,5 +881,5 @@ } else if(status != PS_LOOKUP_SUCCESS) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED); - return 0.0; + return NAN; } @@ -574,5 +887,5 @@ } -struct psSphere* psTimeGetPoleCoords(const psTime* time) +struct psSphere* p_psTimeGetPoleCoords(const psTime* time) { psU32 nTables = 3; @@ -587,5 +900,5 @@ psLookupStatusType yStatus = PS_LOOKUP_SUCCESS; psMetadataItem *tableMetadataItem = NULL; - char *metadataTableNames[3] = {"ser7Table", "eopcTable", "finalsTable"}; + char *metadataTableNames[3] = {"dailyTable", "eopcTable", "finalsTable"}; psVector *xp = NULL; psVector *yp = NULL; @@ -593,5 +906,5 @@ // Error checks PS_ASSERT_PTR_NON_NULL(time,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,1e9-1,NULL); + PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL); if(time->type != PS_TIME_TAI) @@ -603,6 +916,8 @@ // Attempt to find value through table lookup and interpolation mjd = psTimeToMJD(time); - x = searchTables(mjd, 0, &xStatus, metadataTableNames, nTables); - y = searchTables(mjd, 0, &yStatus, metadataTableNames, nTables); + // x = searchTables(mjd, 0, &xStatus, metadataTableNames, nTables); + x = searchTables(mjd, 0, metadataTableNames, nTables,&xStatus); + // y = searchTables(mjd, 0, &yStatus, metadataTableNames, nTables); + y = searchTables(mjd, 0, metadataTableNames, nTables,&yStatus); // Value could not be found through table lookup and interpolation @@ -641,5 +956,4 @@ psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES, mjd, "polar motion"); - // Get predicted MJD tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.mjd"); @@ -692,5 +1006,5 @@ } - // Create output sphere and convert arcsec to radians (i.e. x/60/60*M_PI/180) + // Create output sphere and convert arcsec to radians (i.e. x/60/60*PS_PI/180) output = psAlloc(sizeof(psSphere)); output->r = x * M_PI / 648000.0; @@ -700,7 +1014,6 @@ } -double psTimeGetTAIDelta(const psTime *time) -{ - psU64 i = 0; +psF64 p_psTimeGetTAIDelta(const psTime *time) +{ psF64 jd = 0.0; psF64 mjd = 0.0; @@ -711,10 +1024,9 @@ psLookupTable* table = NULL; psMetadataItem *tableMetadataItem = NULL; - psVector *stats = NULL; psVector *results = NULL; // Error checks PS_ASSERT_PTR_NON_NULL(time,NAN); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,1e9-1,NAN); + PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN); // Check time metadata @@ -739,20 +1051,29 @@ mjd = psTimeToMJD(time); - stats = psVectorAlloc(table->numCols, PS_TYPE_U32); - results = psLookupTableInterpolateAll(table, jd, stats); + // Set ceiling of the julian date to the last entry in the lookup table + if(table->validTo < jd) { + jd = table->validTo; + } + + // Interpolation of look up table + results = psLookupTableInterpolateAll(table, jd); // Check for successful interpolation - for(i=0; in; i++) { - if(stats->data.U32[i] == PS_LOOKUP_ERROR) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED); - } - } - - const1 = results->data.F64[0]; - const2 = results->data.F64[1]; - const3 = results->data.F64[2]; - out = const1 + (mjd - const2) * const3; - - psFree(stats); + if(results == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED); + } + + // Set constants from table + const1 = results->data.F64[1]; + const2 = results->data.F64[2]; + const3 = results->data.F64[3]; + + // If const3 not equal to zero solve for difference else floor of const1 + if(fabs(const3-0.0) > FLT_EPSILON) { + out = const1 + (mjd - const2) * const3; + } else { + out = floor(const1); + } + psFree(results); @@ -760,33 +1081,64 @@ } - psS64 psTimeLeapSecondDelta(const psTime *time1, const psTime *time2) { psS64 diff = 0; - // Error checks PS_ASSERT_PTR_NON_NULL(time1,0); PS_ASSERT_PTR_NON_NULL(time2,0); - PS_ASSERT_INT_WITHIN_RANGE(time1->nsec,0,1e9-1,0); - PS_ASSERT_INT_WITHIN_RANGE(time2->nsec,0,1e9-1,0); - diff = abs((psS64)psTimeGetTAIDelta((psTime*)time1)-(psS64)psTimeGetTAIDelta((psTime*)time2)); + PS_ASSERT_INT_WITHIN_RANGE(time1->nsec,0,(psU32)((1e9)-1),0); + PS_ASSERT_INT_WITHIN_RANGE(time2->nsec,0,(psU32)((1e9)-1),0); + diff = abs((psS64)p_psTimeGetTAIDelta((psTime*)time1)-(psS64)p_psTimeGetTAIDelta((psTime*)time2)); return diff; } -double psTimeToJD(const psTime *time) -{ - double jd = 0.0; +psBool psTimeIsLeapSecond(const psTime* utc) +{ + psTime* prevUtc = NULL; + psBool returnValue = false; + + // Check for valid time + PS_ASSERT_PTR_NON_NULL(utc,false); + + // Verify time is UTC type + if(utc->type != PS_TIME_UTC) { + psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_TYPE_INCORRECT,utc->type); + return false; + } + + // Allocate time to hold utc - 1 second + prevUtc = psTimeAlloc(PS_TIME_UTC); + prevUtc->sec = utc->sec - 1; + + // Check the absolute difference between the two times for leapsecond + if(psTimeLeapSecondDelta(utc,prevUtc) >= 1.0) { + returnValue = true; + } else { + returnValue = false; + } + + // Free prevUtc + psFree(prevUtc); + + return returnValue; +} + +psF64 psTimeToJD(const psTime *time) +{ + psF64 jd = NAN; // Error checks PS_ASSERT_PTR_NON_NULL(time,NAN); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,1e9-1,NAN); + PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN); // Julian date conversion if(time->sec < 0) { - jd = time->sec / SEC_PER_DAY - time->nsec / NSEC_PER_DAY + 2440587.5; // psTime earlier than epoch + // psTime earlier than epoch + jd = time->sec / SEC_PER_DAY - time->nsec / NSEC_PER_DAY + JD_EPOCH_OFFSET; } else { - jd = time->sec / SEC_PER_DAY + time->nsec / NSEC_PER_DAY + 2440587.5; // psTime greater than epoch + // psTime greater than epoch + jd = time->sec / SEC_PER_DAY + time->nsec / NSEC_PER_DAY + JD_EPOCH_OFFSET; } @@ -794,18 +1146,19 @@ } -double psTimeToMJD(const psTime *time) -{ - double mjd = 0.0; - +psF64 psTimeToMJD(const psTime *time) +{ + psF64 mjd = NAN; // Error checks PS_ASSERT_PTR_NON_NULL(time,NAN); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,1e9-1,NAN); + PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN); // Modified Julian date conversion if(time->sec < 0) { - mjd = time->sec / SEC_PER_DAY - time->nsec / NSEC_PER_DAY + 40587.0; // psTime earlier than epoch + // psTime earlier than epoch + mjd = time->sec / SEC_PER_DAY - time->nsec / NSEC_PER_DAY + MJD_EPOCH_OFFSET; } else { - mjd = time->sec / SEC_PER_DAY + time->nsec / NSEC_PER_DAY + 40587.0; // psTime greater than epoch + // psTime greater than epoch + mjd = time->sec / SEC_PER_DAY + time->nsec / NSEC_PER_DAY + MJD_EPOCH_OFFSET; } @@ -815,5 +1168,5 @@ char* psTimeToISO(const psTime *time) { - psS32 ms = 0; + psS32 ds = 0; char *timeString = NULL; char *tempString = NULL; @@ -821,14 +1174,23 @@ time_t sec; - // Error checks PS_ASSERT_PTR_NON_NULL(time,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,1e9-1,NULL); - + PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL); + + // Check valid year range + PS_ASSERT_LONG_WITHIN_RANGE(time->sec,YEAR_0000_SEC,YEAR_9999_SEC,NULL) + + // Allocate temp strings tempString = psAlloc(MAX_TIME_STRING_LENGTH); timeString = psAlloc(MAX_TIME_STRING_LENGTH); - ms = time->nsec / 1000000; + // Convert nanoseconds to decaseconds + ds = time->nsec / 100000000; sec = time->sec; + + // If leapsecond use previous day + if(time->leapsecond) { + sec--; + } // tmTime variable is statically allocated, no need to free @@ -840,5 +1202,13 @@ } - if (snprintf(timeString, MAX_TIME_STRING_LENGTH, "%s.%3.3dZ", tempString, ms) < 0) { + // Check if time is UTC and leapsecond + if(((time->type==PS_TIME_UTC)||(time->type==PS_TIME_UT1)) && (time->leapsecond)) { + // Modify second to be 60 + tempString[17] = '6'; + tempString[18] = '0'; + } + + // Create string with milliseconds + if (snprintf(timeString, MAX_TIME_STRING_LENGTH, "%s,%1dZ", tempString, ds) < 0) { psError(PS_ERR_OS_CALL_FAILED, true, PS_ERRORTEXT_psTime_APPEND_MSEC_FAILED); } @@ -848,22 +1218,25 @@ } -struct timeval psTimeToTimeval(const psTime *time) -{ - struct timeval timevalTime; - - timevalTime.tv_sec = 0; - timevalTime.tv_usec = 0; +struct timeval* psTimeToTimeval(const psTime *time) +{ + struct timeval *timevalTime = NULL; // Error checks - PS_ASSERT_PTR_NON_NULL(time,timevalTime); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,1e9-1,timevalTime); - - timevalTime.tv_sec = time->sec; - timevalTime.tv_usec = time->nsec / 1000; + PS_ASSERT_PTR_NON_NULL(time,NULL); + PS_ASSERT_INT_WITHIN_RANGE(time->sec,0,INT32_MAX,NULL); + PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL); + + // Allocate structure timeval + timevalTime = (struct timeval*)psAlloc(sizeof(struct timeval)); + + // Set structure members + timevalTime->tv_sec = time->sec; + timevalTime->tv_usec = time->nsec / 1000; return timevalTime; } -struct tm* psTimeToTM(const psTime *time) +/* +struct tm* p_psTimeToTM(const psTime *time) { psS64 cent = 0; @@ -876,15 +1249,15 @@ psS64 temp = 0; struct tm* tmTime = NULL; - - + + // Error checks PS_ASSERT_PTR_NON_NULL(time,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,1e9-1,NULL); - + PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL); + seconds = time->sec%60; minute = time->sec/60%60; hour = time->sec/3600%24; day = (time->sec+62135596800)/86400; - + // Add 306 days to make relative to Mar 1, 0; also adjust day to be within a range (1..2**28-1) where our // calculations will work with 32bit ints @@ -898,5 +1271,5 @@ day -= temp * 146097; } - + cent = (day*4-1)/146097; // Calc number of centuries day is after 29 Feb of yr 0 day -= cent*146097/4; // 4 centuries = 146097 days @@ -906,5 +1279,5 @@ day -= (month*367-1094)/12; // February of following year) year += cent*100+temp*400; // Get the real year, which is off by - + // One if month is January or February if(month > 12) @@ -913,8 +1286,8 @@ month -= 12; } - + // Allocate output tmTime = (struct tm*)psAlloc(sizeof(struct tm)); - + tmTime->tm_year = year - 1900; tmTime->tm_mon = month - 1; @@ -924,14 +1297,14 @@ tmTime->tm_sec = seconds; tmTime->tm_isdst = -1; - + return tmTime; } - -psTime* psTimeFromJD(double time) -{ - double days = 0.0; - double seconds = 0.0; +*/ + +psTime* psTimeFromJD(psF64 time) +{ + psF64 days = 0.0; + psF64 seconds = 0.0; psTime *outTime = NULL; - // Allocate psTime struct @@ -949,15 +1322,14 @@ // Error check - PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,1e9-1,outTime); + PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),outTime); return outTime; } -psTime* psTimeFromMJD(double time) -{ - double days = 0.0; - double seconds = 0.0; +psTime* psTimeFromMJD(psF64 time) +{ + psF64 days = 0.0; + psF64 seconds = 0.0; psTime *outTime = NULL; - // Allocate psTime struct @@ -967,4 +1339,5 @@ days = time - 40587.0; seconds = days * SEC_PER_DAY; + if(seconds < 0.0) { outTime->nsec = (seconds - (psS64)seconds) * -1000000000.0; // psTime earlier than epoch @@ -975,5 +1348,5 @@ // Error check - PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,1e9-1,outTime); + PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),NULL); return outTime; @@ -986,6 +1359,9 @@ psTime *outTime = NULL; + // Check for NULL string + PS_ASSERT_PTR_NON_NULL(time,NULL); + // Convert YYYY-MM-DDThh:mm:ss.sss in string form to tm time - if (sscanf(time, "%d-%d-%dT%d:%d:%d.%d", &tmTime.tm_year, &tmTime.tm_mon, &tmTime.tm_mday, + if (sscanf(time, "%d-%d-%dT%d:%d:%d,%d", &tmTime.tm_year, &tmTime.tm_mon, &tmTime.tm_mday, &tmTime.tm_hour, &tmTime.tm_min, &tmTime.tm_sec,&millisecond) < 7) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_ISOTIME_MALFORMED, time); @@ -1006,6 +1382,44 @@ // Convert tm time to psTime - outTime = psTimeFromTM(&tmTime); + outTime = p_psTimeFromTM(&tmTime); outTime->nsec = millisecond * 1000000; + + return outTime; +} + +psTime* psTimeFromTT(psS64 sec, psU32 nsec) +{ + psTime* outTime = NULL; + + // Verify nsec within range + PS_ASSERT_INT_WITHIN_RANGE(nsec,0,(psU32)((1e9)-1),NULL); + + // Allocate psTime data + outTime = psTimeAlloc(PS_TIME_TT); + + // Set data members + outTime->sec = sec; + outTime->nsec = nsec; + + // Return data structure + return outTime; +} + +psTime* psTimeFromUTC(psS64 sec, psU32 nsec, psBool leapsecond) +{ + psTime* outTime = NULL; + + // Verify nsec within range + PS_ASSERT_INT_WITHIN_RANGE(nsec,0,(psU32)((1e9)-1),NULL); + + // Allocate psTime data + outTime = psTimeAlloc(PS_TIME_UTC); + + // Set data members + outTime->sec = sec; + outTime->nsec = nsec; + + // Set leapsecond flag if necessary + outTime->leapsecond = psTimeIsLeapSecond(outTime); return outTime; @@ -1028,10 +1442,10 @@ // Error check - PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,1e9-1,outTime); + PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),outTime); return outTime; } -psTime* psTimeFromTM(const struct tm* time) +psTime* p_psTimeFromTM(const struct tm* time) { psS64 year; @@ -1044,5 +1458,4 @@ psTime *outTime = NULL; - // Error check PS_ASSERT_PTR_NON_NULL(time,NULL); @@ -1062,9 +1475,13 @@ if( month <= 2 ) { - year -= (temp = (14 - month) / 12); + temp = (14 - month) / 12; + // year -= (temp = (14 - month) / 12); + year -= temp; month += 12 * temp; } else if(month > 14) { - year += (temp = (month - 3) / 12); + temp = (month - 3) / 12; + // year += (temp = (month - 3) / 12); + year += temp; month -= 12 * temp; } @@ -1081,5 +1498,6 @@ // days to adjust from Mar 1, year 0-relative to Jan 1, year 1-relative. Add hours, minutes, and seconds. day += (month * 367 - 1094) / 12 + year % 100 * 1461 / 4 + (year/100 * 36524 + year/400) - 306; - outTime->sec = (((day - 1) * SEC_PER_DAY) - 62135596800) + hour*SEC_PER_HOUR + minute*SEC_PER_MINUTE + seconds; + outTime->sec = (((day - 1) * SEC_PER_DAY) - 62135596800) + hour*SEC_PER_HOUR + + minute*SEC_PER_MINUTE + seconds; // C's TM does not define a microsecond field. Microseconds must be manipulated by calling function. @@ -1087,5 +1505,5 @@ // Error check - PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,1e9-1,outTime); + PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),outTime); return outTime; @@ -1098,8 +1516,7 @@ psTime *tempTime = NULL; - // Error checks PS_ASSERT_PTR_NON_NULL(time,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,1e9-1,NULL); + PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL); // Convert time to TAI if necessary, but without changing input arguments @@ -1109,16 +1526,18 @@ tempTime->nsec = time->nsec; tempTime = psTimeConvert(tempTime, PS_TIME_TAI); + outTime = psTimeAlloc(PS_TIME_TAI); } else { tempTime = psMemIncrRefCounter((psTime*)time); + outTime = psTimeAlloc(time->type); } // Create output time - outTime = psTimeAlloc(PS_TIME_TAI); sec = delta + (psF64)tempTime->sec + (psF64)tempTime->nsec/1e9; - outTime->sec = sec; - outTime->nsec = (sec - outTime->sec)*1e9; + PS_ASSERT_LONG_WITHIN_RANGE((psS64)sec,0,PS_MAX_S64,outTime); + outTime->sec = (psS64)sec; + outTime->nsec = (psU32)((sec - (psF64)outTime->sec)*1e9); // Error check - PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,1e9-1,outTime); + PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),outTime); // Convert result to same time type as input @@ -1140,10 +1559,15 @@ psTime *tempTime2 = NULL; - // Error checks PS_ASSERT_PTR_NON_NULL(time1,0.0); - PS_ASSERT_INT_WITHIN_RANGE(time1->nsec,0,1e9-1,0.0); + PS_ASSERT_INT_WITHIN_RANGE(time1->nsec,0,(psU32)((1e9)-1),0.0); PS_ASSERT_PTR_NON_NULL(time2,0.0); - PS_ASSERT_INT_WITHIN_RANGE(time2->nsec,0,1e9-1,0.0); + PS_ASSERT_INT_WITHIN_RANGE(time2->nsec,0,(psU32)((1e9)-1),0.0); + + // Verify both times of the same type + if(time1->type != time2->type) { + psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_TYPE_INCORRECT,time1->type); + return out; + } // Convert time to TAI if necessary, but without changing input arguments