PLplot
5.10.0
|
00001 // 00002 // This software originally contributed under the LGPL in January 2009 to 00003 // PLplot by the 00004 // Cluster Science Centre 00005 // QSAS team, 00006 // Imperial College, London 00007 // Copyright (C) 2009 Imperial College, London 00008 // Copyright (C) 2009 Alan W. Irwin 00009 // 00010 // This file is part of PLplot. 00011 // 00012 // PLplot is free software; you can redistribute it and/or modify 00013 // it under the terms of the GNU Library General Public License as published 00014 // by the Free Software Foundation; either version 2 of the License, or 00015 // (at your option) any later version. 00016 // 00017 // PLplot is distributed in the hope that it will be useful, 00018 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00020 // GNU Library General Public License for more details. 00021 // 00022 // You should have received a copy of the GNU Library General Public License 00023 // along with PLplot; if not, write to the Free Software 00024 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00025 // 00026 // 00027 00028 // MJD measures from the start of 17 Nov 1858 00029 00030 // These utilities use the Gregorian calendar after 4 Oct 1582 (Julian) i.e. from 15 Oct 1582 Gregorian 00031 // Note C libraries use Gregorian only from 14 Sept 1752 00032 // More detailed discussion can be found at http://aa.usno.navy.mil/data/docs/JulianDate.php 00033 // These routines have been compared with the results of the US Naval Observatory online converter. 00034 // Modified Julian Date (MJD) = Julian Date (JD) - 2400000.5 00035 // 00036 // In all routines, specifying a day, hour, minute or second field greater than would be valid is 00037 // handled with modulo arithmetic and safe. 00038 // Thus 2006-12-32 00:62:00.0 will safely, and correctly, be treated as 2007-01-01 01:02:00.0 00039 // 00040 // 00041 #include <ctype.h> 00042 #include <math.h> 00043 #include "qsastimeP.h" 00044 #include "tai-utc.h" 00045 00046 // MJD for 0000-01-01 (correctly Jan 01, BCE 1) 00047 // Julian proleptic calendar value. 00048 #define MJD_0000J -678943 00049 // Gregorian proleptic calendar value. (At MJD_0000J the Gregorian proleptic 00050 // calendar reads two days behind the Julian proleptic calendar, i.e. - 2 days, 00051 // see http://en.wikipedia.org/wiki/Proleptic_Gregorian_calendar, 00052 // so MJD_0000G = MJD_0000J+2) 00053 #define MJD_0000G -678941 00054 // MJD for 0001-01-01 which is 366 days later than previous definitions because 00055 // the year 0 is a leap year in both calendars. 00056 #define MJD_0001J -678577 00057 #define MJD_0001G -678575 00058 // MJD for Jan 01, 1970 00:00:00 Gregorian, the Unix epoch. 00059 #define MJD_1970 40587 00060 00061 static const double SecInDay = 86400; // we ignore leap seconds 00062 static const int MonthStartDOY[] = { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 }; 00063 static const int MonthStartDOY_L[] = { 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 }; 00064 00065 // Static function declarations. 00066 static int geMJDtime_TAI( const MJDtime *number1, const TAI_UTC *number2 ); 00067 static int geMJDtime_UTC( const MJDtime *number1, const TAI_UTC *number2 ); 00068 static double leap_second_TAI( const MJDtime *MJD_TAI, int *inleap, int *index ); 00069 // End of static function declarations. 00070 00071 int setFromUT( int year, int month, int day, int hour, int min, double sec, MJDtime *MJD, int forceJulian ) 00072 { 00073 // convert broken-down time to MJD 00074 // MJD measures from the start of 17 Nov 1858 00075 00076 // If forceJulian is true (non-zero), the Julian proleptic calendar is 00077 // used whatever the year. Otherwise, the Gregorian proleptic calendar 00078 // is used whatever the year. 00079 // Note C libraries use Gregorian only (at least on Linux). In contrast, 00080 // the Linux (and Unix?) cal application uses Julian for earlier dates 00081 // and Gregorian from 14 Sept 1752 onwards. 00082 00083 int leaps, year4, year100, year400; 00084 double dbase_day, non_leaps = 365.; 00085 //int dbase_day, non_leaps = 365; 00086 double time_sec, dextraDays; 00087 int extraDays; 00088 00089 if ( month < 0 || month > 11 ) 00090 { 00091 fprintf( stderr, "setfromUT: invalid month value\n" ); 00092 exit( EXIT_FAILURE ); 00093 } 00094 // As year increases, year4/4 increments by 1 at 00095 // year = -7, -3, 1, 5, 9, etc. 00096 // As year increases, year100/100 increments by 1 at 00097 // year = -299, -199, -99, 1, 101, 201, 301, etc. 00098 // As year increases, year400/400 increments by 1 at 00099 // year = -1199, -799, -399, 1, 401, 801, 1201, etc. 00100 if ( year <= 0 ) 00101 { 00102 year4 = year - 4; 00103 year100 = year - 100; 00104 year400 = year - 400; 00105 } 00106 else 00107 { 00108 year4 = year - 1; 00109 year100 = year - 1; 00110 year400 = year - 1; 00111 } 00112 00113 if ( forceJulian ) 00114 { 00115 // count leap years on proleptic Julian Calendar starting from MJD_0000J 00116 leaps = year4 / 4; 00117 if ( year % 4 == 0 ) 00118 dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000J; 00119 else 00120 dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000J; 00121 } 00122 else 00123 { 00124 // count leap years for proleptic Gregorian Calendar. 00125 // Algorithm below for 1858-11-17 (0 MJD) gives 00126 // leaps = 450 and hence dbase_day of 678941, so subtract that value 00127 // or add MJD_0000G (which is two days different from MJD_0000J, see 00128 // above). 00129 leaps = year4 / 4 - year100 / 100 + year400 / 400; 00130 00131 // left to right associativity means the double value of 00132 // non_leaps propagate to make all calculations be 00133 // done in double precision without the potential of 00134 // integer overflow. The result should be a double which 00135 // stores the expected exact integer results of the 00136 // calculation with exact representation unless the 00137 // result is much larger than the integer overflow limit. 00138 if ( ( year % 4 == 0 && year % 100 != 0 ) || ( year % 4 == 0 && year % 400 == 0 ) ) 00139 dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000G; 00140 else 00141 dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000G; 00142 } 00143 00144 time_sec = sec + ( (double) min + (double) hour * 60. ) * 60.; 00145 00146 if ( time_sec >= SecInDay ) 00147 { 00148 dextraDays = ( time_sec / SecInDay ); 00149 // precaution against overflowing extraDays. 00150 if ( fabs( dextraDays ) > 2.e9 ) 00151 { 00152 return 3; 00153 } 00154 extraDays = (int) ( dextraDays ); 00155 dbase_day += extraDays; 00156 time_sec -= extraDays * SecInDay; 00157 } 00158 // precaution against overflowing MJD->base_day. 00159 if ( fabs( dbase_day ) > 2.e9 ) 00160 { 00161 return 4; 00162 } 00163 else 00164 { 00165 // The exact integer result should be represented exactly in the 00166 // double, dbase_day, and its absolute value should be less than 00167 // the integer overflow limit. So the cast to int should be 00168 // exact. 00169 MJD->base_day = (int) dbase_day; 00170 MJD->time_sec = time_sec; 00171 return 0; 00172 } 00173 } 00174 00175 void getYAD( int *year, int *ifleapyear, int *doy, const MJDtime *MJD, int forceJulian ) 00176 { 00177 // Get year and day of year from normalized MJD 00178 00179 int j, ifcorrect, year4, year100, year400; 00180 00181 j = MJD->base_day; 00182 00183 if ( forceJulian ) 00184 { 00185 // Shift j epoch to 0000-01-01 for the Julian proleptic calendar. 00186 j -= MJD_0000J; 00187 00188 // 365.25 is the exact period of the Julian year so year will be correct 00189 // if the day offset is set exactly right so that years -4, 0, 4 are 00190 // leap years, i.e. years -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 start with 00191 // j = -1826 -1461, -1095, -730, -365, 0, 366, 731, 1096, 1461, 1827 00192 if ( j >= 366 ) 00193 { 00194 *year = (int) ( (double) ( j ) / 365.25 ); 00195 year4 = *year - 1; 00196 } 00197 else 00198 { 00199 *year = (int) ( (double) ( j - 365 ) / 365.25 ); 00200 year4 = *year - 4; 00201 } 00202 00203 *doy = j - *year * 365 - year4 / 4; 00204 00205 *ifleapyear = *year % 4 == 0; 00206 } 00207 else 00208 { 00209 // Shift j epoch to 0000-01-01 for the Gregorian proleptic calendar. 00210 j -= MJD_0000G; 00211 // 365.245 is the exact period of the Gregorian year so year will be correct 00212 // on average, but because the leap year rule is irregular within 00213 // the 400-year Gregorian cycle, the first and last days of the 00214 // year may need further adjustment. 00215 00216 if ( j >= 366 ) 00217 { 00218 *year = (int) ( (double) ( j ) / 365.2425 ); 00219 year4 = *year - 1; 00220 year100 = *year - 1; 00221 year400 = *year - 1; 00222 } 00223 else 00224 { 00225 *year = (int) ( (double) ( j - 365 ) / 365.2425 ); 00226 year4 = *year - 4; 00227 year100 = *year - 100; 00228 year400 = *year - 400; 00229 } 00230 00231 *doy = j - *year * 365 - year4 / 4 + year100 / 100 - year400 / 400; 00232 *ifleapyear = ( *year % 4 == 0 && *year % 100 != 0 ) || ( *year % 4 == 0 && *year % 400 == 0 ); 00233 00234 // Rare corrections to above average Gregorian relations. 00235 if ( *doy < 1 ) 00236 { 00237 ( *year )--; 00238 ifcorrect = 1; 00239 } 00240 else if ( *doy > 365 && ( !*ifleapyear || *doy > 366 ) ) 00241 { 00242 ( *year )++; 00243 ifcorrect = 1; 00244 } 00245 else 00246 { 00247 ifcorrect = 0; 00248 } 00249 if ( ifcorrect ) 00250 { 00251 if ( j >= 366 ) 00252 { 00253 year4 = *year - 1; 00254 year100 = *year - 1; 00255 year400 = *year - 1; 00256 } 00257 else 00258 { 00259 year4 = *year - 4; 00260 year100 = *year - 100; 00261 year400 = *year - 400; 00262 } 00263 00264 *doy = j - *year * 365 - year4 / 4 + year100 / 100 - year400 / 400; 00265 *ifleapyear = ( *year % 4 == 0 && *year % 100 != 0 ) || ( *year % 4 == 0 && *year % 400 == 0 ); 00266 } 00267 } 00268 } 00269 00270 void normalize_MJD( MJDtime *MJD ) 00271 { 00272 int extra_days; 00273 // Calculate MJDout as normalized version 00274 // (i.e., 0. <= MJDout->time_sec < 86400.) of MJDin. 00275 if ( MJD->time_sec >= 0 ) 00276 { 00277 extra_days = (int) ( MJD->time_sec / SecInDay ); 00278 } 00279 else 00280 { 00281 // allow for negative seconds push into previous day even if less than 1 day 00282 extra_days = (int) ( MJD->time_sec / SecInDay ) - 1; 00283 } 00284 00285 MJD->base_day += extra_days; 00286 MJD->time_sec -= extra_days * SecInDay; 00287 } 00288 00289 void breakDownMJD( int *year, int *month, int *day, int *hour, int *min, double *sec, const MJDtime *MJD, int forceJulian ) 00290 { 00291 // Convert MJD struct into date/time elements 00292 // Note year 0 CE (AD) [1 BCE (BC)] is a leap year 00293 00294 int doy, ifleapyear; 00295 MJDtime nMJD_value, *nMJD = &nMJD_value; 00296 00297 *nMJD = *MJD; 00298 normalize_MJD( nMJD ); 00299 // Time part 00300 00301 *sec = nMJD->time_sec; 00302 *hour = (int) ( *sec / 3600. ); 00303 *sec -= (double) *hour * 3600.; 00304 *min = (int) ( *sec / 60. ); 00305 *sec -= (double) *min * 60.; 00306 00307 getYAD( year, &ifleapyear, &doy, nMJD, forceJulian ); 00308 00309 // calculate month part with doy set to be the day within 00310 // the year in the range from 1 to 366 00311 *month = -1; 00312 if ( ifleapyear ) 00313 { 00314 while ( doy > MonthStartDOY_L[*month + 1] ) 00315 { 00316 ( *month )++; 00317 if ( *month == 11 ) 00318 break; 00319 } 00320 *day = doy - MonthStartDOY_L[*month]; 00321 } 00322 else 00323 { 00324 while ( doy > MonthStartDOY[*month + 1] ) 00325 { 00326 ( *month )++; 00327 if ( *month == 11 ) 00328 break; 00329 } 00330 *day = doy - MonthStartDOY[*month]; 00331 } 00332 } 00333 00334 const char * getDayOfWeek( const MJDtime *MJD ) 00335 { 00336 static const char *dow = { "Wed\0Thu\0Fri\0Sat\0Sun\0Mon\0Tue" }; 00337 int d = MJD->base_day % 7; 00338 if ( d < 0 ) 00339 d += 7; 00340 return &( dow[d * 4] ); 00341 } 00342 00343 const char * getLongDayOfWeek( const MJDtime *MJD ) 00344 { 00345 static const char *dow = { "Wednesday\0Thursday\0\0Friday\0\0\0\0Saturday\0\0Sunday\0\0\0\0Monday\0\0\0\0Tuesday" }; 00346 int d = MJD->base_day % 7; 00347 if ( d < 0 ) 00348 d += 7; 00349 return &( dow[d * 10] ); 00350 } 00351 00352 const char * getMonth( int m ) 00353 { 00354 static const char *months = { "Jan\0Feb\0Mar\0Apr\0May\0Jun\0Jul\0Aug\0Sep\0Oct\0Nov\0Dec" }; 00355 return &( months[( m ) * 4] ); 00356 } 00357 00358 const char * getLongMonth( int m ) 00359 { 00360 static const char *months = { "January\0\0\0February\0\0March\0\0\0\0\0April\0\0\0\0\0May\0\0\0\0\0\0\0June\0\0\0\0\0\0July\0\0\0\0\0\0August\0\0\0\0September\0October\0\0\0November\0\0December" }; 00361 return &( months[( m ) * 10] ); 00362 } 00363 00364 00365 size_t strfMJD( char * buf, size_t len, const char *format, const MJDtime *MJD, int forceJulian, int inleap ) 00366 { 00367 // Format a text string according to the format string. 00368 // Uses the same syntax as strftime() but does not use current locale. 00369 // The null terminator is included in len for safety. 00370 00371 // if inleap is true (non-zero) then renormalize so that (int) sec 00372 // is 60 to mark results as a flag that a leap increment (recently 00373 // always a second, but historically it was sometimes smaller than 00374 // that) was in the process of being inserted for this particular 00375 // epoch just prior to a positive discontinuity in TAI-UTC. 00376 00377 int year, month, day, hour, min, ysign, second, d, y; 00378 int y1, ifleapyear; 00379 int i, secsSince1970; 00380 int nplaces, fmtlen, slen; 00381 int resolution; 00382 double shiftPlaces; 00383 char * ptr; 00384 double sec, sec_fraction; 00385 int w, doy, days_in_wk1; 00386 const char *dayText; 00387 const char *monthText; 00388 char DateTime[80]; 00389 size_t posn = 0; 00390 size_t last = len - 1; 00391 MJDtime nMJD_value, *nMJD = &nMJD_value; 00392 char dynamic_format[10]; 00393 00394 // Find required resolution 00395 resolution = 0; 00396 fmtlen = (int) strlen( format ); 00397 i = 0; 00398 while ( i < fmtlen ) 00399 { 00400 char next = format[i]; 00401 if ( next == '%' ) 00402 { 00403 // find seconds format if used 00404 i++; 00405 next = format[i]; 00406 if ( isdigit( next ) != 0 ) 00407 { 00408 nplaces = (int) strtol( &( format[i] ), NULL, 10 ); 00409 if ( nplaces > resolution ) 00410 resolution = nplaces; 00411 } 00412 else if ( next == '.' ) 00413 { 00414 resolution = 9; // maximum resolution allowed 00415 } 00416 } 00417 i++; 00418 } 00419 00420 // ensure rounding is done before breakdown 00421 shiftPlaces = pow( 10, (double) resolution ); 00422 *nMJD = *MJD; 00423 nMJD->time_sec += 0.5 / shiftPlaces; 00424 00425 buf[last] = '\0'; 00426 buf[0] = '\0'; // force overwrite of old buffer since strnctat() used hereafter 00427 00428 if ( inleap ) 00429 nMJD->time_sec -= 1.; 00430 00431 breakDownMJD( &year, &month, &day, &hour, &min, &sec, nMJD, forceJulian ); 00432 if ( inleap ) 00433 sec += 1.; 00434 00435 if ( year < 0 ) 00436 { 00437 ysign = 1; 00438 year = -year; 00439 } 00440 else 00441 ysign = 0; 00442 00443 //truncate seconds to resolution to stop formatting rounding up 00444 sec = floor( sec * shiftPlaces ) / shiftPlaces; 00445 second = (int) sec; 00446 00447 // Read format string, character at a time 00448 i = 0; 00449 while ( i < fmtlen ) 00450 { 00451 char next = format[i]; 00452 if ( next == '%' ) 00453 { 00454 // format character or escape 00455 i++; 00456 next = format[i]; 00457 if ( next == '%' ) 00458 { 00459 // escaped %, pass it on 00460 buf[posn] = next; 00461 posn++; 00462 if ( posn >= last ) 00463 return posn; 00464 } 00465 else if ( next == 'a' ) 00466 { 00467 // short day name 00468 dayText = getDayOfWeek( nMJD ); 00469 strncat( &( buf[posn] ), dayText, last - posn ); 00470 posn = strlen( buf ); 00471 if ( posn >= last ) 00472 return posn; 00473 } 00474 else if ( next == 'A' ) 00475 { 00476 // long day name 00477 dayText = getLongDayOfWeek( nMJD ); 00478 strncat( &( buf[posn] ), dayText, last - posn ); 00479 posn = strlen( buf ); 00480 if ( posn >= last ) 00481 return posn; 00482 } 00483 else if ( next == 'b' || next == 'h' ) 00484 { 00485 // short month name 00486 monthText = getMonth( month ); 00487 strncat( &( buf[posn] ), monthText, last - posn ); 00488 posn = strlen( buf ); 00489 if ( posn >= last ) 00490 return posn; 00491 } 00492 else if ( next == 'B' ) 00493 { 00494 // long month name 00495 monthText = getLongMonth( month ); 00496 strncat( &( buf[posn] ), monthText, last - posn ); 00497 posn = strlen( buf ); 00498 if ( posn >= last ) 00499 return posn; 00500 } 00501 else if ( next == 'c' ) 00502 { 00503 // Date and Time with day of week 00504 dayText = getDayOfWeek( nMJD ); 00505 monthText = getMonth( month ); 00506 if ( ysign == 0 ) 00507 sprintf( DateTime, "%s %s %02d %02d:%02d:%02d %04d", dayText, monthText, day, hour, min, second, year ); 00508 else 00509 sprintf( DateTime, "%s %s %02d %02d:%02d:%02d -%04d", dayText, monthText, day, hour, min, second, year ); 00510 00511 strncat( &( buf[posn] ), DateTime, last - posn ); 00512 posn = strlen( buf ); 00513 if ( posn >= last ) 00514 return posn; 00515 } 00516 else if ( next == 'C' ) 00517 { 00518 // year / 100 so, e.g. 1989 is 20th century but comes out as 19 00519 int century = year / 100; 00520 if ( ysign == 0 ) 00521 sprintf( DateTime, "%02d", century ); 00522 else 00523 sprintf( DateTime, "-%02d", century + 1 ); 00524 00525 strncat( &( buf[posn] ), DateTime, last - posn ); 00526 posn = strlen( buf ); 00527 if ( posn >= last ) 00528 return posn; 00529 } 00530 else if ( next == 'd' ) 00531 { 00532 // day of month (01 - 31) 00533 sprintf( DateTime, "%02d", day ); 00534 00535 strncat( &( buf[posn] ), DateTime, last - posn ); 00536 posn = strlen( buf ); 00537 if ( posn >= last ) 00538 return posn; 00539 } 00540 else if ( next == 'D' ) 00541 { 00542 // month/day/year 00543 y = year % 100; 00544 if ( ysign == 0 ) 00545 sprintf( DateTime, "%02d/%02d/%02d", month + 1, day, y ); 00546 else 00547 sprintf( DateTime, "%02d/%02d/-%02d", month + 1, day, y ); 00548 00549 strncat( &( buf[posn] ), DateTime, last - posn ); 00550 posn = strlen( buf ); 00551 if ( posn >= last ) 00552 return posn; 00553 } 00554 else if ( next == 'e' ) 00555 { 00556 // day of month ( 1 - 31) 00557 if ( day < 10 ) 00558 sprintf( DateTime, " %01d", day ); 00559 else 00560 sprintf( DateTime, "%02d", day ); 00561 00562 strncat( &( buf[posn] ), DateTime, last - posn ); 00563 posn = strlen( buf ); 00564 if ( posn >= last ) 00565 return posn; 00566 } 00567 else if ( next == 'F' ) 00568 { 00569 // year-month-day 00570 if ( ysign == 0 ) 00571 sprintf( DateTime, "%04d-%02d-%02d", year, month + 1, day ); 00572 else 00573 sprintf( DateTime, "-%04d-%02d-%02d", year, month + 1, day ); 00574 00575 strncat( &( buf[posn] ), DateTime, last - posn ); 00576 posn = strlen( buf ); 00577 if ( posn >= last ) 00578 return posn; 00579 } 00580 else if ( next == 'H' ) 00581 { 00582 // hour, 24 hour clock (00 - 23) 00583 sprintf( DateTime, "%02d", hour ); 00584 00585 strncat( &( buf[posn] ), DateTime, last - posn ); 00586 posn = strlen( buf ); 00587 if ( posn >= last ) 00588 return posn; 00589 } 00590 else if ( next == 'I' ) 00591 { 00592 // hour, 12 hour clock (01 - 12) 00593 if ( hour == 0 ) 00594 sprintf( DateTime, "%02d", hour + 12 ); 00595 else if ( hour > 12 ) 00596 sprintf( DateTime, "%02d", hour - 12 ); 00597 else 00598 sprintf( DateTime, "%02d", hour ); 00599 00600 strncat( &( buf[posn] ), DateTime, last - posn ); 00601 posn = strlen( buf ); 00602 if ( posn >= last ) 00603 return posn; 00604 } 00605 else if ( next == 'j' ) 00606 { 00607 // day of year 00608 getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian ); 00609 sprintf( DateTime, "%03d", doy ); 00610 00611 strncat( &( buf[posn] ), DateTime, last - posn ); 00612 posn = strlen( buf ); 00613 if ( posn >= last ) 00614 return posn; 00615 } 00616 else if ( next == 'k' ) 00617 { 00618 // hour, 24 hour clock ( 0 - 23) 00619 if ( hour < 10 ) 00620 sprintf( DateTime, " %01d", hour ); 00621 else 00622 sprintf( DateTime, "%02d", hour ); 00623 00624 strncat( &( buf[posn] ), DateTime, last - posn ); 00625 posn = strlen( buf ); 00626 if ( posn >= last ) 00627 return posn; 00628 } 00629 else if ( next == 'l' ) 00630 { 00631 // hour, 12 hour clock ( 1 - 12) 00632 if ( hour == 0 ) 00633 sprintf( DateTime, "%02d", hour + 12 ); 00634 else if ( hour < 10 ) 00635 sprintf( DateTime, " %01d", hour ); 00636 else if ( hour <= 12 ) 00637 sprintf( DateTime, "%02d", hour ); 00638 else if ( hour < 22 ) 00639 sprintf( DateTime, " %01d", hour - 12 ); 00640 else 00641 sprintf( DateTime, "%02d", hour - 12 ); 00642 00643 strncat( &( buf[posn] ), DateTime, last - posn ); 00644 posn = strlen( buf ); 00645 if ( posn >= last ) 00646 return posn; 00647 } 00648 else if ( next == 'm' ) 00649 { 00650 // month (01 - 12) 00651 sprintf( DateTime, "%02d", month + 1 ); 00652 00653 strncat( &( buf[posn] ), DateTime, last - posn ); 00654 posn = strlen( buf ); 00655 if ( posn >= last ) 00656 return posn; 00657 } 00658 else if ( next == 'M' ) 00659 { 00660 // minute (00 - 59) 00661 sprintf( DateTime, "%02d", min ); 00662 00663 strncat( &( buf[posn] ), DateTime, last - posn ); 00664 posn = strlen( buf ); 00665 if ( posn >= last ) 00666 return posn; 00667 } 00668 else if ( next == 'n' ) 00669 { 00670 // newline 00671 buf[posn] = '\n'; 00672 posn++; 00673 if ( posn >= last ) 00674 return posn; 00675 } 00676 else if ( next == 'p' ) 00677 { 00678 // am/pm on12 hour clock 00679 if ( hour < 0 ) 00680 sprintf( DateTime, "AM" ); 00681 else 00682 sprintf( DateTime, "PM" ); 00683 00684 strncat( &( buf[posn] ), DateTime, last - posn ); 00685 posn = strlen( buf ); 00686 if ( posn >= last ) 00687 return posn; 00688 } 00689 else if ( next == 'r' ) 00690 { 00691 // hour:min:sec AM, 12 hour clock (01 - 12):(00 - 59):(00 - 59) (AM - PM) 00692 if ( hour == 0 ) 00693 sprintf( DateTime, "%02d:%02d:%02d AM", hour + 12, min, second ); 00694 else if ( hour > 12 ) 00695 sprintf( DateTime, "%02d:%02d:%02d PM", hour - 12, min, second ); 00696 else if ( hour == 12 ) 00697 sprintf( DateTime, "%02d:%02d:%02d PM", hour, min, second ); 00698 else 00699 sprintf( DateTime, "%02d:%02d:%02d AM", hour, min, second ); 00700 00701 strncat( &( buf[posn] ), DateTime, last - posn ); 00702 posn = strlen( buf ); 00703 if ( posn >= last ) 00704 return posn; 00705 } 00706 else if ( next == 'R' ) 00707 { 00708 // hour:min, 24 hour clock (00 - 23):(00 - 59) 00709 sprintf( DateTime, "%02d:%02d", hour, min ); 00710 00711 strncat( &( buf[posn] ), DateTime, last - posn ); 00712 posn = strlen( buf ); 00713 if ( posn >= last ) 00714 return posn; 00715 } 00716 else if ( next == 'S' ) 00717 { 00718 // second (00 - 59 with optional decimal point and numbers after the decimal point.) 00719 if ( i + 2 < fmtlen && format[i + 1] == '%' && ( format[i + 2] == '.' || isdigit( format[i + 2] ) != 0 ) ) 00720 { 00721 // nplaces is number of decimal places ( 0 < nplaces <= 9 ) 00722 if ( format[i + 2] == '.' ) 00723 // maximum number of places 00724 nplaces = 9; 00725 else 00726 nplaces = (int) strtol( &( format[i + 2] ), NULL, 10 ); 00727 i += 2; 00728 } 00729 else 00730 { 00731 nplaces = 0; 00732 } 00733 00734 if ( nplaces == 0 ) 00735 { 00736 sprintf( DateTime, "%02d", (int) ( sec + 0.5 ) ); 00737 } 00738 else 00739 { 00740 sprintf( dynamic_format, "%%0%d.%df", nplaces + 3, nplaces ); 00741 sprintf( DateTime, dynamic_format, sec ); 00742 if ( format[i] == '.' ) 00743 { 00744 slen = (int) strlen( DateTime ) - 1; 00745 while ( DateTime[slen] == '0' && DateTime[slen - 1] != '.' ) 00746 { 00747 DateTime[slen] = '\0'; // remove trailing zeros 00748 slen--; 00749 } 00750 } 00751 } 00752 00753 strncat( &( buf[posn] ), DateTime, last - posn ); 00754 posn = strlen( buf ); 00755 if ( posn >= last ) 00756 return posn; 00757 } 00758 else if ( next == 's' ) 00759 { 00760 // seconds since 01 Jan 1970 Gregorian 00761 secsSince1970 = (int) ( nMJD->time_sec + ( nMJD->base_day - MJD_1970 ) * SecInDay ); 00762 sprintf( DateTime, "%d", secsSince1970 ); 00763 00764 strncat( &( buf[posn] ), DateTime, last - posn ); 00765 posn = strlen( buf ); 00766 if ( posn >= last ) 00767 return posn; 00768 } 00769 else if ( next == 't' ) 00770 { 00771 // tab 00772 buf[posn] = '\t'; 00773 posn++; 00774 if ( posn >= last ) 00775 return posn; 00776 } 00777 else if ( next == 'T' ) 00778 { 00779 // hour:min:sec, 24 hour clock (00 - 23):(00 - 59):(00 - 59) 00780 sprintf( DateTime, "%02d:%02d:%02d", hour, min, second ); 00781 00782 strncat( &( buf[posn] ), DateTime, last - posn ); 00783 posn = strlen( buf ); 00784 if ( posn >= last ) 00785 return posn; 00786 } 00787 else if ( next == 'U' ) 00788 { 00789 // week of year as a number, (00 - 53) start of week is Sunday 00790 getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian ); 00791 days_in_wk1 = ( nMJD->base_day - doy - 4 ) % 7; 00792 00793 w = ( doy + 6 - days_in_wk1 ) / 7; 00794 00795 sprintf( DateTime, "%02d", w ); 00796 00797 strncat( &( buf[posn] ), DateTime, last - posn ); 00798 posn = strlen( buf ); 00799 if ( posn >= last ) 00800 return posn; 00801 } 00802 else if ( next == 'u' ) 00803 { 00804 // weekday as a number, 0 = Monday 00805 d = 1 + ( nMJD->base_day - 5 ) % 7; 00806 00807 sprintf( DateTime, "%01d", d ); 00808 00809 strncat( &( buf[posn] ), DateTime, last - posn ); 00810 posn = strlen( buf ); 00811 if ( posn >= last ) 00812 return posn; 00813 } 00814 else if ( next == 'v' ) 00815 { 00816 // day-MonthName-year day of month ( 1 - 31) - month (Jan ... Dec) - year (yyyy) 00817 00818 monthText = getMonth( month ); 00819 00820 if ( ysign == 0 ) 00821 { 00822 if ( day < 10 ) 00823 sprintf( DateTime, " %01d-%s-%04d", day, monthText, year ); 00824 else 00825 sprintf( DateTime, "%02d-%s-%04d", day, monthText, year ); 00826 } 00827 else 00828 { 00829 if ( day < 10 ) 00830 sprintf( DateTime, " %01d-%s-(-)%04d", day, monthText, year ); 00831 else 00832 sprintf( DateTime, "%02d-%s-(-)%04d", day, monthText, year ); 00833 } 00834 00835 strncat( &( buf[posn] ), DateTime, last - posn ); 00836 posn = strlen( buf ); 00837 if ( posn >= last ) 00838 return posn; 00839 } 00840 else if ( next == 'V' ) 00841 { 00842 // week of year as a number, (01 - 53) start of week is Monday and first week has at least 3 days in year 00843 getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian ); 00844 days_in_wk1 = ( nMJD->base_day - doy - 3 ) % 7; 00845 00846 if ( days_in_wk1 <= 3 ) 00847 w = ( doy + 6 - days_in_wk1 ) / 7; // ensure first week has at least 3 days in this year 00848 else 00849 w = 1 + ( doy + 6 - days_in_wk1 ) / 7; 00850 00851 if ( w == 0 ) 00852 w = 53; 00853 sprintf( DateTime, "%02d", w ); 00854 00855 strncat( &( buf[posn] ), DateTime, last - posn ); 00856 posn = strlen( buf ); 00857 if ( posn >= last ) 00858 return posn; 00859 } 00860 else if ( next == 'w' ) 00861 { 00862 // weekday as a number, 0 = Sunday 00863 d = ( nMJD->base_day - 4 ) % 7; 00864 00865 sprintf( DateTime, "%01d", d ); 00866 00867 strncat( &( buf[posn] ), DateTime, last - posn ); 00868 posn = strlen( buf ); 00869 if ( posn >= last ) 00870 return posn; 00871 } 00872 else if ( next == 'W' ) 00873 { 00874 // week of year as a number, (00 - 53) start of week is Monday 00875 getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian ); 00876 days_in_wk1 = ( nMJD->base_day - doy - 3 ) % 7; 00877 00878 w = ( doy + 6 - days_in_wk1 ) / 7; 00879 00880 sprintf( DateTime, "%02d", w ); 00881 00882 strncat( &( buf[posn] ), DateTime, last - posn ); 00883 posn = strlen( buf ); 00884 if ( posn >= last ) 00885 return posn; 00886 } 00887 else if ( next == 'x' ) 00888 { 00889 // date string 00890 dayText = getDayOfWeek( nMJD ); 00891 monthText = getMonth( month ); 00892 if ( ysign == 0 ) 00893 sprintf( DateTime, "%s %s %02d, %04d", dayText, monthText, day, year ); 00894 else 00895 sprintf( DateTime, "%s %s %02d, -%04d", dayText, monthText, day, year ); 00896 00897 strncat( &( buf[posn] ), DateTime, last - posn ); 00898 posn = strlen( buf ); 00899 if ( posn >= last ) 00900 return posn; 00901 } 00902 else if ( next == 'X' ) 00903 { 00904 // time string 00905 sprintf( DateTime, "%02d:%02d:%02d", hour, min, second ); 00906 00907 strncat( &( buf[posn] ), DateTime, last - posn ); 00908 posn = strlen( buf ); 00909 if ( posn >= last ) 00910 return posn; 00911 } 00912 else if ( next == 'y' ) 00913 { 00914 // 2 digit year 00915 y = year % 100; 00916 00917 if ( ysign == 0 ) 00918 sprintf( DateTime, "%02d", y ); 00919 else 00920 sprintf( DateTime, "-%02d", y ); 00921 00922 strncat( &( buf[posn] ), DateTime, last - posn ); 00923 posn = strlen( buf ); 00924 if ( posn >= last ) 00925 return posn; 00926 } 00927 else if ( next == 'Y' ) 00928 { 00929 // 4 digit year 00930 if ( ysign == 0 ) 00931 sprintf( DateTime, "%04d", year ); 00932 else 00933 sprintf( DateTime, "-%04d", year ); 00934 00935 strncat( &( buf[posn] ), DateTime, last - posn ); 00936 posn = strlen( buf ); 00937 if ( posn >= last ) 00938 return posn; 00939 } 00940 else if ( next == 'Z' ) 00941 { 00942 // time zone and calendar, always UTC 00943 if ( forceJulian ) 00944 strncat( &( buf[posn] ), "UTC Julian", last - posn ); 00945 else 00946 strncat( &( buf[posn] ), "UTC Gregorian", last - posn ); 00947 00948 posn = strlen( buf ); 00949 if ( posn >= last ) 00950 return posn; 00951 } 00952 else if ( next == 'z' ) 00953 { 00954 // time zone, always UTC 00955 strncat( &( buf[posn] ), "+0000", last - posn ); 00956 posn = strlen( buf ); 00957 if ( posn >= last ) 00958 return posn; 00959 } 00960 else if ( next == '+' ) 00961 { 00962 // date and time 00963 dayText = getDayOfWeek( nMJD ); 00964 monthText = getMonth( month ); 00965 if ( ysign == 0 ) 00966 sprintf( DateTime, "%s %s %02d %02d:%02d:%02d UTC %04d", dayText, monthText, day, hour, min, second, year ); 00967 else 00968 sprintf( DateTime, "%s %s %02d %02d:%02d:%02d UTC -%04d", dayText, monthText, day, hour, min, second, year ); 00969 00970 strncat( &( buf[posn] ), DateTime, last - posn ); 00971 posn = strlen( buf ); 00972 if ( posn >= last ) 00973 return posn; 00974 } 00975 else if ( next == '.' || isdigit( next ) != 0 ) 00976 { 00977 // nplaces is number of decimal places ( 0 < nplaces <= 9 ) 00978 if ( next == '.' ) 00979 // maximum number of places 00980 nplaces = 9; 00981 else 00982 nplaces = (int) strtol( &( format[i] ), NULL, 10 ); 00983 00984 // fractional part of seconds to maximum available accuracy 00985 sec_fraction = sec - (int) sec; 00986 sprintf( dynamic_format, "%%-%d.%df", nplaces + 2, nplaces ); 00987 //sprintf(DateTime, "%-11.9f", sec_fraction); 00988 sprintf( DateTime, dynamic_format, sec_fraction ); 00989 while ( ( ptr = strrchr( &( DateTime[0] ), ' ' ) ) != NULL ) 00990 ptr[0] = '\0'; // remove trailing white space 00991 00992 if ( next == '.' ) 00993 { 00994 slen = (int) strlen( DateTime ) - 1; 00995 while ( DateTime[slen] == '0' && DateTime[slen - 1] != '.' ) 00996 { 00997 DateTime[slen] = '\0'; // remove trailing zeros 00998 slen--; 00999 } 01000 } 01001 01002 ptr = strchr( DateTime, '.' ); 01003 // remove everything in front of the decimal point, and 01004 // ignore case (%0) where no decimal point exists at all. 01005 if ( ptr != NULL ) 01006 strncat( &( buf[posn] ), ptr, last - posn ); 01007 posn = strlen( buf ); 01008 if ( posn >= last ) 01009 return posn; 01010 } 01011 } 01012 else 01013 { 01014 // regular multi-byte character, pass it on 01015 buf[posn] = next; 01016 posn++; 01017 if ( posn >= last ) 01018 return posn; 01019 } 01020 buf[posn] = '\0'; 01021 i++; 01022 } 01023 return posn; 01024 } 01025 01026 int geMJDtime_TAI( const MJDtime *number1, const TAI_UTC *number2 ) 01027 { 01028 // Returns true if number1 >= number2. 01029 // N.B. both number1 and number2 must be normalized. 01030 if ( number1->base_day > number2->base_day ) 01031 { 01032 return 1; 01033 } 01034 else if ( number1->base_day < number2->base_day ) 01035 { 01036 return 0; 01037 } 01038 else 01039 { 01040 return ( number1->time_sec >= number2->time_sec_tai ); 01041 } 01042 } 01043 01044 int geMJDtime_UTC( const MJDtime *number1, const TAI_UTC *number2 ) 01045 { 01046 // Returns true if number1 >= number2. 01047 // N.B. both number1 and number2 must be normalized. 01048 if ( number1->base_day > number2->base_day ) 01049 { 01050 return 1; 01051 } 01052 else if ( number1->base_day < number2->base_day ) 01053 { 01054 return 0; 01055 } 01056 else 01057 { 01058 return ( number1->time_sec >= number2->time_sec_utc ); 01059 } 01060 } 01061 01062 double leap_second_TAI( const MJDtime *MJD_TAI, int *inleap, int *index ) 01063 { 01064 // Logic assumes input MJD_TAI is in TAI 01065 // *inleap lets the calling routine know whether MJD_TAI corresponds 01066 // to an epoch when a positive leap increment is being inserted. 01067 01068 MJDtime MJD_value, *MJD = &MJD_value; 01069 double leap; 01070 int debug = 0; 01071 // N.B. geMJDtime_TAI only works for normalized values. 01072 *MJD = *MJD_TAI; 01073 normalize_MJD( MJD ); 01074 // Search for index such that TAI_UTC_lookup_table[*index] <= MJD(TAI) < TAI_UTC_lookup_table[*index+1] 01075 bhunt_search( MJD, TAI_UTC_lookup_table, number_of_entries_in_tai_utc_table, sizeof ( TAI_UTC ), index, ( int ( * )( const void *, const void * ) )geMJDtime_TAI ); 01076 if ( debug == 2 ) 01077 fprintf( stderr, "*index = %d\n", *index ); 01078 if ( *index == -1 ) 01079 { 01080 // MJD is less than first table entry. 01081 // Debug: check that condition is met 01082 if ( debug && geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index + 1] ) ) 01083 { 01084 fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index ); 01085 exit( EXIT_FAILURE ); 01086 } 01087 // There is (by assertion) no discontinuity at the start of the table. 01088 // Therefore, *inleap cannot be true. 01089 *inleap = 0; 01090 // Use initial offset for MJD values before first table entry. 01091 // Calculate this offset strictly from offset1. The slope term 01092 // doesn't enter because offset2 is the same as the UTC of the 01093 // first epoch of the table. 01094 return -TAI_UTC_lookup_table[*index + 1].offset1; 01095 } 01096 else if ( *index == number_of_entries_in_tai_utc_table - 1 ) 01097 { 01098 // MJD is greater than or equal to last table entry. 01099 // Debug: check that condition is met 01100 if ( debug && !geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index] ) ) 01101 { 01102 fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index ); 01103 exit( EXIT_FAILURE ); 01104 } 01105 // If beyond end of table, cannot be in middle of leap second insertion. 01106 *inleap = 0; 01107 // Use final offset for MJD values after last table entry. 01108 // The slope term doesn't enter because modern values of the slope 01109 // are zero. 01110 return -TAI_UTC_lookup_table[*index].offset1; 01111 } 01112 else if ( *index >= 0 && *index < number_of_entries_in_tai_utc_table ) 01113 { 01114 // table[*index] <= MJD < table[*index+1]. 01115 // Debug: check that condition is met 01116 if ( debug && !( geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index] ) && !geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index + 1] ) ) ) 01117 { 01118 fprintf( stderr, "MJD = {%d, %f}\n", MJD->base_day, MJD->time_sec ); 01119 fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index ); 01120 exit( EXIT_FAILURE ); 01121 } 01122 leap = -( TAI_UTC_lookup_table[*index].offset1 + ( ( MJD->base_day - TAI_UTC_lookup_table[*index].offset2 ) + MJD->time_sec / SecInDay ) * TAI_UTC_lookup_table[*index].slope ) / ( 1. + TAI_UTC_lookup_table[*index].slope / SecInDay ); 01123 // Convert MJD(TAI) to normalized MJD(UTC). 01124 MJD->time_sec += leap; 01125 normalize_MJD( MJD ); 01126 // If MJD(UT) is in the next interval of the corresponding 01127 // TAI_UTC_lookup_table, then we are right in the middle of a 01128 // leap interval (recently a second but for earlier epochs it could be 01129 // less) insertion. Note this logic even works when leap intervals 01130 // are taken away from UTC (i.e., leap is positive) since in that 01131 // case the UTC *index always corresponds to the TAI *index. 01132 *inleap = geMJDtime_UTC( MJD, &TAI_UTC_lookup_table[*index + 1] ); 01133 return leap; 01134 } 01135 else 01136 { 01137 fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad *index = %d\n", *index ); 01138 exit( EXIT_FAILURE ); 01139 } 01140 } 01141 01142 void configqsas( double scale, double offset1, double offset2, int ccontrol, int ifbtime_offset, int year, int month, int day, int hour, int min, double sec, QSASConfig **qsasconfig ) 01143 { 01144 // Configure the transformation between continuous time and broken-down time 01145 // that is used for ctimeqsas, btimeqsas, and strfqsas. 01146 int forceJulian, ret; 01147 MJDtime MJD_value, *MJD = &MJD_value; 01148 01149 // Allocate memory for *qsasconfig if that hasn't been done by a 01150 // previous call. 01151 if ( *qsasconfig == NULL ) 01152 { 01153 *qsasconfig = (QSASConfig *) malloc( (size_t) sizeof ( QSASConfig ) ); 01154 if ( *qsasconfig == NULL ) 01155 { 01156 fprintf( stderr, "configqsas: out of memory\n" ); 01157 exit( EXIT_FAILURE ); 01158 } 01159 } 01160 01161 // Set bhunt_search index to a definite value less than -1 to insure no 01162 // initial hunt phase from some random index value is done. 01163 ( *qsasconfig )->index = -40; 01164 01165 if ( scale != 0. ) 01166 { 01167 if ( ifbtime_offset ) 01168 { 01169 if ( ccontrol & 0x1 ) 01170 forceJulian = 1; 01171 else 01172 forceJulian = 0; 01173 ret = setFromUT( year, month, day, hour, min, sec, MJD, forceJulian ); 01174 if ( ret ) 01175 { 01176 fprintf( stderr, "configqsas: some problem with broken-down arguments\n" ); 01177 exit( EXIT_FAILURE ); 01178 } 01179 offset1 = (double) MJD->base_day; 01180 offset2 = MJD->time_sec / (double) SecInDay; 01181 } 01182 ( *qsasconfig )->scale = scale; 01183 ( *qsasconfig )->offset1 = offset1; 01184 ( *qsasconfig )->offset2 = offset2; 01185 ( *qsasconfig )->ccontrol = ccontrol; 01186 } 01187 else 01188 { 01189 // if scale is 0., then use default values. Currently, that 01190 // default is continuous time (stored as a double) is seconds since 01191 // 1970-01-01 while broken-down time is Gregorian with no other 01192 // additional corrections. 01193 ( *qsasconfig )->scale = 1. / (double) SecInDay; 01194 ( *qsasconfig )->offset1 = (double) MJD_1970; 01195 ( *qsasconfig )->offset2 = 0.; 01196 ( *qsasconfig )->ccontrol = 0x0; 01197 } 01198 } 01199 01200 void closeqsas( QSASConfig **qsasconfig ) 01201 { 01202 // Close library if it has been opened. 01203 if ( *qsasconfig != NULL ) 01204 { 01205 free( (void *) *qsasconfig ); 01206 *qsasconfig = NULL; 01207 } 01208 } 01209 01210 int ctimeqsas( int year, int month, int day, int hour, int min, double sec, double * ctime, QSASConfig *qsasconfig ) 01211 { 01212 MJDtime MJD_value, *MJD = &MJD_value; 01213 int forceJulian, ret; 01214 01215 if ( qsasconfig == NULL ) 01216 { 01217 fprintf( stderr, "libqsastime (ctimeqsas) ERROR: configqsas must be called first.\n" ); 01218 exit( EXIT_FAILURE ); 01219 } 01220 01221 if ( qsasconfig->ccontrol & 0x1 ) 01222 forceJulian = 1; 01223 else 01224 forceJulian = 0; 01225 01226 ret = setFromUT( year, month, day, hour, min, sec, MJD, forceJulian ); 01227 if ( ret ) 01228 return ret; 01229 *ctime = ( ( (double) ( MJD->base_day ) - qsasconfig->offset1 ) - qsasconfig->offset2 + MJD->time_sec / (double) SecInDay ) / qsasconfig->scale; 01230 return 0; 01231 } 01232 01233 void btimeqsas( int *year, int *month, int *day, int *hour, int *min, double *sec, double ctime, QSASConfig *qsasconfig ) 01234 { 01235 MJDtime MJD_value, *MJD = &MJD_value; 01236 int forceJulian; 01237 double integral_offset1, integral_offset2, integral_scaled_ctime; 01238 int inleap; 01239 01240 if ( qsasconfig == NULL ) 01241 { 01242 fprintf( stderr, "libqsastime (btimeqsas) ERROR: configqsas must be called first.\n" ); 01243 exit( EXIT_FAILURE ); 01244 } 01245 01246 MJD->time_sec = SecInDay * ( modf( qsasconfig->offset1, &integral_offset1 ) + modf( qsasconfig->offset2, &integral_offset2 ) + modf( ctime * qsasconfig->scale, &integral_scaled_ctime ) ); 01247 MJD->base_day = (int) ( integral_offset1 + integral_offset2 + integral_scaled_ctime ); 01248 01249 if ( qsasconfig->ccontrol & 0x1 ) 01250 forceJulian = 1; 01251 else 01252 forceJulian = 0; 01253 01254 if ( qsasconfig->ccontrol & 0x2 ) 01255 MJD->time_sec += leap_second_TAI( MJD, &inleap, &( qsasconfig->index ) ); 01256 else 01257 inleap = 0; 01258 01259 // If in the middle of a positive leap increment insertion, normalize the 01260 // broken-down result so that *sec exceeds 60 to mark the insertion 01261 // (similar to the way February 29 marks a leap day). 01262 01263 if ( inleap ) 01264 MJD->time_sec -= 1.; 01265 01266 breakDownMJD( year, month, day, hour, min, sec, MJD, forceJulian ); 01267 if ( inleap ) 01268 *sec += 1.; 01269 } 01270 01271 size_t strfqsas( char * buf, size_t len, const char *format, double ctime, QSASConfig *qsasconfig ) 01272 { 01273 MJDtime MJD_value, *MJD = &MJD_value; 01274 int forceJulian; 01275 double integral_offset1, integral_offset2, integral_scaled_ctime; 01276 int inleap; 01277 01278 if ( qsasconfig == NULL ) 01279 { 01280 fprintf( stderr, "libqsastime (strfqsas) ERROR: configqsas must be called first.\n" ); 01281 exit( EXIT_FAILURE ); 01282 } 01283 MJD->time_sec = SecInDay * ( modf( qsasconfig->offset1, &integral_offset1 ) + modf( qsasconfig->offset2, &integral_offset2 ) + modf( ctime * qsasconfig->scale, &integral_scaled_ctime ) ); 01284 MJD->base_day = (int) ( integral_offset1 + integral_offset2 + integral_scaled_ctime ); 01285 01286 if ( qsasconfig->ccontrol & 0x1 ) 01287 forceJulian = 1; 01288 else 01289 forceJulian = 0; 01290 01291 if ( qsasconfig->ccontrol & 0x2 ) 01292 MJD->time_sec += leap_second_TAI( MJD, &inleap, &( qsasconfig->index ) ); 01293 else 01294 inleap = 0; 01295 01296 return strfMJD( buf, len, format, MJD, forceJulian, inleap ); 01297 } 01298 01299 // bhunt_search. Search an ordered table with a binary hunt phase and a 01300 // binary search phase. 01301 // 01302 // On entry *low is used to help the hunt phase speed up the binary 01303 // search when consecutive calls to bhunt_search are made with 01304 // similar key values. On exit, *low is adjusted such that 01305 // base[*low] <= key < base[(*low+1)] with the special cases of 01306 //low set to -1 to indicate the key < base[0] and *low set to n-1 01307 // to indicate base[n-1] <= key. The function *ge must return true (1) 01308 // if its first argument (the search key) is greater than or equal to 01309 // its second argument (a table entry). Otherwise it returns false 01310 // (0). Items in the array base must be in ascending order. 01311 01312 void bhunt_search( const void *key, const void *base, int n, size_t size, int *low, int ( *ge )( const void *keyval, const void *datum ) ) 01313 { 01314 const void *indexbase; 01315 int mid, high, hunt_inc = 1; 01316 // If previous search found below range, then assure one hunt cycle 01317 // just in case new key is also below range. 01318 if ( *low == -1 ) 01319 *low = 0; 01320 // Protect against invalid or undefined *low values where hunt 01321 // is waste of time. 01322 if ( *low < 0 || *low >= n ) 01323 { 01324 *low = -1; 01325 high = n; 01326 } 01327 else 01328 { 01329 // binary hunt phase where we are assured 0 <= *low < n 01330 indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) ( *low ) ) ); 01331 if ( ( *ge )( key, indexbase ) ) 01332 { 01333 high = ( *low ) + hunt_inc; 01334 indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) high ) ); 01335 // indexbase is valid if high < n. 01336 while ( ( high < n ) && ( ( *ge )( key, indexbase ) ) ) 01337 { 01338 *low = high; 01339 hunt_inc += hunt_inc; 01340 high = high + hunt_inc; 01341 indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) high ) ); 01342 } 01343 if ( high >= n ) 01344 high = n; 01345 // At this point, low is valid and base[low] <= key 01346 // and either key < base[high] for valid high or high = n. 01347 } 01348 else 01349 { 01350 high = *low; 01351 *low = high - hunt_inc; 01352 indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) ( *low ) ) ); 01353 // indexbase is valid if(*low) >= 0 01354 while ( ( ( *low ) >= 0 ) && !( ( *ge )( key, indexbase ) ) ) 01355 { 01356 high = *low; 01357 hunt_inc += hunt_inc; 01358 *low = ( *low ) - hunt_inc; 01359 indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) ( *low ) ) ); 01360 } 01361 if ( ( *low ) < 0 ) 01362 *low = -1; 01363 // At this point high is valid and key < base[high] 01364 // and either base[low] <= key for valid low or low = -1. 01365 } 01366 } 01367 // binary search phase where we are assured base[low] <= key < base[high] 01368 // when both low and high are valid with obvious special cases signalled 01369 // by low = -1 or high = n. 01370 while ( high - *low > 1 ) 01371 { 01372 mid = *low + ( high - *low ) / 2; 01373 indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) mid ) ); 01374 if ( ( *ge )( key, indexbase ) ) 01375 *low = mid; 01376 else 01377 high = mid; 01378 } 01379 }