Index: trunk/psLib/src/astro/psTime.c =================================================================== --- trunk/psLib/src/astro/psTime.c (revision 2273) +++ trunk/psLib/src/astro/psTime.c (revision 2301) @@ -5,12 +5,11 @@ * * A collection of functions are required by psLib to manipulate time data. These functions primarily consist - * of conversions between specific time formats. PSLib currently uses the UNIX timeval time system as the - * base upon which International Atomic Time (TAI) and Universal Time Coordinated (UTC) are calculated. TAI - * time varies over time due to the earth's rotation and the movement of the continental plates. + * of conversions between specific time formats. They use the UNIX timeval time system as the + * base upon which International Atomic Time (TAI) and Universal Time Coordinated (UTC) are calculated. * * @author Ross Harman, MHPCC * - * @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-11-04 01:04:57 $ + * @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-11-09 00:36:13 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -21,4 +20,5 @@ #include #include +#include #include "psTime.h" @@ -31,17 +31,12 @@ #include "psMetadata.h" #include "psMetadataIO.h" +#include "psLookupTable.h" #include "psConstants.h" #include "psAstronomyErrors.h" -#ifndef SER7_FILE -#define SER7_FILE "../../data/ser7.dat" -#pragma warning SER7_FILE was not defined in the makefile. +#ifndef TIME_CONFIG_FILE +#define TIME_CONFIG_FILE "../../config/time.config" +#pragma warning TIME_CONFIG_FILE was not defined in the makefile. #endif - -#ifndef TAIUTC_FILE -#define TAIUTC_FILE "../../data/tai-utc.dat" -#pragma warning TAIUTC_FILE was not defined in the makefile. -#endif - /** Sidereal angular conversion from seconds to radians for GMST in seconds (i.e. pi/(180*240)) */ @@ -72,361 +67,245 @@ #define USEC_PER_DAY 86400000000.0 - /** Max size of a line read */ - #define LINESIZE 1024 - - /** Time metadata read from config file */ static psMetadata *timeMetadata = NULL; - -/** Read TAI/UTC data file. - * - * Reads TAI/UTC data file. - * - * @return psImage*: Two-dimensional array with table data. - */ -static psImage* readTaiUtcFile( - const char *fileName ///< Name of file to read. -); - - -/** Read SER7 data file. - * - * Reads SER7 data file. - * - * @return psImage*: Two-dimensional array with table data. - */ -static psImage* readSer7File( - const char *fileName ///< Name of file to read. -); - - -/** Read finals data file. - * - * Reads finals data file. - * - * @return psImage*: Two-dimensional array with table data. - */ -/* not yet used -static psImage* readFinalsFile( - char *fileName ///< Name of file to read. -); -*/ - -/** Read EOPC data file. - * - * Reads EOPC data file. - * - * @return psImage*: Two-dimensional array with table data. - */ -/* not yet used -static psImage* readEopcFile( - char *fileName ///< Name of file to read. -); -*/ - -/** Lookup and interpolate TAI-UTC data - * - * Interpolates TAI-UTC data, given time in MJD format. - * - * @return double: Interpolated value. - */ -static double lookupTaiUtcTable( - const psTime *time ///< time to lookup. -); - -/** Lookup and interpolate polar or UT1-UTC data - * - * Interpolates TAI-UT1 or polar coordinate data, given time in MJD format. - * - * @return double: Interpolated value. - */ -static double lookupSer7Table( - const psTime *time, ///< time to lookup. - psS32 col ///< Column to lookup. -); - -static psImage* readTaiUtcFile(const char *fileName) -{ - char line[LINESIZE]; - psS32 j = 0; - psS32 maxLines = 1000; - psF64 *ptr = NULL; - psImage *table = NULL; - FILE *fd = NULL; - - - fd=fopen(fileName, "r"); - if(fd==NULL) { - psError(PS_ERR_IO,true, - PS_ERRORTEXT_psTime_FILE_NOT_FOUND, - fileName,"Tai-Utc"); - return NULL; - } - - // Table shouldn't be larger than 100 rows. Select columns are read based on indices. - table = psImageAlloc(4, maxLines, PS_TYPE_F64); - - while(fgets(line, LINESIZE, fd) != NULL) { - if(j>=maxLines) { - psError(PS_ERR_IO,true, - PS_ERRORTEXT_psTime_FILE_TOO_MANY_ROWS, - fileName,maxLines); - return NULL; - psFree(table); - } - - ptr = table->data.F64[j++]; - *ptr = atof(line+16); // MJD - *(ptr+1) = atof(line+38); // TAI-UTC - *(ptr+2) = atof(line+59); // Contant #1 - *(ptr+3) = atof(line+69); // Contant #2 - } - - *(psU32 *)&table->numRows = j; - - p_psMemSetPersistent(table,true); - p_psMemSetPersistent(table->data.V,true); - p_psMemSetPersistent(table->rawDataBuffer,true); - - return table; -} - -static psImage* readSer7File(const char *fileName) -{ - psBool beginRecord = false; - char line[LINESIZE]; - psS32 j = 0; - psS32 maxLines = 400; - psImage *table = NULL; - FILE *fd = NULL; - - - fd=fopen(fileName, "r"); - if(fd==NULL) { - psError(PS_ERR_IO,true, - PS_ERRORTEXT_psTime_FILE_NOT_FOUND, - fileName,"Ser7"); - return NULL ; - } - - // Table shouldn't be larger than 400 rows. All columns in area of interest are read. - table = psImageAlloc(7, maxLines, PS_TYPE_F64); - - while(fgets(line, LINESIZE, fd) != NULL) { - - // Stop parsing after finding "These" in data file after read has began - if(beginRecord) { - if(strstr(line, "These")!=NULL) { - beginRecord = false; - break; - } - - if(j>=maxLines) { - psError(PS_ERR_IO,true, - PS_ERRORTEXT_psTime_FILE_TOO_MANY_ROWS, - fileName,maxLines); - return NULL; - psFree(table); - } - psF64 *ptr = table->data.F64[j++]; - sscanf(line, "%lf %lf %lf %lf %lf %lf %lf", ptr, ptr+1, ptr+2, ptr+3, ptr+4, ptr+5, ptr+6); - } - - // Start parsing after finding "x(arcsec)" in data file - if(strstr(line, "x(arcsec)") != NULL) { - beginRecord = true; - } - } - *(psU32 *)&table->numRows = j; - - p_psMemSetPersistent(table,true); - p_psMemSetPersistent(table->data.V,true); - p_psMemSetPersistent(table->rawDataBuffer,true); - - return table; -} - - -static double lookupTaiUtcTable(const psTime *time) -{ - psS32 hiIdx = 0; - psS32 loIdx = 0; - psS32 numRows = 0; - double jd = 0.0; - double mjd = 0.0; - double out = 0.0; - double denom = 0.0; - double const1 = 0.0; - double const2 = 0.0; - double const3 = 0.0; - static psImage* table = NULL; - - if (table == NULL) { - table = readTaiUtcFile(TAIUTC_FILE); - if (table == NULL) { - psAbort(__func__,"Failed to read tai-utc.dat file"); - } - } - - // Number of rows is the number of rows of data read from the data file - numRows = table->numRows; - - // Determine Julian and modified Julian dates used in table lookup and time delta calculation - jd = psTimeToJD(time); - mjd = psTimeToMJD(time); - - // Variable in should be in MJD format - if(jd < table->data.F64[0][0]) { // MJD time before start of table - return 0.0; - } else if(jd > table->data.F64[numRows-1][0]) { // MJD time after end of table - return table->data.F64[numRows-1][1]; - } else { // All other times..interpolate - while(jd > table->data.F64[hiIdx][0]) { - hiIdx++; - if(hiIdx >= numRows) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true, - PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLE, - jd, "TAI-UTC"); - return 0.0; - } - } - loIdx = hiIdx--; - if(loIdx < 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true, - PS_ERRORTEXT_psTime_TIME_PREDATES_TABLE, - jd, "TAI-UTC"); +/** 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); + +/** Removes leading and trailing whitespace and # characters from a string. The cleaned string is a new null + * terminated copy of the original input string. */ +static char *cleanString(char *inString, int sLen) +{ + char *ptrB = NULL; + char *ptrE = NULL; + char *cleaned = NULL; + + + ptrB = inString; + + /* Skip over leading # or whitespace */ + while (isspace(*ptrB) || *ptrB=='#') { + ptrB++; + } + + /* Skip over trailing whitespace, null terminators, and # characters */ + ptrE = inString + sLen; + while(isspace(*ptrE) || *ptrE=='\0' || *ptrE=='#') { + ptrE--; + } + + // Length, sLen, does not include '\0' + sLen = ptrE - ptrB + 1; + + // Adds '\0' to end of string and +1 to sLen + cleaned = psStringNCopy(ptrB, sLen); + + return cleaned; +} + +/** Returns cleaned token based on delimiter, but not including delimiter. Also changes the pointer location + * the beginning of the string. Tokens are newly allocated null terminated strings. */ +static char* getToken(char **inString, char *delimiter, psParseErrorType *status) +{ + char *cleanToken = NULL; + int sLen = 0; + + + // Skip over leading whitespace + while(isspace(**inString)) { + (*inString)++; + } + + // Length of token, not including delimiter + sLen = strcspn(*inString, delimiter); + if(sLen) { + + // Create new, cleaned, and null terminated token + cleanToken = cleanString(*inString, sLen); + + // Move to end of token + (*inString) += sLen; + } else if(**inString!='\0' && sLen==0) { + *status = PS_PARSE_ERROR_GENERAL; + } + + return cleanToken; +} + +/** 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; + + + + // Check time metadata. Function call reports errors. + if(timeMetadata == NULL) { + if(!psTimeInit(TIME_CONFIG_FILE)) return 0.0; - } - - // since rows of table can not have same jd, this check is probably not needed. - denom = table->data.F64[hiIdx][0]-table->data.F64[loIdx][0]; - if(fabs(denom) < FLT_EPSILON) { - psError(PS_ERR_UNKNOWN,true, - PS_ERRORTEXT_psTime_TABLE_DUPLICATE_ROWS - "Ser7"); + } + + // Search each table in priority order: ser7, eopc,finals + for(i=0; idata.F64[loIdx][1]+(table->data.F64[hiIdx][1]-table->data.F64[loIdx][1]) - *(jd-table->data.F64[loIdx][0])/denom; - - const2 = table->data.F64[loIdx][2]+(table->data.F64[hiIdx][2]-table->data.F64[loIdx][2]) - *(jd-table->data.F64[loIdx][0])/denom; - - const3 = table->data.F64[loIdx][3]+(table->data.F64[hiIdx][3]-table->data.F64[loIdx][3]) - *(jd-table->data.F64[loIdx][0])/denom; - } - } - - out = const1 + (mjd - const2) * const3; - - return out; -} - -static double lookupSer7Table(const psTime *time, psS32 col) -{ - psS32 hiIdx = 0; - psS32 loIdx = 0; - psS32 numRows = 0; - double out = 0.0; - double denom = 0.0; - double mjdUtc = 0.0; - static psImage* table = NULL; - - if (table == NULL) { - table = readSer7File(SER7_FILE); - if (table == NULL) { - psAbort(__func__,"Failed to read the ser7.dat file."); - } - } - - // Number of rows is the number of rows of data read from the data file - numRows = table->numRows; - - // Interpolation value in should be in MJD format - mjdUtc = psTimeToMJD(time); - if(mjdUtc < table->data.F64[0][3]) { // MJD time before start of table - psError(PS_ERR_BAD_PARAMETER_VALUE,true, - PS_ERRORTEXT_psTime_TIME_PREDATES_TABLE, - mjdUtc, "Ser7"); - - } else if(mjdUtc > table->data.F64[numRows-1][3]) { // MJD time after end of table - psError(PS_ERR_BAD_PARAMETER_VALUE,true, - PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLE, - mjdUtc, "Ser7"); - } else { // All other times..interpolate - while(mjdUtc > table->data.F64[hiIdx][3]) { - hiIdx++; - if(hiIdx >= numRows) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true, - PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLE, - mjdUtc, "Ser7"); - return 0.0; - } - } - - loIdx = hiIdx--; - if(loIdx < 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true, - PS_ERRORTEXT_psTime_TIME_PREDATES_TABLE, - mjdUtc, "Ser7"); - } - - denom = table->data.F64[hiIdx][3] - table->data.F64[loIdx][3]; - if(fabs(denom) < FLT_EPSILON) { - psError(PS_ERR_UNKNOWN,true, - PS_ERRORTEXT_psTime_TABLE_DUPLICATE_ROWS, - "Ser7"); - } else { - out = table->data.F64[loIdx][col]+(table->data.F64[hiIdx][col]-table->data.F64[loIdx][col]) - *(mjdUtc-table->data.F64[loIdx][3])/denom; - } - } - - return out; -} - -void psTimeInit(char *fileName) -{ + } + table = (psLookupTable*)tableMetadataItem->data.V; + PS_PTR_CHECK_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; + } + } + + return result; +} + +bool psTimeInit(char *fileName) +{ + bool foundTable = false; char *tableDir = NULL; char *tableNames = NULL; - psU8 numTables = 0; + char *namesPtr = NULL; + char *metadataNamesPtr = NULL; + char *tableName = NULL; + char *fullTableName = NULL; + psS32 i = 0; + psS32 j = 0; + psS32 numTables = 0; psVector *tablesFrom = NULL; psVector *tablesTo = NULL; - psF64 beforeXp = 0.0; - psF64 beforeYp = 0.0; - psF64 beforeDut = 0.0; psMetadataItem *metadataItem = NULL; + psLookupTable *table = NULL; + psLookupStatusType status = PS_LOOKUP_SUCCESS; + char *metadataTableNames[4] = {"ser7", "eopc", "finals", "tai"}; // Read time config file - psMetadataParseConfig(&timeMetadata, fileName, true); + if(psMetadataParseConfig(&timeMetadata, fileName, true) < 0) { + return false; + } // Get path to time data files metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.dir"); - tableDir = (char*)metadataItem->data.V; + if(metadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.tables.dir"); + return false; + } + tableDir = psStringCopy(metadataItem->data.V); // Get number of tables metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.n"); - numTables = (psU8)metadataItem->data.S32; - - // Get table names + if(metadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.tables.n"); + return false; + } + numTables = (psS32)metadataItem->data.S32; + + // Table file names metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.files"); - tableNames = (char*)metadataItem->data.V; - - // Get ranges of tables + if(metadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.tables.files"); + return false; + } + tableNames = psStringCopy(metadataItem->data.V); + + // Get lower range of tables metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.from"); - tablesFrom = (psVector*)metadataItem->data.V; + if(metadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.tables.from"); + return false; + } + tablesFrom = psVectorCopy(tablesFrom, metadataItem->data.V, PS_TYPE_F64); + if(tablesFrom->n != numTables) { + return false; + } + + // Get upper range of tables metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.to"); - tablesTo = (psVector*)metadataItem->data.V; - - // Out of range values - metadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.xp"); - beforeXp =(psF64)metadataItem->data.F64; - metadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.yp"); - beforeYp =(psF64)metadataItem->data.F64; - metadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.dut"); - beforeDut = (psF64)metadataItem->data.F64; + if(metadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.tables.to"); + return false; + } + tablesTo = psVectorCopy(tablesTo, metadataItem->data.V, PS_TYPE_F64); + if(tablesTo->n != numTables) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_VECTOR, tablesTo->n, numTables); + return false; + } + + // Read time tables + namesPtr = tableNames; + while((tableName=getToken(&namesPtr, " ", &status)) != NULL) { + fullTableName = (char*)psAlloc(strlen(tableDir)+strlen(tableName)+1+1); // Add one for last dir "/" + fullTableName[0]='\0'; // Old strings may come back from psAlloc, so set init pos to EOL + strcat(fullTableName, tableDir); + strcat(fullTableName, "/"); + strcat(fullTableName, tableName); + + if(i < numTables) { + table = psLookupTableAlloc(fullTableName, tablesFrom->data.F64[i], tablesTo->data.F64[i]); + psLookupTableRead(table); + } else { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_TABLE_COUNT, i+1, numTables); + } + + // Place tables into metadata slightly altered names as keys to create consistent naming conventions + foundTable = false; + for(j=0; j i+1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_TABLE_COUNT, i+1, numTables); + } + + psFree(tableDir); + psFree(tableNames); + psFree(tablesFrom); + psFree(tablesTo); + + return true; +} + +void psTimeFinalize(void) +{ + psFree(timeMetadata); } @@ -552,6 +431,5 @@ // Convert Universal Time (UTC) to UT1 ut1UtcDelta = psTimeAlloc(PS_TIME_UTC); - - ut1UtcDbl = psTimeGetUT1Delta(utcTime)*1e6; + ut1UtcDbl = psTimeGetUT1Delta(taiTime)*1e6; ut1UtcDelta->usec = (psU32)fabs(ut1UtcDbl); if(ut1UtcDbl<0) { @@ -561,5 +439,4 @@ } - // Calculate UT1 as Julian Centuries since J2000.0 jdUt1Days = psTimeToJD(ut1Time); @@ -580,5 +457,5 @@ fracDays = fmod(mjdUt1Days, 1.0); - // Calculate Greenwich Mean Sidereal Time (GMST) in radians. Equation set up to minimize multiplications + // 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; @@ -602,4 +479,14 @@ 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] = {"ser7", "eopc", "finals"}; + // Error checks @@ -607,13 +494,72 @@ PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NAN); - return lookupSer7Table(time, 6); + if(time->type != PS_TIME_TAI) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_INCORRECT, time->type); + return 0.0; + } + + // Attempt to find value through table lookup and interpolation + mjd = psTimeToMJD(time); + result = searchTables(mjd, 0, &status, metadataTableNames, nTables); + + // Value could not be found through table lookup and interpolation + if(status == PS_LOOKUP_PAST_TOP) { + + // 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"); + + 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; + } + result = tableMetadataItem->data.F64; + + } else if(status == PS_LOOKUP_PAST_BOTTOM) { + /* Date too late for tables. Issue warning and use following formulae for predicting + ahead of the most recent available table entry. + ut1-utc = [0] + [1]*(MJD - [2]) - (ut2-ut1) + [0, 1, 2] = @psLib.time.predict.dut + ut2-ut1 = 0.022 sin(2*pi*t) - 0.012 cos(2*pi*t) - 0.006 sin(4*pi*t) + 0.007 cos(4*pi*t) + t = 2000.0 + (MJD - 51544.03)/365.2422 + */ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES, mjd, "UT1-UTC"); + + 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; + } + dut = (psVector*)tableMetadataItem->data.V; + PS_PTR_CHECK_NULL(dut,0.0); + + 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); + result = dut->data.F64[0] + dut->data.F64[1]*(mjd - dut->data.F64[2]) - dut2ut1; + + } else if(status != PS_LOOKUP_SUCCESS) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED); + return 0.0; + } + + return result; } struct psSphere* psTimeGetPoleCoords(const psTime* time) { - double x = 0.0; - double y = 0.0; + psU32 nTables = 3; + psF64 x = 0.0; + psF64 y = 0.0; + psF64 mjd = 0.0; + psF64 a = 0.0; + psF64 c = 0.0; + psF64 mjdPred = 0.0; struct psSphere* output = NULL; - + psLookupStatusType xStatus = PS_LOOKUP_SUCCESS; + psLookupStatusType yStatus = PS_LOOKUP_SUCCESS; + psMetadataItem *tableMetadataItem = NULL; + char *metadataTableNames[3] = {"ser7", "eopc", "finals"}; + psVector *xp = NULL; + psVector *yp = NULL; // Error checks @@ -621,9 +567,104 @@ PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NULL); - x = lookupSer7Table((psTime*)time, 4); - y = lookupSer7Table((psTime*)time, 5); - + if(time->type != PS_TIME_TAI) + { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_INCORRECT, time->type); + return NULL; + } + + // 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); + + // Value could not be found through table lookup and interpolation + if(xStatus==PS_LOOKUP_PAST_TOP && yStatus==PS_LOOKUP_PAST_TOP) + { + + // Date too earlier for tables. Get default polar coodinate values from metadata, and issue warning. + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES, mjd, "polar motion"); + + tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.xp"); + if(tableMetadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.before.xp"); + return NULL; + } + x = tableMetadataItem->data.F64; + + tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.yp"); + if(tableMetadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.before.yp"); + return NULL; + } + y = tableMetadataItem->data.F64; + + } else if(xStatus==PS_LOOKUP_PAST_BOTTOM && yStatus==PS_LOOKUP_PAST_BOTTOM) + { + + /* Date too late for tables. Issue warning and use following formulae for predicting + ahead of the most recent available table entry. + x = [0] + [1]*cos a + [2]*sin a + [3]*cos c + [4]*sin c + [0], [1], [2], [3] = @psLib.time.predict.xp + y = [0] + [1]*cos a + [2]*sin a + [3]*cos c + [4]*sin c + [0], [1], [2], [3] = @psLib.time.predict.yp + a = 2*pi*(mjd - pslib.time.predict.mjd)/365.25 + c = 2*pi*(mjd - pslib.time.predict.mjd)/435.0 + */ + 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"); + if(tableMetadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.predict.mjd"); + return NULL; + } + mjdPred = tableMetadataItem->data.F64; + + // Get xp + tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.xp"); + if(tableMetadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.xp"); + return NULL; + } + xp = (psVector*)tableMetadataItem->data.V; + PS_PTR_CHECK_NULL(xp,NULL); + + // Get yp + tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.yp"); + if(tableMetadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.yp"); + return NULL; + } + yp = (psVector*)tableMetadataItem->data.V; + PS_PTR_CHECK_NULL(yp,NULL); + + // Calculate "a" and "c" constants + a = TWOPI*(mjd - mjdPred)/365.25; + c = TWOPI*(mjd - mjdPred)/435.0; + + // Calculate x and y polar coordinates + x = xp->data.F64[0] + + xp->data.F64[1]*cos(a) + + xp->data.F64[2]*sin(a) + + xp->data.F64[3]*cos(c) + + xp->data.F64[4]*sin(c); + + y = yp->data.F64[0] + + yp->data.F64[1]*cos(a) + + yp->data.F64[2]*sin(a) + + yp->data.F64[3]*cos(c) + + yp->data.F64[4]*sin(c); + + } else if(xStatus!=PS_LOOKUP_SUCCESS || yStatus!=PS_LOOKUP_SUCCESS) + { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED); + return NULL; + } + + // Create output sphere and convert arcsec to radians (i.e. x/60/60*M_PI/180) output = psAlloc(sizeof(psSphere)); - output->r = x * M_PI / 648000.0; // x in arcsec converted to radians, i.e., x/60/60*M_PI/180 + output->r = x * M_PI / 648000.0; output->d = y * M_PI / 648000.0; @@ -633,4 +674,15 @@ double psTimeGetTAIDelta(const psTime *time) { + psU64 i = 0; + psF64 jd = 0.0; + psF64 mjd = 0.0; + psF64 out = 0.0; + psF64 const1 = 0.0; + psF64 const2 = 0.0; + psF64 const3 = 0.0; + psLookupTable* table = NULL; + psMetadataItem *tableMetadataItem = NULL; + psVector *stats = NULL; + psVector *results = NULL; // Error checks @@ -638,6 +690,45 @@ PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NAN); - return lookupTaiUtcTable(time); -} + // Check time metadata + if(timeMetadata == NULL) { + if(!psTimeInit(TIME_CONFIG_FILE)) + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_FILE_NOT_FOUND, "psTime.config"); + return 0.0; + } + + // Get table from metadata + tableMetadataItem = psMetadataLookup(timeMetadata, "taiTable"); + if(tableMetadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "taiTable"); + return 0.0; + } + table = (psLookupTable*)tableMetadataItem->data.V; + PS_PTR_CHECK_NULL(table,0); + + // Determine Julian and modified Julian dates used in table lookup and time delta calculation + jd = psTimeToJD(time); + mjd = psTimeToMJD(time); + + stats = psVectorAlloc(table->numCols, PS_TYPE_U32); + results = psLookupTableInterpolateAll(table, jd, stats); + + // Check for successful interpolation + for(i=0; in; i++) { + if(stats->data.U32[i]!=PS_LOOKUP_PAST_BOTTOM && stats->data.U32[i]!=PS_LOOKUP_PAST_TOP) { + 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); + psFree(results); + + return out; +} + psS64 psTimeLeapseconds(const psTime *time1, const psTime *time2)