Changeset 1788 for trunk/psLib/src/astro/psTime.c
- Timestamp:
- Sep 10, 2004, 2:57:31 PM (22 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/astro/psTime.c (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/astro/psTime.c
r1614 r1788 7 7 * of conversions between specific time formats. PSLib currently uses the UNIX timeval time system as the 8 8 * base upon which International Atomic Time (TAI) and Universal Time Coordinated (UTC) are calculated. TAI 9 * time varies over time due to the earth's rotation and the movement of the continental plates. It is 10 * currentls 32 seconds faster than UTC/timeval time. The conversion between TAI and UTC time for other dates 11 * is given in the table below, which comes from: 12 * 13 * ftp://maia.usno.navy.mil/ser7/tai-utc.dat 14 * 15 * This table is placed into a static array and used by the psTime functions for conversions. 16 * 17 * 1961 JAN 1 =JD 2437300.5 TAI-UTC= 1.4228180 S + (MJD - 37300.) X 0.001296 S 18 * 1961 AUG 1 =JD 2437512.5 TAI-UTC= 1.3728180 S + (MJD - 37300.) X 0.001296 S 19 * 1962 JAN 1 =JD 2437665.5 TAI-UTC= 1.8458580 S + (MJD - 37665.) X 0.0011232S 20 * 1963 NOV 1 =JD 2438334.5 TAI-UTC= 1.9458580 S + (MJD - 37665.) X 0.0011232S 21 * 1964 JAN 1 =JD 2438395.5 TAI-UTC= 3.2401300 S + (MJD - 38761.) X 0.001296 S 22 * 1964 APR 1 =JD 2438486.5 TAI-UTC= 3.3401300 S + (MJD - 38761.) X 0.001296 S 23 * 1964 SEP 1 =JD 2438639.5 TAI-UTC= 3.4401300 S + (MJD - 38761.) X 0.001296 S 24 * 1965 JAN 1 =JD 2438761.5 TAI-UTC= 3.5401300 S + (MJD - 38761.) X 0.001296 S 25 * 1965 MAR 1 =JD 2438820.5 TAI-UTC= 3.6401300 S + (MJD - 38761.) X 0.001296 S 26 * 1965 JUL 1 =JD 2438942.5 TAI-UTC= 3.7401300 S + (MJD - 38761.) X 0.001296 S 27 * 1965 SEP 1 =JD 2439004.5 TAI-UTC= 3.8401300 S + (MJD - 38761.) X 0.001296 S 28 * 1966 JAN 1 =JD 2439126.5 TAI-UTC= 4.3131700 S + (MJD - 39126.) X 0.002592 S 29 * 1968 FEB 1 =JD 2439887.5 TAI-UTC= 4.2131700 S + (MJD - 39126.) X 0.002592 S 30 * 1972 JAN 1 =JD 2441317.5 TAI-UTC= 10.0 S + (MJD - 41317.) X 0.0 S 31 * 1972 JUL 1 =JD 2441499.5 TAI-UTC= 11.0 S + (MJD - 41317.) X 0.0 S 32 * 1973 JAN 1 =JD 2441683.5 TAI-UTC= 12.0 S + (MJD - 41317.) X 0.0 S 33 * 1974 JAN 1 =JD 2442048.5 TAI-UTC= 13.0 S + (MJD - 41317.) X 0.0 S 34 * 1975 JAN 1 =JD 2442413.5 TAI-UTC= 14.0 S + (MJD - 41317.) X 0.0 S 35 * 1976 JAN 1 =JD 2442778.5 TAI-UTC= 15.0 S + (MJD - 41317.) X 0.0 S 36 * 1977 JAN 1 =JD 2443144.5 TAI-UTC= 16.0 S + (MJD - 41317.) X 0.0 S 37 * 1978 JAN 1 =JD 2443509.5 TAI-UTC= 17.0 S + (MJD - 41317.) X 0.0 S 38 * 1979 JAN 1 =JD 2443874.5 TAI-UTC= 18.0 S + (MJD - 41317.) X 0.0 S 39 * 1980 JAN 1 =JD 2444239.5 TAI-UTC= 19.0 S + (MJD - 41317.) X 0.0 S 40 * 1981 JUL 1 =JD 2444786.5 TAI-UTC= 20.0 S + (MJD - 41317.) X 0.0 S 41 * 1982 JUL 1 =JD 2445151.5 TAI-UTC= 21.0 S + (MJD - 41317.) X 0.0 S 42 * 1983 JUL 1 =JD 2445516.5 TAI-UTC= 22.0 S + (MJD - 41317.) X 0.0 S 43 * 1985 JUL 1 =JD 2446247.5 TAI-UTC= 23.0 S + (MJD - 41317.) X 0.0 S 44 * 1988 JAN 1 =JD 2447161.5 TAI-UTC= 24.0 S + (MJD - 41317.) X 0.0 S 45 * 1990 JAN 1 =JD 2447892.5 TAI-UTC= 25.0 S + (MJD - 41317.) X 0.0 S 46 * 1991 JAN 1 =JD 2448257.5 TAI-UTC= 26.0 S + (MJD - 41317.) X 0.0 S 47 * 1992 JUL 1 =JD 2448804.5 TAI-UTC= 27.0 S + (MJD - 41317.) X 0.0 S 48 * 1993 JUL 1 =JD 2449169.5 TAI-UTC= 28.0 S + (MJD - 41317.) X 0.0 S 49 * 1994 JUL 1 =JD 2449534.5 TAI-UTC= 29.0 S + (MJD - 41317.) X 0.0 S 50 * 1996 JAN 1 =JD 2450083.5 TAI-UTC= 30.0 S + (MJD - 41317.) X 0.0 S 51 * 1997 JUL 1 =JD 2450630.5 TAI-UTC= 31.0 S + (MJD - 41317.) X 0.0 S 52 * 1999 JAN 1 =JD 2451179.5 TAI-UTC= 32.0 S + (MJD - 41317.) X 0.0 S 9 * time varies over time due to the earth's rotation and the movement of the continental plates. 53 10 * 54 11 * @author Ross Harman, MHPCC 55 12 * 56 * @version $Revision: 1.2 1$ $Name: not supported by cvs2svn $57 * @date $Date: 2004-0 8-25 01:37:34 $13 * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ 14 * @date $Date: 2004-09-11 00:55:24 $ 58 15 * 59 16 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 63 20 #include <stdlib.h> 64 21 #include <string.h> 22 #include <math.h> 65 23 66 24 #include "psTime.h" … … 68 26 #include "psMemory.h" 69 27 #include "psAbort.h" 70 71 72 /** Number of available leapsecond updates */ 73 #define NUM_UPDATES 36 74 75 /** Maximum length of time string */ 76 #define MAX_TIME_STRING_LENGTH 256 77 78 /** Seconds per minute */ 79 #define SEC_PER_MINUTE 60.0 80 81 /** Seconds per hour */ 82 #define SEC_PER_HOUR (60.0*SEC_PER_MINUTE) 83 84 /** Seconds per day */ 85 #define SEC_PER_DAY (24.0*SEC_PER_HOUR) 86 87 /** Seconds per year */ 88 #define SEC_PER_YEAR (365.0*SEC_PER_DAY) 89 90 /** Microseconds per day */ 91 #define USEC_PER_DAY 86400000000.0 92 93 /** Preprocessor macro to check for null psTime struct */ 94 #define CHECK_NULL_TIME(TIME,RETURN) \ 95 if(TIME == NULL) { \ 96 psError(__func__,"NULL value not allowed"); \ 97 return RETURN; \ 98 } 99 100 /** Preprocessor macro to allocate psTime strict */ 101 #define ALLOC_TIME(TIME) \ 102 TIME = (psTime*)psAlloc(sizeof(psTime)); \ 103 if (TIME == NULL) { \ 104 psAbort(__func__, " : Line %d - Failed to allocate memory", __LINE__); \ 105 } 106 107 /* 108 The table below comes from ftp://maia.usno.navy.mil/ser7/tai-utc.dat, which reports the current state of 109 the TAI/UTC leapsecond conversion times. 110 */ 111 static double lsTable[NUM_UPDATES][4] = { 112 {2437300.5, 1.422818, 37300.0, 0.0012960}, 113 {2437512.5, 1.372818, 37300.0, 0.0012960}, 114 {2437665.5, 1.845858, 37665.0, 0.0011232}, 115 {2438334.5, 1.945858, 37665.0, 0.0011232}, 116 {2438395.5, 3.240130, 38761.0, 0.0012960}, 117 {2438486.5, 3.340130, 38761.0, 0.0012960}, 118 {2438639.5, 3.440130, 38761.0, 0.0012960}, 119 {2438761.5, 3.540130, 38761.0, 0.0012960}, 120 {2438820.5, 3.640130, 38761.0, 0.0012960}, 121 {2438942.5, 3.740130, 38761.0, 0.0012960}, 122 {2439004.5, 3.840130, 38761.0, 0.0012960}, 123 {2439126.5, 4.313170, 39126.0, 0.0025920}, 124 {2439887.5, 4.213170, 39126.0, 0.0025920}, 125 {2441317.5, 10.000000, 41317.0, 0.0000000}, 126 {2441499.5, 11.000000, 41317.0, 0.0000000}, 127 {2441683.5, 12.000000, 41317.0, 0.0000000}, 128 {2442048.5, 13.000000, 41317.0, 0.0000000}, 129 {2442413.5, 14.000000, 41317.0, 0.0000000}, 130 {2442778.5, 15.000000, 41317.0, 0.0000000}, 131 {2443144.5, 16.000000, 41317.0, 0.0000000}, 132 {2443509.5, 17.000000, 41317.0, 0.0000000}, 133 {2443874.5, 18.000000, 41317.0, 0.0000000}, 134 {2444239.5, 19.000000, 41317.0, 0.0000000}, 135 {2444786.5, 20.000000, 41317.0, 0.0000000}, 136 {2445151.5, 21.000000, 41317.0, 0.0000000}, 137 {2445516.5, 22.000000, 41317.0, 0.0000000}, 138 {2446247.5, 23.000000, 41317.0, 0.0000000}, 139 {2447161.5, 24.000000, 41317.0, 0.0000000}, 140 {2447892.5, 25.000000, 41317.0, 0.0000000}, 141 {2448257.5, 26.000000, 41317.0, 0.0000000}, 142 {2448804.5, 27.000000, 41317.0, 0.0000000}, 143 {2449169.5, 28.000000, 41317.0, 0.0000000}, 144 {2449534.5, 29.000000, 41317.0, 0.0000000}, 145 {2450083.5, 30.000000, 41317.0, 0.0000000}, 146 {2450630.5, 31.000000, 41317.0, 0.0000000}, 147 {2451179.5, 32.000000, 41317.0, 0.0000000} 148 }; 149 28 #include "psImage.h" 29 30 31 /** Sidereal angular conversion from seconds to radians for GMST in seconds (i.e. pi/(180*240)) */ 32 #define S2R (7.272205216643039903848711535369e-5) 33 34 /** Two times pi with double precision accuracy */ 35 #define TWOPI (2.0*M_PIl) 36 37 /** Conversion from radians to degrees */ 38 #define R2DEG = (180.0/M_PIl) 39 40 /** Maximum length of time string */ 41 #define MAX_TIME_STRING_LENGTH 256 42 43 /** Seconds per minute */ 44 #define SEC_PER_MINUTE 60.0 45 46 /** Seconds per hour */ 47 #define SEC_PER_HOUR (60.0*SEC_PER_MINUTE) 48 49 /** Seconds per day */ 50 #define SEC_PER_DAY (24.0*SEC_PER_HOUR) 51 52 /** Seconds per year */ 53 #define SEC_PER_YEAR (365.0*SEC_PER_DAY) 54 55 /** Microseconds per day */ 56 #define USEC_PER_DAY 86400000000.0 57 58 /** Max size of a line read */ 59 #define LINESIZE 1024 60 61 /** Preprocessor macro to check for null psTime struct */ 62 #define CHECK_NULL_TIME(TIME,RETURN) \ 63 if(TIME == NULL) { \ 64 psError(__func__,"NULL value not allowed"); \ 65 return RETURN; \ 66 } 67 68 /** Preprocessor macro to allocate psTime strict */ 69 #define ALLOC_TIME(TIME) \ 70 TIME = (psTime*)psAlloc(sizeof(psTime)); \ 71 if (TIME == NULL) { \ 72 psAbort(__func__, " : Line %d - Failed to allocate memory", __LINE__); \ 73 } \ 74 TIME->sec = 0; \ 75 TIME->usec = 0; 76 77 78 psImage* readTaiUtcFile(char *fileName) 79 { 80 char line[LINESIZE]; 81 int j = 0; 82 int maxLines = 100; 83 psF64 *ptr = NULL; 84 psImage *table = NULL; 85 FILE *fd = NULL; 86 87 88 fd=fopen(fileName, "r"); 89 if(fd==NULL) { 90 psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName); 91 return NULL; 92 } 93 94 // Table shouldn't be larger than 100 rows. Select columns are read based on indices. 95 table = psImageAlloc(4, maxLines, PS_TYPE_F64); 96 97 while(fgets(line, LINESIZE, fd) != NULL) { 98 if(j>=maxLines) { 99 psError(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, maxLines); 100 return NULL; 101 psFree(table); 102 } 103 104 ptr = table->data.F64[j++]; 105 *ptr = atof(line+16); // MJD 106 *(ptr+1) = atof(line+38); // TAI-UTC 107 *(ptr+2) = atof(line+59); // Contant #1 108 *(ptr+3) = atof(line+69); // Contant #2 109 } 110 111 *(unsigned int *)&table->numRows = j; 112 113 return table; 114 } 115 116 psImage* readSer7File(char *fileName) 117 { 118 bool beginRecord = false; 119 char line[LINESIZE]; 120 int j = 0; 121 int maxLines = 400; 122 psImage *table = NULL; 123 FILE *fd = NULL; 124 125 126 fd=fopen(fileName, "r"); 127 if(fd==NULL) { 128 psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName); 129 return NULL ; 130 } 131 132 // Table shouldn't be larger than 400 rows. All columns in area of interest are read. 133 table = psImageAlloc(7, maxLines, PS_TYPE_F64); 134 135 while(fgets(line, LINESIZE, fd) != NULL) { 136 137 // Stop parsing after finding "These" in data file after read has began 138 if(beginRecord) { 139 if(strstr(line, "These")!=NULL) { 140 beginRecord = false; 141 break; 142 } 143 144 if(j>=maxLines) { 145 psError(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, maxLines); 146 return NULL; 147 psFree(table); 148 } 149 psF64 *ptr = table->data.F64[j++]; 150 sscanf(line, "%lf %lf %lf %lf %lf %lf %lf", ptr, ptr+1, ptr+2, ptr+3, ptr+4, ptr+5, ptr+6); 151 } 152 153 // Start parsing after finding "x(arcsec)" in data file 154 if(strstr(line, "x(arcsec)") != NULL) { 155 beginRecord = true; 156 } 157 } 158 *(unsigned int *)&table->numRows = j; 159 160 return table; 161 } 162 163 psImage* readEopcFile(char *fileName) 164 { 165 char line[LINESIZE]; 166 int i = 0; 167 int j = 0; 168 int maxLines = 2500; 169 psImage *table = NULL; 170 FILE *fd = NULL; 171 172 173 fd=fopen(fileName, "r"); 174 if(fd==NULL) { 175 psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName); 176 return NULL ; 177 } 178 179 // Table shouldn't be larger than 2500 rows. All columns are read. 180 table = psImageAlloc(11, maxLines, PS_TYPE_F64); 181 for(i=0; i<3; i++) { 182 if((fgets(line, LINESIZE, fd) != NULL)) { 183 psError(__func__, " : Line %d - Error reading %s", __LINE__, fileName); 184 } 185 } 186 187 while(fgets(line, LINESIZE, fd) != NULL) { 188 if(j>=maxLines) { 189 psError(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, maxLines); 190 return NULL; 191 psFree(table); 192 } 193 psF64 *ptr = table->data.F64[j++]; 194 sscanf(line, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", 195 ptr, ptr+1, ptr+2, ptr+3, ptr+4, ptr+5, ptr+6, ptr+7, ptr+8, ptr+9, ptr+10); 196 } 197 198 *(unsigned int *)&table->numRows = j; 199 200 return table; 201 } 202 203 psImage* readFinalsFile(char *fileName) 204 { 205 char line[LINESIZE]; 206 int j = 0; 207 int maxLines = 15000; 208 psImage *table = NULL; 209 FILE *fd = NULL; 210 211 212 fd=fopen(fileName, "r"); 213 if(fd==NULL) { 214 psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName); 215 return NULL ; 216 } 217 218 // Table shouldn't be larger than 15000 rows. Select columns are read. 219 table = psImageAlloc(6, maxLines, PS_TYPE_F64); 220 while(fgets(line, LINESIZE, fd) != NULL) { 221 if(j>=maxLines) { 222 psError(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, maxLines); 223 return NULL; 224 psFree(table); 225 } 226 227 psF64 *ptr = table->data.F64[j++]; 228 ptr = table->data.F64[j++]; 229 230 // Sometimes there are blank entries in this table, so check for a space to determine if there is no entry 231 if(*(line+7) == ' ') { // MJD 232 *ptr = 0.0; 233 } else { 234 *ptr = atof(line+7); 235 } 236 237 if(*(line+18) == ' ') { // Polar motion X 238 *(ptr+1) = 0.0; 239 } else { 240 *(ptr+1) = atof(line+18); 241 } 242 243 if(*(line+27) == ' ') { // Polar motion X error 244 *(ptr+2) = 0.0; 245 } else { 246 *(ptr+2) = atof(line+27); 247 } 248 249 if(*(line+37) == ' ') { // Polar motion Y 250 *(ptr+3) = 0.0; 251 } else { 252 *(ptr+3) = atof(line+37); 253 } 254 255 if(*(line+46) == ' ') { // Polar motion Y error 256 *(ptr+4) = 0.0; 257 } else { 258 *(ptr+4) = atof(line+46); 259 } 260 261 if(*(line+58) == ' ') { // UT1-UTC 262 *(ptr+5) = 0.0; 263 } else { 264 *(ptr+5) = atof(line+58); 265 } 266 } 267 268 *(unsigned int *)&table->numRows = j; 269 270 return table; 271 } 272 273 274 double lookupTaiUtcTable(psImage *table, psTime *time) 275 { 276 int hiIdx = 0; 277 int loIdx = 0; 278 int numRows = 0; 279 double jd = 0.0; 280 double mjd = 0.0; 281 double out = 0.0; 282 double denom = 0.0; 283 double const1 = 0.0; 284 double const2 = 0.0; 285 double const3 = 0.0; 286 287 288 // Number of rows is the number of rows of data read from the data file 289 numRows = table->numRows; 290 291 // Determine Julian and modified Julian dates used in table lookup and time delta calculation 292 jd = psTimeToJD(time); 293 mjd = psTimeToMJD(time); 294 295 // Variable in should be in MJD format 296 if(jd < table->data.F64[0][0]) { // MJD time before start of table 297 return 0.0; 298 } else if(jd > table->data.F64[numRows-1][0]) { // MJD time after end of table 299 return table->data.F64[numRows-1][1]; 300 } else { // All other times..interpolate 301 while(jd > table->data.F64[hiIdx][0]) { 302 hiIdx++; 303 if(hiIdx >= numRows) { 304 psError(__func__, " : Line %d - High index too big. Value: %d", __LINE__, hiIdx); 305 return 0.0; 306 } 307 } 308 loIdx = hiIdx--; 309 if(loIdx < 0) { 310 psError(__func__, " : Line %d - Low index too small. Value: %d", __LINE__, loIdx); 311 return 0.0; 312 } 313 314 denom = table->data.F64[hiIdx][0]-table->data.F64[loIdx][0]; 315 if(fabs(denom) < FLT_EPSILON) { 316 psError(__func__, " : Line %d - Divide by zero error during interpolation", __LINE__); 317 return 0.0; 318 } else { 319 320 // Interpolate three constants required to calculate TAI-UTC 321 const1 = table->data.F64[loIdx][1]+(table->data.F64[hiIdx][1]-table->data.F64[loIdx][1]) 322 *(jd-table->data.F64[loIdx][0])/denom; 323 324 const2 = table->data.F64[loIdx][2]+(table->data.F64[hiIdx][2]-table->data.F64[loIdx][2]) 325 *(jd-table->data.F64[loIdx][0])/denom; 326 327 const3 = table->data.F64[loIdx][3]+(table->data.F64[hiIdx][3]-table->data.F64[loIdx][3]) 328 *(jd-table->data.F64[loIdx][0])/denom; 329 } 330 } 331 332 out = const1 + (mjd - const2) * const3; 333 334 return out; 335 } 336 337 double lookupSer7Table(psImage *table, psTime *time, int col) 338 { 339 int hiIdx = 0; 340 int loIdx = 0; 341 int numRows = 0; 342 double out = 0.0; 343 double denom = 0.0; 344 double mjdUtc = 0.0; 345 346 347 // Number of rows is the number of rows of data read from the data file 348 numRows = table->numRows; 349 350 // Interpolation value in should be in MJD format 351 mjdUtc = psTimeToMJD(time); 352 if(mjdUtc < table->data.F64[0][3]) { // MJD time before start of table 353 psError(__func__, " : Line %d - MJD too small. Value: %d", __LINE__, mjdUtc); 354 355 } else if(mjdUtc > table->data.F64[numRows-1][3]) { // MJD time after end of table 356 psError(__func__, " : Line %d - MJD too big. Value: %d", __LINE__, mjdUtc); 357 } else { // All other times..interpolate 358 while(mjdUtc > table->data.F64[hiIdx][3]) { 359 hiIdx++; 360 if(hiIdx >= numRows) { 361 psError(__func__, " : Line %d - High index too big. Value: %d", __LINE__, hiIdx); 362 return 0.0; 363 } 364 } 365 366 loIdx = hiIdx--; 367 if(loIdx < 0) { 368 psError(__func__, " : Line %d - Low index too small. Value: %d", __LINE__, loIdx); 369 } 370 371 denom = table->data.F64[hiIdx][3] - table->data.F64[loIdx][3]; 372 if(fabs(denom) < FLT_EPSILON) { 373 psError(__func__, " : Line %d - Dividie by zero error during interpolation", __LINE__); 374 } else { 375 out = table->data.F64[loIdx][col]+(table->data.F64[hiIdx][col]-table->data.F64[loIdx][col]) 376 *(mjdUtc-table->data.F64[loIdx][3])/denom; 377 } 378 } 379 380 return out; 381 } 150 382 151 383 psTime* psTimeGetTime(psTimeType type) … … 153 385 struct timeval now; 154 386 psTime *time = NULL; 387 155 388 156 389 // Allocate psTime struct … … 178 411 { 179 412 double delta = 0.0; 413 180 414 181 415 if(time->type == type) { … … 199 433 double psTimeToLST(psTime *time, double longitude) 200 434 { 201 double delta = 0; 435 psF64 jdTdtDays = 0.0; 436 psF64 jdUt1Days = 0.0; 437 psF64 mjdUt1Days = 0.0; 438 psF64 lmstRad = 0.0; 439 psF64 fracDays = 0.0; 440 psF64 gmstRad = 0.0; 441 psF64 t = 0.0; 442 psF64 tu = 0.0; 443 psF64 const1 = 24110.5493771; 444 psF64 const2 = 8639877.3173760; 445 psF64 const3 = 307.4771600; 446 psF64 const4 = 0.0931118; 447 psF64 const5 = -0.0000062; 448 psF64 const6 = 0.0000013; 449 psTime *ut1UtcDelta = NULL; 450 psTime *tdtDelta = NULL; 451 psTime *tdtTime = NULL; 452 psTime *taiTime = NULL; 453 psTime *utcTime = NULL; 454 psTime *ut1Time = NULL; 455 456 457 // Calculate TAI or UTC time based on type of time user passes 458 if(time->type == PS_TIME_TAI) { 459 taiTime = psMemIncrRefCounter(time); 460 ALLOC_TIME(utcTime); 461 utcTime->sec = taiTime->sec - psGetTAIDelta(time); 462 } else if(time->type == PS_TIME_UTC) { 463 utcTime = psMemIncrRefCounter(time); 464 ALLOC_TIME(taiTime); 465 taiTime->sec = utcTime->sec + psGetTAIDelta(time); 466 } 467 468 // Convert Universal Time (UTC) to UT1 469 ALLOC_TIME(ut1UtcDelta); 470 ALLOC_TIME(ut1Time); 471 ut1UtcDelta->usec = psGetUT1Delta(utcTime)*1e6; 472 ut1Time = psTimeAdd(utcTime, ut1UtcDelta); 473 474 // Calculate UT1 as Julian Centuries since J2000.0 475 jdUt1Days = psTimeToJD(ut1Time); 476 mjdUt1Days = psTimeToMJD(ut1Time); 477 t = (jdUt1Days - 2451545.0)/36525.0; 478 479 // Calculate Terrestial Dynamical Time (TDT) 480 ALLOC_TIME(tdtTime); 481 ALLOC_TIME(tdtDelta); 482 tdtDelta->sec = 32; 483 tdtDelta->usec = 184000; 484 tdtTime = psTimeAdd(taiTime, tdtDelta); 485 486 // Calculate TDT as Julian centuries since J2000.0 487 jdTdtDays = psTimeToJD(tdtTime); 488 tu = (jdTdtDays - 2451545.0)/36525.0; 489 490 // Calculate fractional part of MJD 491 fracDays = fmod(mjdUt1Days, 1.0); 492 493 // Calculate Greenwich Mean Sidereal Time (GMST) in radians. Equation set up to minimize multiplications 494 gmstRad = fracDays*TWOPI+(const1+const2*tu+t*(const3+t*(const4+t*(const5+const6*t))))*S2R; 495 496 // Place GMST between 0 and 2*pi 497 gmstRad = fmod(gmstRad, TWOPI); 498 499 // Calculate Local Mean Sidereal Time (LMST) in radians 500 lmstRad = gmstRad + longitude; 501 502 // Free temporary structs 503 psFree(ut1UtcDelta); 504 psFree(ut1Time); 505 psFree(tdtTime); 506 psFree(tdtDelta); 507 psFree(utcTime); 508 psFree(taiTime); 509 510 return lmstRad; 511 } 512 513 double psGetUT1Delta(psTime *time) 514 { 515 psImage *iersTable = NULL; 516 double ut1UtcDeltaUsec = 0.0; 517 518 519 iersTable = readSer7File("../../data/ser7.dat"); 520 ut1UtcDeltaUsec = lookupSer7Table(iersTable, time, 6)*1.0e6; 521 522 return ut1UtcDeltaUsec; 523 } 524 525 double psGetTAIDelta(psTime *time) 526 { 527 double delta = 0.0; 528 psImage *taiUtcTable = NULL; 529 530 531 taiUtcTable = readTaiUtcFile("../../data/tai-utc.dat"); 532 delta = lookupTaiUtcTable(taiUtcTable, time); 202 533 203 534 return delta; 204 535 } 205 536 206 psTime* psLSTToTime(double time, double longitude) 207 { 208 psTime *outTime = NULL; 209 210 return outTime; 211 } 212 213 double psGetUT1Delta(psTime *time) 214 { 215 double delta = 0; 216 217 return delta; 218 } 219 220 double psGetTAIDelta(psTime *time) 221 { 222 int i = 0; 223 double delta = 0.0; 537 psF64 psTimeLeapseconds(const psTime *time1, const psTime *time2) 538 { 539 psF64 diff = 0.0; 540 541 // NULL error checks 542 CHECK_NULL_TIME(time1,0); 543 CHECK_NULL_TIME(time2,0); 544 545 diff = fabs(psGetTAIDelta((psTime*)time1)-psGetTAIDelta((psTime*)time2)); 546 547 return diff; 548 } 549 550 double psTimeToJD(psTime *time) 551 { 224 552 double jd = 0.0; 225 double mjd = 0.0; 226 227 // Determine Julian and modified Julian dates used in table lookup and time delta calculation 228 jd = psTimeToJD(time); 229 mjd = psTimeToMJD(time); 230 231 // Perform table lookup 232 if(jd < lsTable[0][1]) { // Times earlier than 1961 233 delta = 0.0; 234 } else if(jd > lsTable[NUM_UPDATES-1][1]) { // Times greater than latest table entry 235 delta = lsTable[NUM_UPDATES-1][1] + (mjd - lsTable[NUM_UPDATES-1][2]) * lsTable[NUM_UPDATES-1][3]; 236 } else { // All other times 237 for (i = 0; i<NUM_UPDATES-1; i++) { 238 if (jd>lsTable[i][0] && jd<lsTable[i+1][0]) { 239 delta = lsTable[i][1] + (mjd - lsTable[i][2]) * lsTable[i][3]; 240 } 241 } 242 } 243 244 return delta; 245 } 246 247 double psTimeToJD(psTime *time) 248 { 249 double jd = 0.0; 553 250 554 251 555 // NULL error check 252 556 CHECK_NULL_TIME(time,0.0); 253 557 254 // Julian date conversion courtesy of Eugene Magnier558 // Julian date conversion 255 559 if(time->sec < 0) { 256 560 jd = time->sec / SEC_PER_DAY - time->usec / USEC_PER_DAY + 2440587.5; // psTime earlier than epoch … … 266 570 double mjd = 0.0; 267 571 572 268 573 // NULL error check 269 574 CHECK_NULL_TIME(time,0.0); 270 575 271 // Modified Julian date conversion courtesy of Eugene Magnier576 // Modified Julian date conversion 272 577 if(time->sec < 0) { 273 578 mjd = time->sec / SEC_PER_DAY - time->usec / USEC_PER_DAY + 40587.0; // psTime earlier than epoch … … 287 592 time_t sec; 288 593 594 289 595 // NULL error check 290 596 CHECK_NULL_TIME(time,NULL); … … 316 622 struct timeval timevalTime; 317 623 624 318 625 // NULL error check 319 626 CHECK_NULL_TIME(time,timevalTime); … … 327 634 struct tm* psTimeToTM(psTime *time) 328 635 { 329 struct tm *tmTime = NULL; 330 time_t sec; 331 332 sec = time->sec; 333 tmTime = gmtime(&sec); 636 long cent = 0; 637 long year = 0; 638 long month = 0; 639 long day = 0; 640 long hour = 0; 641 long minute = 0; 642 long seconds = 0; 643 long temp = 0; 644 struct tm* tmTime = NULL; 645 646 647 seconds = time->sec%60; 648 minute = time->sec/60%60; 649 hour = time->sec/3600%24; 650 day = (time->sec+62135596800)/86400; 651 652 // Add 306 days to make relative to Mar 1, 0; also adjust day to be within a range (1..2**28-1) where our 653 // calculations will work with 32bit ints 654 if(day > (pow(2, 28)-307)) 655 { 656 temp = (day - 146097+306)/146097+1; // Avoid overflow if day close to maxint 657 day -= temp * 146097-306; 658 } else if((day += 306 ) <= 0) 659 { 660 temp = -( -day / 146097 + 1); // Avoid ambiguity in C division of negatives 661 day -= temp * 146097; 662 } 663 664 cent = (day*4-1)/146097; // Calc number of centuries day is after 29 Feb of yr 0 665 day -= cent*146097/4; // 4 centuries = 146097 days 666 year = (day*4-1)/1461; // Calc number of years into the century 667 day -= year*1461/4; // Again March-based (4 yrs =\u02dc 146[01] days) 668 month = (day*12+1093)/367; // Get the month (3..14 represent March through 669 day -= (month*367-1094)/12; // February of following year) 670 year += cent*100+temp*400; // Get the real year, which is off by 671 672 // One if month is January or February 673 if(month > 12) 674 { 675 year++; 676 month -= 12; 677 } 678 679 // Allocate output 680 tmTime = (struct tm*)psAlloc(sizeof(struct tm)); 681 if (tmTime == NULL) 682 { 683 psAbort(__func__, " : Line %d - Failed to allocate memory", __LINE__); 684 } 685 686 tmTime->tm_year = year - 1900; 687 tmTime->tm_mon = month - 1; 688 tmTime->tm_mday = day + 1; 689 tmTime->tm_hour = hour; 690 tmTime->tm_min = minute; 691 tmTime->tm_sec = seconds; 692 tmTime->tm_isdst = -1; 334 693 335 694 return tmTime; … … 341 700 double seconds = 0.0; 342 701 psTime *outTime = NULL; 702 343 703 344 704 // Allocate psTime struct … … 363 723 double seconds = 0.0; 364 724 psTime *outTime = NULL; 725 365 726 366 727 // Allocate psTime struct … … 393 754 psTime *outTime = NULL; 394 755 756 395 757 // Preserve input string by creating duplicate 396 758 strncpy(tempString, time, MAX_TIME_STRING_LENGTH); … … 461 823 psTime *outTime = NULL; 462 824 825 463 826 // NULL error check 464 827 CHECK_NULL_TIME(time,NULL); 465 828 466 829 // Allocate psTime struct 467 ALLOC_TIME(outTime) ;830 ALLOC_TIME(outTime) 468 831 469 832 // Convert to psTime … … 476 839 psTime* psTMToTime(struct tm* time) 477 840 { 478 time_t localTimeInSeconds = 0; 479 time_t temp = 0; 480 struct tm *localTMTime = NULL; 841 long year; 842 long month; 843 long day; 844 long hour; 845 long minute; 846 long seconds; 847 long temp; 481 848 psTime *outTime = NULL; 482 849 … … 487 854 ALLOC_TIME(outTime); 488 855 489 // The mktime function returns seconds since last epoch in local time, so get local time and offset 490 localTimeInSeconds = mktime(time); 491 localTMTime = localtime(&temp); 492 493 // Convert to psTime 856 // Extract data from TM struct 857 year = time->tm_year + 1900; 858 month = time->tm_mon + 1; 859 day = time->tm_mday; 860 hour = time->tm_hour; 861 minute = time->tm_min; 862 seconds = time->tm_sec; 863 864 // Make month in range 3..14 (treat Jan & Feb as months 13..14 of prev year) 865 if( month <= 2 ) 866 { 867 year -= (temp = (14 - month) / 12); 868 month += 12 * temp; 869 } else if(month > 14) 870 { 871 year += (temp = (month - 3) / 12); 872 month -= 12 * temp; 873 } 874 875 // Make year positive 876 if (year < 0 ) 877 { 878 day -= 146097 * (temp = (399 - year) / 400); 879 year += 400 * temp; 880 } 881 882 // Add day of month, days of previous 0-11 month period that began w/March, days of previous 0-399 year 883 // period that began w/March of a 400-multiple year), days of any 400-year periods before that, and 306 884 // days to adjust from Mar 1, year 0-relative to Jan 1, year 1-relative. Add hours, minutes, and seconds. 885 day += (month * 367 - 1094) / 12 + year % 100 * 1461 / 4 + (year/100 * 36524 + year/400) - 306; 886 outTime->sec = (((day - 1) * SEC_PER_DAY) - 62135596800) + hour*SEC_PER_HOUR + minute*SEC_PER_MINUTE + seconds; 887 888 // C's TM does not define a microsecond field. Microseconds must be manipulated by calling function. 494 889 outTime->usec = 0; 495 outTime->sec = localTMTime->tm_gmtoff + localTimeInSeconds;496 890 497 891 return outTime; 498 892 } 499 893 500 psTime* psT AIAdd(psTime *tai1, psTime *tai2)894 psTime* psTimeAdd(psTime *tai1, psTime *tai2) 501 895 { 502 896 psTime *outTime = NULL; 503 897 504 if(tai1->type != PS_TIME_TAI) { 505 psError(__func__, " : Line %d - First argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai1->type); 506 return NULL; 507 } 508 509 if(tai2->type != PS_TIME_TAI) { 510 psError(__func__, " : Line %d - Second argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai2->type); 898 899 if(tai1->type != tai2->type) { 900 psError(__func__, " : Line %d - Argument types inconsistent for time addition", __LINE__); 511 901 return NULL; 512 902 } … … 525 915 } 526 916 527 psTime* psT AISub(psTime *tai1, psTime *tai2)917 psTime* psTimeSub(psTime *tai1, psTime *tai2) 528 918 { 529 919 psTime *outTime = NULL; 530 920 531 if(tai1->type != PS_TIME_TAI) { 532 psError(__func__, " : Line %d - First argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai1->type); 533 return NULL; 534 } 535 536 if(tai2->type != PS_TIME_TAI) { 537 psError(__func__, " : Line %d - Second argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai2->type); 921 922 if(tai1->type != tai2->type) { 923 psError(__func__, " : Line %d - Argument types inconsistent for time subtraction", __LINE__); 538 924 return NULL; 539 925 } … … 552 938 } 553 939 554 psTime* psT AIDelta(psTime *tai1, psTime *tai2)940 psTime* psTimeDelta(psTime *tai1, psTime *tai2) 555 941 { 556 942 psTime *outTime = NULL; 557 943 558 if(tai1->type != PS_TIME_TAI) { 559 psError(__func__, " : Line %d - First argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai1->type); 560 return NULL; 561 } 562 563 if(tai2->type != PS_TIME_TAI) { 564 psError(__func__, " : Line %d - Second argument is not PS_TIME_TAI type. Type: %d", __LINE__, tai2->type); 944 945 if(tai1->type != tai2->type) { 946 psError(__func__, " : Line %d - Argument types inconsistent for time delta", __LINE__); 565 947 return NULL; 566 948 } … … 578 960 outTime->sec = abs(outTime->sec); 579 961 580 581 962 return outTime; 582 963 }
Note:
See TracChangeset
for help on using the changeset viewer.
