Index: trunk/psLib/src/astronomy/psTime.c =================================================================== --- trunk/psLib/src/astronomy/psTime.c (revision 1614) +++ trunk/psLib/src/astronomy/psTime.c (revision 1788) @@ -7,53 +7,10 @@ * 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. It is - * currentls 32 seconds faster than UTC/timeval time. The conversion between TAI and UTC time for other dates - * is given in the table below, which comes from: - * - * ftp://maia.usno.navy.mil/ser7/tai-utc.dat - * - * This table is placed into a static array and used by the psTime functions for conversions. - * - * 1961 JAN 1 =JD 2437300.5 TAI-UTC= 1.4228180 S + (MJD - 37300.) X 0.001296 S - * 1961 AUG 1 =JD 2437512.5 TAI-UTC= 1.3728180 S + (MJD - 37300.) X 0.001296 S - * 1962 JAN 1 =JD 2437665.5 TAI-UTC= 1.8458580 S + (MJD - 37665.) X 0.0011232S - * 1963 NOV 1 =JD 2438334.5 TAI-UTC= 1.9458580 S + (MJD - 37665.) X 0.0011232S - * 1964 JAN 1 =JD 2438395.5 TAI-UTC= 3.2401300 S + (MJD - 38761.) X 0.001296 S - * 1964 APR 1 =JD 2438486.5 TAI-UTC= 3.3401300 S + (MJD - 38761.) X 0.001296 S - * 1964 SEP 1 =JD 2438639.5 TAI-UTC= 3.4401300 S + (MJD - 38761.) X 0.001296 S - * 1965 JAN 1 =JD 2438761.5 TAI-UTC= 3.5401300 S + (MJD - 38761.) X 0.001296 S - * 1965 MAR 1 =JD 2438820.5 TAI-UTC= 3.6401300 S + (MJD - 38761.) X 0.001296 S - * 1965 JUL 1 =JD 2438942.5 TAI-UTC= 3.7401300 S + (MJD - 38761.) X 0.001296 S - * 1965 SEP 1 =JD 2439004.5 TAI-UTC= 3.8401300 S + (MJD - 38761.) X 0.001296 S - * 1966 JAN 1 =JD 2439126.5 TAI-UTC= 4.3131700 S + (MJD - 39126.) X 0.002592 S - * 1968 FEB 1 =JD 2439887.5 TAI-UTC= 4.2131700 S + (MJD - 39126.) X 0.002592 S - * 1972 JAN 1 =JD 2441317.5 TAI-UTC= 10.0 S + (MJD - 41317.) X 0.0 S - * 1972 JUL 1 =JD 2441499.5 TAI-UTC= 11.0 S + (MJD - 41317.) X 0.0 S - * 1973 JAN 1 =JD 2441683.5 TAI-UTC= 12.0 S + (MJD - 41317.) X 0.0 S - * 1974 JAN 1 =JD 2442048.5 TAI-UTC= 13.0 S + (MJD - 41317.) X 0.0 S - * 1975 JAN 1 =JD 2442413.5 TAI-UTC= 14.0 S + (MJD - 41317.) X 0.0 S - * 1976 JAN 1 =JD 2442778.5 TAI-UTC= 15.0 S + (MJD - 41317.) X 0.0 S - * 1977 JAN 1 =JD 2443144.5 TAI-UTC= 16.0 S + (MJD - 41317.) X 0.0 S - * 1978 JAN 1 =JD 2443509.5 TAI-UTC= 17.0 S + (MJD - 41317.) X 0.0 S - * 1979 JAN 1 =JD 2443874.5 TAI-UTC= 18.0 S + (MJD - 41317.) X 0.0 S - * 1980 JAN 1 =JD 2444239.5 TAI-UTC= 19.0 S + (MJD - 41317.) X 0.0 S - * 1981 JUL 1 =JD 2444786.5 TAI-UTC= 20.0 S + (MJD - 41317.) X 0.0 S - * 1982 JUL 1 =JD 2445151.5 TAI-UTC= 21.0 S + (MJD - 41317.) X 0.0 S - * 1983 JUL 1 =JD 2445516.5 TAI-UTC= 22.0 S + (MJD - 41317.) X 0.0 S - * 1985 JUL 1 =JD 2446247.5 TAI-UTC= 23.0 S + (MJD - 41317.) X 0.0 S - * 1988 JAN 1 =JD 2447161.5 TAI-UTC= 24.0 S + (MJD - 41317.) X 0.0 S - * 1990 JAN 1 =JD 2447892.5 TAI-UTC= 25.0 S + (MJD - 41317.) X 0.0 S - * 1991 JAN 1 =JD 2448257.5 TAI-UTC= 26.0 S + (MJD - 41317.) X 0.0 S - * 1992 JUL 1 =JD 2448804.5 TAI-UTC= 27.0 S + (MJD - 41317.) X 0.0 S - * 1993 JUL 1 =JD 2449169.5 TAI-UTC= 28.0 S + (MJD - 41317.) X 0.0 S - * 1994 JUL 1 =JD 2449534.5 TAI-UTC= 29.0 S + (MJD - 41317.) X 0.0 S - * 1996 JAN 1 =JD 2450083.5 TAI-UTC= 30.0 S + (MJD - 41317.) X 0.0 S - * 1997 JUL 1 =JD 2450630.5 TAI-UTC= 31.0 S + (MJD - 41317.) X 0.0 S - * 1999 JAN 1 =JD 2451179.5 TAI-UTC= 32.0 S + (MJD - 41317.) X 0.0 S + * time varies over time due to the earth's rotation and the movement of the continental plates. * * @author Ross Harman, MHPCC * - * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-25 01:37:34 $ + * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-11 00:55:24 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -63,4 +20,5 @@ #include #include +#include #include "psTime.h" @@ -68,84 +26,358 @@ #include "psMemory.h" #include "psAbort.h" - - -/** Number of available leapsecond updates */ -#define NUM_UPDATES 36 - -/** Maximum length of time string */ -#define MAX_TIME_STRING_LENGTH 256 - -/** Seconds per minute */ -#define SEC_PER_MINUTE 60.0 - -/** Seconds per hour */ -#define SEC_PER_HOUR (60.0*SEC_PER_MINUTE) - -/** Seconds per day */ -#define SEC_PER_DAY (24.0*SEC_PER_HOUR) - -/** Seconds per year */ -#define SEC_PER_YEAR (365.0*SEC_PER_DAY) - -/** Microseconds per day */ -#define USEC_PER_DAY 86400000000.0 - -/** Preprocessor macro to check for null psTime struct */ -#define CHECK_NULL_TIME(TIME,RETURN) \ -if(TIME == NULL) { \ - psError(__func__,"NULL value not allowed"); \ - return RETURN; \ -} - -/** Preprocessor macro to allocate psTime strict */ -#define ALLOC_TIME(TIME) \ -TIME = (psTime*)psAlloc(sizeof(psTime)); \ -if (TIME == NULL) { \ - psAbort(__func__, " : Line %d - Failed to allocate memory", __LINE__); \ -} - -/* - The table below comes from ftp://maia.usno.navy.mil/ser7/tai-utc.dat, which reports the current state of - the TAI/UTC leapsecond conversion times. -*/ -static double lsTable[NUM_UPDATES][4] = { - {2437300.5, 1.422818, 37300.0, 0.0012960}, - {2437512.5, 1.372818, 37300.0, 0.0012960}, - {2437665.5, 1.845858, 37665.0, 0.0011232}, - {2438334.5, 1.945858, 37665.0, 0.0011232}, - {2438395.5, 3.240130, 38761.0, 0.0012960}, - {2438486.5, 3.340130, 38761.0, 0.0012960}, - {2438639.5, 3.440130, 38761.0, 0.0012960}, - {2438761.5, 3.540130, 38761.0, 0.0012960}, - {2438820.5, 3.640130, 38761.0, 0.0012960}, - {2438942.5, 3.740130, 38761.0, 0.0012960}, - {2439004.5, 3.840130, 38761.0, 0.0012960}, - {2439126.5, 4.313170, 39126.0, 0.0025920}, - {2439887.5, 4.213170, 39126.0, 0.0025920}, - {2441317.5, 10.000000, 41317.0, 0.0000000}, - {2441499.5, 11.000000, 41317.0, 0.0000000}, - {2441683.5, 12.000000, 41317.0, 0.0000000}, - {2442048.5, 13.000000, 41317.0, 0.0000000}, - {2442413.5, 14.000000, 41317.0, 0.0000000}, - {2442778.5, 15.000000, 41317.0, 0.0000000}, - {2443144.5, 16.000000, 41317.0, 0.0000000}, - {2443509.5, 17.000000, 41317.0, 0.0000000}, - {2443874.5, 18.000000, 41317.0, 0.0000000}, - {2444239.5, 19.000000, 41317.0, 0.0000000}, - {2444786.5, 20.000000, 41317.0, 0.0000000}, - {2445151.5, 21.000000, 41317.0, 0.0000000}, - {2445516.5, 22.000000, 41317.0, 0.0000000}, - {2446247.5, 23.000000, 41317.0, 0.0000000}, - {2447161.5, 24.000000, 41317.0, 0.0000000}, - {2447892.5, 25.000000, 41317.0, 0.0000000}, - {2448257.5, 26.000000, 41317.0, 0.0000000}, - {2448804.5, 27.000000, 41317.0, 0.0000000}, - {2449169.5, 28.000000, 41317.0, 0.0000000}, - {2449534.5, 29.000000, 41317.0, 0.0000000}, - {2450083.5, 30.000000, 41317.0, 0.0000000}, - {2450630.5, 31.000000, 41317.0, 0.0000000}, - {2451179.5, 32.000000, 41317.0, 0.0000000} - }; - +#include "psImage.h" + + +/** Sidereal angular conversion from seconds to radians for GMST in seconds (i.e. pi/(180*240)) */ +#define S2R (7.272205216643039903848711535369e-5) + +/** Two times pi with double precision accuracy */ +#define TWOPI (2.0*M_PIl) + +/** Conversion from radians to degrees */ +#define R2DEG = (180.0/M_PIl) + + /** Maximum length of time string */ + #define MAX_TIME_STRING_LENGTH 256 + + /** Seconds per minute */ + #define SEC_PER_MINUTE 60.0 + + /** Seconds per hour */ + #define SEC_PER_HOUR (60.0*SEC_PER_MINUTE) + + /** Seconds per day */ + #define SEC_PER_DAY (24.0*SEC_PER_HOUR) + + /** Seconds per year */ + #define SEC_PER_YEAR (365.0*SEC_PER_DAY) + + /** Microseconds per day */ + #define USEC_PER_DAY 86400000000.0 + + /** Max size of a line read */ + #define LINESIZE 1024 + + /** Preprocessor macro to check for null psTime struct */ + #define CHECK_NULL_TIME(TIME,RETURN) \ + if(TIME == NULL) { \ + psError(__func__,"NULL value not allowed"); \ + return RETURN; \ + } + + /** Preprocessor macro to allocate psTime strict */ + #define ALLOC_TIME(TIME) \ + TIME = (psTime*)psAlloc(sizeof(psTime)); \ +if (TIME == NULL) { \ + psAbort(__func__, " : Line %d - Failed to allocate memory", __LINE__); \ +} \ +TIME->sec = 0; \ +TIME->usec = 0; + + +psImage* readTaiUtcFile(char *fileName) +{ + char line[LINESIZE]; + int j = 0; + int maxLines = 100; + psF64 *ptr = NULL; + psImage *table = NULL; + FILE *fd = NULL; + + + fd=fopen(fileName, "r"); + if(fd==NULL) { + psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName); + 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(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, 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 + } + + *(unsigned int *)&table->numRows = j; + + return table; +} + +psImage* readSer7File(char *fileName) +{ + bool beginRecord = false; + char line[LINESIZE]; + int j = 0; + int maxLines = 400; + psImage *table = NULL; + FILE *fd = NULL; + + + fd=fopen(fileName, "r"); + if(fd==NULL) { + psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName); + 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(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, 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; + } + } + *(unsigned int *)&table->numRows = j; + + return table; +} + +psImage* readEopcFile(char *fileName) +{ + char line[LINESIZE]; + int i = 0; + int j = 0; + int maxLines = 2500; + psImage *table = NULL; + FILE *fd = NULL; + + + fd=fopen(fileName, "r"); + if(fd==NULL) { + psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName); + return NULL ; + } + + // Table shouldn't be larger than 2500 rows. All columns are read. + table = psImageAlloc(11, maxLines, PS_TYPE_F64); + for(i=0; i<3; i++) { + if((fgets(line, LINESIZE, fd) != NULL)) { + psError(__func__, " : Line %d - Error reading %s", __LINE__, fileName); + } + } + + while(fgets(line, LINESIZE, fd) != NULL) { + if(j>=maxLines) { + psError(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, maxLines); + return NULL; + psFree(table); + } + psF64 *ptr = table->data.F64[j++]; + sscanf(line, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", + ptr, ptr+1, ptr+2, ptr+3, ptr+4, ptr+5, ptr+6, ptr+7, ptr+8, ptr+9, ptr+10); + } + + *(unsigned int *)&table->numRows = j; + + return table; +} + +psImage* readFinalsFile(char *fileName) +{ + char line[LINESIZE]; + int j = 0; + int maxLines = 15000; + psImage *table = NULL; + FILE *fd = NULL; + + + fd=fopen(fileName, "r"); + if(fd==NULL) { + psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName); + return NULL ; + } + + // Table shouldn't be larger than 15000 rows. Select columns are read. + table = psImageAlloc(6, maxLines, PS_TYPE_F64); + while(fgets(line, LINESIZE, fd) != NULL) { + if(j>=maxLines) { + psError(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, maxLines); + return NULL; + psFree(table); + } + + psF64 *ptr = table->data.F64[j++]; + ptr = table->data.F64[j++]; + + // Sometimes there are blank entries in this table, so check for a space to determine if there is no entry + if(*(line+7) == ' ') { // MJD + *ptr = 0.0; + } else { + *ptr = atof(line+7); + } + + if(*(line+18) == ' ') { // Polar motion X + *(ptr+1) = 0.0; + } else { + *(ptr+1) = atof(line+18); + } + + if(*(line+27) == ' ') { // Polar motion X error + *(ptr+2) = 0.0; + } else { + *(ptr+2) = atof(line+27); + } + + if(*(line+37) == ' ') { // Polar motion Y + *(ptr+3) = 0.0; + } else { + *(ptr+3) = atof(line+37); + } + + if(*(line+46) == ' ') { // Polar motion Y error + *(ptr+4) = 0.0; + } else { + *(ptr+4) = atof(line+46); + } + + if(*(line+58) == ' ') { // UT1-UTC + *(ptr+5) = 0.0; + } else { + *(ptr+5) = atof(line+58); + } + } + + *(unsigned int *)&table->numRows = j; + + return table; +} + + +double lookupTaiUtcTable(psImage *table, psTime *time) +{ + int hiIdx = 0; + int loIdx = 0; + int 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; + + + // 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(__func__, " : Line %d - High index too big. Value: %d", __LINE__, hiIdx); + return 0.0; + } + } + loIdx = hiIdx--; + if(loIdx < 0) { + psError(__func__, " : Line %d - Low index too small. Value: %d", __LINE__, loIdx); + return 0.0; + } + + denom = table->data.F64[hiIdx][0]-table->data.F64[loIdx][0]; + if(fabs(denom) < FLT_EPSILON) { + psError(__func__, " : Line %d - Divide by zero error during interpolation", __LINE__); + return 0.0; + } else { + + // Interpolate three constants required to calculate TAI-UTC + const1 = table->data.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; +} + +double lookupSer7Table(psImage *table, psTime *time, int col) +{ + int hiIdx = 0; + int loIdx = 0; + int numRows = 0; + double out = 0.0; + double denom = 0.0; + double mjdUtc = 0.0; + + + // 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(__func__, " : Line %d - MJD too small. Value: %d", __LINE__, mjdUtc); + + } else if(mjdUtc > table->data.F64[numRows-1][3]) { // MJD time after end of table + psError(__func__, " : Line %d - MJD too big. Value: %d", __LINE__, mjdUtc); + } else { // All other times..interpolate + while(mjdUtc > table->data.F64[hiIdx][3]) { + hiIdx++; + if(hiIdx >= numRows) { + psError(__func__, " : Line %d - High index too big. Value: %d", __LINE__, hiIdx); + return 0.0; + } + } + + loIdx = hiIdx--; + if(loIdx < 0) { + psError(__func__, " : Line %d - Low index too small. Value: %d", __LINE__, loIdx); + } + + denom = table->data.F64[hiIdx][3] - table->data.F64[loIdx][3]; + if(fabs(denom) < FLT_EPSILON) { + psError(__func__, " : Line %d - Dividie by zero error during interpolation", __LINE__); + } 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; +} psTime* psTimeGetTime(psTimeType type) @@ -153,4 +385,5 @@ struct timeval now; psTime *time = NULL; + // Allocate psTime struct @@ -178,4 +411,5 @@ { double delta = 0.0; + if(time->type == type) { @@ -199,58 +433,128 @@ double psTimeToLST(psTime *time, double longitude) { - double delta = 0; + psF64 jdTdtDays = 0.0; + psF64 jdUt1Days = 0.0; + psF64 mjdUt1Days = 0.0; + psF64 lmstRad = 0.0; + psF64 fracDays = 0.0; + psF64 gmstRad = 0.0; + psF64 t = 0.0; + psF64 tu = 0.0; + psF64 const1 = 24110.5493771; + psF64 const2 = 8639877.3173760; + psF64 const3 = 307.4771600; + psF64 const4 = 0.0931118; + psF64 const5 = -0.0000062; + psF64 const6 = 0.0000013; + psTime *ut1UtcDelta = NULL; + psTime *tdtDelta = NULL; + psTime *tdtTime = NULL; + psTime *taiTime = NULL; + psTime *utcTime = NULL; + psTime *ut1Time = NULL; + + + // Calculate TAI or UTC time based on type of time user passes + if(time->type == PS_TIME_TAI) { + taiTime = psMemIncrRefCounter(time); + ALLOC_TIME(utcTime); + utcTime->sec = taiTime->sec - psGetTAIDelta(time); + } else if(time->type == PS_TIME_UTC) { + utcTime = psMemIncrRefCounter(time); + ALLOC_TIME(taiTime); + taiTime->sec = utcTime->sec + psGetTAIDelta(time); + } + + // Convert Universal Time (UTC) to UT1 + ALLOC_TIME(ut1UtcDelta); + ALLOC_TIME(ut1Time); + ut1UtcDelta->usec = psGetUT1Delta(utcTime)*1e6; + ut1Time = psTimeAdd(utcTime, ut1UtcDelta); + + // Calculate UT1 as Julian Centuries since J2000.0 + jdUt1Days = psTimeToJD(ut1Time); + mjdUt1Days = psTimeToMJD(ut1Time); + t = (jdUt1Days - 2451545.0)/36525.0; + + // Calculate Terrestial Dynamical Time (TDT) + ALLOC_TIME(tdtTime); + ALLOC_TIME(tdtDelta); + tdtDelta->sec = 32; + tdtDelta->usec = 184000; + tdtTime = psTimeAdd(taiTime, tdtDelta); + + // Calculate TDT as Julian centuries since J2000.0 + jdTdtDays = psTimeToJD(tdtTime); + tu = (jdTdtDays - 2451545.0)/36525.0; + + // Calculate fractional part of MJD + 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; + + // Place GMST between 0 and 2*pi + gmstRad = fmod(gmstRad, TWOPI); + + // Calculate Local Mean Sidereal Time (LMST) in radians + lmstRad = gmstRad + longitude; + + // Free temporary structs + psFree(ut1UtcDelta); + psFree(ut1Time); + psFree(tdtTime); + psFree(tdtDelta); + psFree(utcTime); + psFree(taiTime); + + return lmstRad; +} + +double psGetUT1Delta(psTime *time) +{ + psImage *iersTable = NULL; + double ut1UtcDeltaUsec = 0.0; + + + iersTable = readSer7File("../../data/ser7.dat"); + ut1UtcDeltaUsec = lookupSer7Table(iersTable, time, 6)*1.0e6; + + return ut1UtcDeltaUsec; +} + +double psGetTAIDelta(psTime *time) +{ + double delta = 0.0; + psImage *taiUtcTable = NULL; + + + taiUtcTable = readTaiUtcFile("../../data/tai-utc.dat"); + delta = lookupTaiUtcTable(taiUtcTable, time); return delta; } -psTime* psLSTToTime(double time, double longitude) -{ - psTime *outTime = NULL; - - return outTime; -} - -double psGetUT1Delta(psTime *time) -{ - double delta = 0; - - return delta; -} - -double psGetTAIDelta(psTime *time) -{ - int i = 0; - double delta = 0.0; +psF64 psTimeLeapseconds(const psTime *time1, const psTime *time2) +{ + psF64 diff = 0.0; + + // NULL error checks + CHECK_NULL_TIME(time1,0); + CHECK_NULL_TIME(time2,0); + + diff = fabs(psGetTAIDelta((psTime*)time1)-psGetTAIDelta((psTime*)time2)); + + return diff; +} + +double psTimeToJD(psTime *time) +{ double jd = 0.0; - double mjd = 0.0; - - // Determine Julian and modified Julian dates used in table lookup and time delta calculation - jd = psTimeToJD(time); - mjd = psTimeToMJD(time); - - // Perform table lookup - if(jd < lsTable[0][1]) { // Times earlier than 1961 - delta = 0.0; - } else if(jd > lsTable[NUM_UPDATES-1][1]) { // Times greater than latest table entry - delta = lsTable[NUM_UPDATES-1][1] + (mjd - lsTable[NUM_UPDATES-1][2]) * lsTable[NUM_UPDATES-1][3]; - } else { // All other times - for (i = 0; ilsTable[i][0] && jdsec < 0) { jd = time->sec / SEC_PER_DAY - time->usec / USEC_PER_DAY + 2440587.5; // psTime earlier than epoch @@ -266,8 +570,9 @@ double mjd = 0.0; + // NULL error check CHECK_NULL_TIME(time,0.0); - // Modified Julian date conversion courtesy of Eugene Magnier + // Modified Julian date conversion if(time->sec < 0) { mjd = time->sec / SEC_PER_DAY - time->usec / USEC_PER_DAY + 40587.0; // psTime earlier than epoch @@ -287,4 +592,5 @@ time_t sec; + // NULL error check CHECK_NULL_TIME(time,NULL); @@ -316,4 +622,5 @@ struct timeval timevalTime; + // NULL error check CHECK_NULL_TIME(time,timevalTime); @@ -327,9 +634,61 @@ struct tm* psTimeToTM(psTime *time) { - struct tm *tmTime = NULL; - time_t sec; - - sec = time->sec; - tmTime = gmtime(&sec); + long cent = 0; + long year = 0; + long month = 0; + long day = 0; + long hour = 0; + long minute = 0; + long seconds = 0; + long temp = 0; + struct tm* tmTime = 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 + if(day > (pow(2, 28)-307)) + { + temp = (day - 146097+306)/146097+1; // Avoid overflow if day close to maxint + day -= temp * 146097-306; + } else if((day += 306 ) <= 0) + { + temp = -( -day / 146097 + 1); // Avoid ambiguity in C division of negatives + 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 + year = (day*4-1)/1461; // Calc number of years into the century + day -= year*1461/4; // Again March-based (4 yrs =\u02dc 146[01] days) + month = (day*12+1093)/367; // Get the month (3..14 represent March through + 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) + { + year++; + month -= 12; + } + + // Allocate output + tmTime = (struct tm*)psAlloc(sizeof(struct tm)); + if (tmTime == NULL) + { + psAbort(__func__, " : Line %d - Failed to allocate memory", __LINE__); + } + + tmTime->tm_year = year - 1900; + tmTime->tm_mon = month - 1; + tmTime->tm_mday = day + 1; + tmTime->tm_hour = hour; + tmTime->tm_min = minute; + tmTime->tm_sec = seconds; + tmTime->tm_isdst = -1; return tmTime; @@ -341,4 +700,5 @@ double seconds = 0.0; psTime *outTime = NULL; + // Allocate psTime struct @@ -363,4 +723,5 @@ double seconds = 0.0; psTime *outTime = NULL; + // Allocate psTime struct @@ -393,4 +754,5 @@ psTime *outTime = NULL; + // Preserve input string by creating duplicate strncpy(tempString, time, MAX_TIME_STRING_LENGTH); @@ -461,9 +823,10 @@ psTime *outTime = NULL; + // NULL error check CHECK_NULL_TIME(time,NULL); // Allocate psTime struct - ALLOC_TIME(outTime); + ALLOC_TIME(outTime) // Convert to psTime @@ -476,7 +839,11 @@ psTime* psTMToTime(struct tm* time) { - time_t localTimeInSeconds = 0; - time_t temp = 0; - struct tm *localTMTime = NULL; + long year; + long month; + long day; + long hour; + long minute; + long seconds; + long temp; psTime *outTime = NULL; @@ -487,26 +854,49 @@ ALLOC_TIME(outTime); - // The mktime function returns seconds since last epoch in local time, so get local time and offset - localTimeInSeconds = mktime(time); - localTMTime = localtime(&temp); - - // Convert to psTime + // Extract data from TM struct + year = time->tm_year + 1900; + month = time->tm_mon + 1; + day = time->tm_mday; + hour = time->tm_hour; + minute = time->tm_min; + seconds = time->tm_sec; + + // Make month in range 3..14 (treat Jan & Feb as months 13..14 of prev year) + if( month <= 2 ) + { + year -= (temp = (14 - month) / 12); + month += 12 * temp; + } else if(month > 14) + { + year += (temp = (month - 3) / 12); + month -= 12 * temp; + } + + // Make year positive + if (year < 0 ) + { + day -= 146097 * (temp = (399 - year) / 400); + year += 400 * temp; + } + + // Add day of month, days of previous 0-11 month period that began w/March, days of previous 0-399 year + // period that began w/March of a 400-multiple year), days of any 400-year periods before that, and 306 + // 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; + + // C's TM does not define a microsecond field. Microseconds must be manipulated by calling function. outTime->usec = 0; - outTime->sec = localTMTime->tm_gmtoff + localTimeInSeconds; return outTime; } -psTime* psTAIAdd(psTime *tai1, psTime *tai2) +psTime* psTimeAdd(psTime *tai1, psTime *tai2) { psTime *outTime = NULL; - if(tai1->type != PS_TIME_TAI) { - psError(__func__, " : Line %d - First argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai1->type); - return NULL; - } - - if(tai2->type != PS_TIME_TAI) { - psError(__func__, " : Line %d - Second argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai2->type); + + if(tai1->type != tai2->type) { + psError(__func__, " : Line %d - Argument types inconsistent for time addition", __LINE__); return NULL; } @@ -525,15 +915,11 @@ } -psTime* psTAISub(psTime *tai1, psTime *tai2) +psTime* psTimeSub(psTime *tai1, psTime *tai2) { psTime *outTime = NULL; - if(tai1->type != PS_TIME_TAI) { - psError(__func__, " : Line %d - First argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai1->type); - return NULL; - } - - if(tai2->type != PS_TIME_TAI) { - psError(__func__, " : Line %d - Second argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai2->type); + + if(tai1->type != tai2->type) { + psError(__func__, " : Line %d - Argument types inconsistent for time subtraction", __LINE__); return NULL; } @@ -552,15 +938,11 @@ } -psTime* psTAIDelta(psTime *tai1, psTime *tai2) +psTime* psTimeDelta(psTime *tai1, psTime *tai2) { psTime *outTime = NULL; - if(tai1->type != PS_TIME_TAI) { - psError(__func__, " : Line %d - First argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai1->type); - return NULL; - } - - if(tai2->type != PS_TIME_TAI) { - psError(__func__, " : Line %d - Second argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai2->type); + + if(tai1->type != tai2->type) { + psError(__func__, " : Line %d - Argument types inconsistent for time delta", __LINE__); return NULL; } @@ -578,5 +960,4 @@ outTime->sec = abs(outTime->sec); - return outTime; }