// // This software originally contributed under the LGPL in January 2009 to // PLplot by the // Cluster Science Centre // QSAS team, // Imperial College, London // Copyright (C) 2009 Imperial College, London // Copyright (C) 2009 Alan W. Irwin // // This file is part of PLplot. // // PLplot is free software; you can redistribute it and/or modify // it under the terms of the GNU Library General Public License as published // by the Free Software Foundation; either version 2 of the License, or // (at your option) any later version. // // PLplot is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU Library General Public License for more details. // // You should have received a copy of the GNU Library General Public License // along with PLplot; if not, write to the Free Software // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA // // // MJD measures from the start of 17 Nov 1858 // These utilities use the Gregorian calendar after 4 Oct 1582 (Julian) i.e. from 15 Oct 1582 Gregorian // Note C libraries use Gregorian only from 14 Sept 1752 // More detailed discussion can be found at http://aa.usno.navy.mil/data/docs/JulianDate.php // These routines have been compared with the results of the US Naval Observatory online converter. // Modified Julian Date (MJD) = Julian Date (JD) - 2400000.5 // // In all routines, specifying a day, hour, minute or second field greater than would be valid is // handled with modulo arithmetic and safe. // Thus 2006-12-32 00:62:00.0 will safely, and correctly, be treated as 2007-01-01 01:02:00.0 // // #include #include #include "qsastimeP.h" #include "tai-utc.h" // MJD for 0000-01-01 (correctly Jan 01, BCE 1) // Julian proleptic calendar value. #define MJD_0000J -678943 // Gregorian proleptic calendar value. (At MJD_0000J the Gregorian proleptic // calendar reads two days behind the Julian proleptic calendar, i.e. - 2 days, // see http://en.wikipedia.org/wiki/Proleptic_Gregorian_calendar, // so MJD_0000G = MJD_0000J+2) #define MJD_0000G -678941 // MJD for 0001-01-01 which is 366 days later than previous definitions because // the year 0 is a leap year in both calendars. #define MJD_0001J -678577 #define MJD_0001G -678575 // MJD for Jan 01, 1970 00:00:00 Gregorian, the Unix epoch. #define MJD_1970 40587 static const double SecInDay = 86400; // we ignore leap seconds static const int MonthStartDOY[] = { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 }; static const int MonthStartDOY_L[] = { 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 }; // Static function declarations. static int geMJDtime_TAI( const MJDtime *number1, const TAI_UTC *number2 ); static int geMJDtime_UTC( const MJDtime *number1, const TAI_UTC *number2 ); static double leap_second_TAI( const MJDtime *MJD_TAI, int *inleap, int *index ); // End of static function declarations. int setFromUT( int year, int month, int day, int hour, int min, double sec, MJDtime *MJD, int forceJulian ) { // convert broken-down time to MJD // MJD measures from the start of 17 Nov 1858 // If forceJulian is true (non-zero), the Julian proleptic calendar is // used whatever the year. Otherwise, the Gregorian proleptic calendar // is used whatever the year. // Note C libraries use Gregorian only (at least on Linux). In contrast, // the Linux (and Unix?) cal application uses Julian for earlier dates // and Gregorian from 14 Sept 1752 onwards. int leaps, year4, year100, year400; double dbase_day, non_leaps = 365.; //int dbase_day, non_leaps = 365; double time_sec, dextraDays; int extraDays; if ( month < 0 || month > 11 ) { fprintf( stderr, "setfromUT: invalid month value\n" ); exit( EXIT_FAILURE ); } // As year increases, year4/4 increments by 1 at // year = -7, -3, 1, 5, 9, etc. // As year increases, year100/100 increments by 1 at // year = -299, -199, -99, 1, 101, 201, 301, etc. // As year increases, year400/400 increments by 1 at // year = -1199, -799, -399, 1, 401, 801, 1201, etc. if ( year <= 0 ) { year4 = year - 4; year100 = year - 100; year400 = year - 400; } else { year4 = year - 1; year100 = year - 1; year400 = year - 1; } if ( forceJulian ) { // count leap years on proleptic Julian Calendar starting from MJD_0000J leaps = year4 / 4; if ( year % 4 == 0 ) dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000J; else dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000J; } else { // count leap years for proleptic Gregorian Calendar. // Algorithm below for 1858-11-17 (0 MJD) gives // leaps = 450 and hence dbase_day of 678941, so subtract that value // or add MJD_0000G (which is two days different from MJD_0000J, see // above). leaps = year4 / 4 - year100 / 100 + year400 / 400; // left to right associativity means the double value of // non_leaps propagate to make all calculations be // done in double precision without the potential of // integer overflow. The result should be a double which // stores the expected exact integer results of the // calculation with exact representation unless the // result is much larger than the integer overflow limit. if ( ( year % 4 == 0 && year % 100 != 0 ) || ( year % 4 == 0 && year % 400 == 0 ) ) dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000G; else dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000G; } time_sec = sec + ( (double) min + (double) hour * 60. ) * 60.; if ( time_sec >= SecInDay ) { dextraDays = ( time_sec / SecInDay ); // precaution against overflowing extraDays. if ( fabs( dextraDays ) > 2.e9 ) { return 3; } extraDays = (int) ( dextraDays ); dbase_day += extraDays; time_sec -= extraDays * SecInDay; } // precaution against overflowing MJD->base_day. if ( fabs( dbase_day ) > 2.e9 ) { return 4; } else { // The exact integer result should be represented exactly in the // double, dbase_day, and its absolute value should be less than // the integer overflow limit. So the cast to int should be // exact. MJD->base_day = (int) dbase_day; MJD->time_sec = time_sec; return 0; } } void getYAD( int *year, int *ifleapyear, int *doy, const MJDtime *MJD, int forceJulian ) { // Get year and day of year from normalized MJD int j, ifcorrect, year4, year100, year400; j = MJD->base_day; if ( forceJulian ) { // Shift j epoch to 0000-01-01 for the Julian proleptic calendar. j -= MJD_0000J; // 365.25 is the exact period of the Julian year so year will be correct // if the day offset is set exactly right so that years -4, 0, 4 are // leap years, i.e. years -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 start with // j = -1826 -1461, -1095, -730, -365, 0, 366, 731, 1096, 1461, 1827 if ( j >= 366 ) { *year = (int) ( (double) ( j ) / 365.25 ); year4 = *year - 1; } else { *year = (int) ( (double) ( j - 365 ) / 365.25 ); year4 = *year - 4; } *doy = j - *year * 365 - year4 / 4; *ifleapyear = *year % 4 == 0; } else { // Shift j epoch to 0000-01-01 for the Gregorian proleptic calendar. j -= MJD_0000G; // 365.245 is the exact period of the Gregorian year so year will be correct // on average, but because the leap year rule is irregular within // the 400-year Gregorian cycle, the first and last days of the // year may need further adjustment. if ( j >= 366 ) { *year = (int) ( (double) ( j ) / 365.2425 ); year4 = *year - 1; year100 = *year - 1; year400 = *year - 1; } else { *year = (int) ( (double) ( j - 365 ) / 365.2425 ); year4 = *year - 4; year100 = *year - 100; year400 = *year - 400; } *doy = j - *year * 365 - year4 / 4 + year100 / 100 - year400 / 400; *ifleapyear = ( *year % 4 == 0 && *year % 100 != 0 ) || ( *year % 4 == 0 && *year % 400 == 0 ); // Rare corrections to above average Gregorian relations. if ( *doy < 1 ) { ( *year )--; ifcorrect = 1; } else if ( *doy > 365 && ( !*ifleapyear || *doy > 366 ) ) { ( *year )++; ifcorrect = 1; } else { ifcorrect = 0; } if ( ifcorrect ) { if ( j >= 366 ) { year4 = *year - 1; year100 = *year - 1; year400 = *year - 1; } else { year4 = *year - 4; year100 = *year - 100; year400 = *year - 400; } *doy = j - *year * 365 - year4 / 4 + year100 / 100 - year400 / 400; *ifleapyear = ( *year % 4 == 0 && *year % 100 != 0 ) || ( *year % 4 == 0 && *year % 400 == 0 ); } } } void normalize_MJD( MJDtime *MJD ) { int extra_days; // Calculate MJDout as normalized version // (i.e., 0. <= MJDout->time_sec < 86400.) of MJDin. if ( MJD->time_sec >= 0 ) { extra_days = (int) ( MJD->time_sec / SecInDay ); } else { // allow for negative seconds push into previous day even if less than 1 day extra_days = (int) ( MJD->time_sec / SecInDay ) - 1; } MJD->base_day += extra_days; MJD->time_sec -= extra_days * SecInDay; } void breakDownMJD( int *year, int *month, int *day, int *hour, int *min, double *sec, const MJDtime *MJD, int forceJulian ) { // Convert MJD struct into date/time elements // Note year 0 CE (AD) [1 BCE (BC)] is a leap year int doy, ifleapyear; MJDtime nMJD_value, *nMJD = &nMJD_value; *nMJD = *MJD; normalize_MJD( nMJD ); // Time part *sec = nMJD->time_sec; *hour = (int) ( *sec / 3600. ); *sec -= (double) *hour * 3600.; *min = (int) ( *sec / 60. ); *sec -= (double) *min * 60.; getYAD( year, &ifleapyear, &doy, nMJD, forceJulian ); // calculate month part with doy set to be the day within // the year in the range from 1 to 366 *month = -1; if ( ifleapyear ) { while ( doy > MonthStartDOY_L[*month + 1] ) { ( *month )++; if ( *month == 11 ) break; } *day = doy - MonthStartDOY_L[*month]; } else { while ( doy > MonthStartDOY[*month + 1] ) { ( *month )++; if ( *month == 11 ) break; } *day = doy - MonthStartDOY[*month]; } } const char * getDayOfWeek( const MJDtime *MJD ) { static const char *dow = { "Wed\0Thu\0Fri\0Sat\0Sun\0Mon\0Tue" }; int d = MJD->base_day % 7; if ( d < 0 ) d += 7; return &( dow[d * 4] ); } const char * getLongDayOfWeek( const MJDtime *MJD ) { static const char *dow = { "Wednesday\0Thursday\0\0Friday\0\0\0\0Saturday\0\0Sunday\0\0\0\0Monday\0\0\0\0Tuesday" }; int d = MJD->base_day % 7; if ( d < 0 ) d += 7; return &( dow[d * 10] ); } const char * getMonth( int m ) { static const char *months = { "Jan\0Feb\0Mar\0Apr\0May\0Jun\0Jul\0Aug\0Sep\0Oct\0Nov\0Dec" }; return &( months[( m ) * 4] ); } const char * getLongMonth( int m ) { 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" }; return &( months[( m ) * 10] ); } size_t strfMJD( char * buf, size_t len, const char *format, const MJDtime *MJD, int forceJulian, int inleap ) { // Format a text string according to the format string. // Uses the same syntax as strftime() but does not use current locale. // The null terminator is included in len for safety. // if inleap is true (non-zero) then renormalize so that (int) sec // is 60 to mark results as a flag that a leap increment (recently // always a second, but historically it was sometimes smaller than // that) was in the process of being inserted for this particular // epoch just prior to a positive discontinuity in TAI-UTC. int year, month, day, hour, min, ysign, second, d, y; int y1, ifleapyear; int i, secsSince1970; int nplaces, fmtlen, slen; int resolution; double shiftPlaces; char * ptr; double sec, sec_fraction; int w, doy, days_in_wk1; const char *dayText; const char *monthText; char DateTime[80]; size_t posn = 0; size_t last = len - 1; MJDtime nMJD_value, *nMJD = &nMJD_value; char dynamic_format[10]; // Find required resolution resolution = 0; fmtlen = (int) strlen( format ); i = 0; while ( i < fmtlen ) { char next = format[i]; if ( next == '%' ) { // find seconds format if used i++; next = format[i]; if ( isdigit( next ) != 0 ) { nplaces = (int) strtol( &( format[i] ), NULL, 10 ); if ( nplaces > resolution ) resolution = nplaces; } else if ( next == '.' ) { resolution = 9; // maximum resolution allowed } } i++; } // ensure rounding is done before breakdown shiftPlaces = pow( 10, (double) resolution ); *nMJD = *MJD; nMJD->time_sec += 0.5 / shiftPlaces; buf[last] = '\0'; buf[0] = '\0'; // force overwrite of old buffer since strnctat() used hereafter if ( inleap ) nMJD->time_sec -= 1.; breakDownMJD( &year, &month, &day, &hour, &min, &sec, nMJD, forceJulian ); if ( inleap ) sec += 1.; if ( year < 0 ) { ysign = 1; year = -year; } else ysign = 0; //truncate seconds to resolution to stop formatting rounding up sec = floor( sec * shiftPlaces ) / shiftPlaces; second = (int) sec; // Read format string, character at a time i = 0; while ( i < fmtlen ) { char next = format[i]; if ( next == '%' ) { // format character or escape i++; next = format[i]; if ( next == '%' ) { // escaped %, pass it on buf[posn] = next; posn++; if ( posn >= last ) return posn; } else if ( next == 'a' ) { // short day name dayText = getDayOfWeek( nMJD ); strncat( &( buf[posn] ), dayText, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'A' ) { // long day name dayText = getLongDayOfWeek( nMJD ); strncat( &( buf[posn] ), dayText, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'b' || next == 'h' ) { // short month name monthText = getMonth( month ); strncat( &( buf[posn] ), monthText, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'B' ) { // long month name monthText = getLongMonth( month ); strncat( &( buf[posn] ), monthText, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'c' ) { // Date and Time with day of week dayText = getDayOfWeek( nMJD ); monthText = getMonth( month ); if ( ysign == 0 ) sprintf( DateTime, "%s %s %02d %02d:%02d:%02d %04d", dayText, monthText, day, hour, min, second, year ); else sprintf( DateTime, "%s %s %02d %02d:%02d:%02d -%04d", dayText, monthText, day, hour, min, second, year ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'C' ) { // year / 100 so, e.g. 1989 is 20th century but comes out as 19 int century = year / 100; if ( ysign == 0 ) sprintf( DateTime, "%02d", century ); else sprintf( DateTime, "-%02d", century + 1 ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'd' ) { // day of month (01 - 31) sprintf( DateTime, "%02d", day ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'D' ) { // month/day/year y = year % 100; if ( ysign == 0 ) sprintf( DateTime, "%02d/%02d/%02d", month + 1, day, y ); else sprintf( DateTime, "%02d/%02d/-%02d", month + 1, day, y ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'e' ) { // day of month ( 1 - 31) if ( day < 10 ) sprintf( DateTime, " %01d", day ); else sprintf( DateTime, "%02d", day ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'F' ) { // year-month-day if ( ysign == 0 ) sprintf( DateTime, "%04d-%02d-%02d", year, month + 1, day ); else sprintf( DateTime, "-%04d-%02d-%02d", year, month + 1, day ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'H' ) { // hour, 24 hour clock (00 - 23) sprintf( DateTime, "%02d", hour ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'I' ) { // hour, 12 hour clock (01 - 12) if ( hour == 0 ) sprintf( DateTime, "%02d", hour + 12 ); else if ( hour > 12 ) sprintf( DateTime, "%02d", hour - 12 ); else sprintf( DateTime, "%02d", hour ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'j' ) { // day of year getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian ); sprintf( DateTime, "%03d", doy ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'k' ) { // hour, 24 hour clock ( 0 - 23) if ( hour < 10 ) sprintf( DateTime, " %01d", hour ); else sprintf( DateTime, "%02d", hour ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'l' ) { // hour, 12 hour clock ( 1 - 12) if ( hour == 0 ) sprintf( DateTime, "%02d", hour + 12 ); else if ( hour < 10 ) sprintf( DateTime, " %01d", hour ); else if ( hour <= 12 ) sprintf( DateTime, "%02d", hour ); else if ( hour < 22 ) sprintf( DateTime, " %01d", hour - 12 ); else sprintf( DateTime, "%02d", hour - 12 ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'm' ) { // month (01 - 12) sprintf( DateTime, "%02d", month + 1 ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'M' ) { // minute (00 - 59) sprintf( DateTime, "%02d", min ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'n' ) { // newline buf[posn] = '\n'; posn++; if ( posn >= last ) return posn; } else if ( next == 'p' ) { // am/pm on12 hour clock if ( hour < 0 ) sprintf( DateTime, "AM" ); else sprintf( DateTime, "PM" ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'r' ) { // hour:min:sec AM, 12 hour clock (01 - 12):(00 - 59):(00 - 59) (AM - PM) if ( hour == 0 ) sprintf( DateTime, "%02d:%02d:%02d AM", hour + 12, min, second ); else if ( hour > 12 ) sprintf( DateTime, "%02d:%02d:%02d PM", hour - 12, min, second ); else if ( hour == 12 ) sprintf( DateTime, "%02d:%02d:%02d PM", hour, min, second ); else sprintf( DateTime, "%02d:%02d:%02d AM", hour, min, second ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'R' ) { // hour:min, 24 hour clock (00 - 23):(00 - 59) sprintf( DateTime, "%02d:%02d", hour, min ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'S' ) { // second (00 - 59 with optional decimal point and numbers after the decimal point.) if ( i + 2 < fmtlen && format[i + 1] == '%' && ( format[i + 2] == '.' || isdigit( format[i + 2] ) != 0 ) ) { // nplaces is number of decimal places ( 0 < nplaces <= 9 ) if ( format[i + 2] == '.' ) // maximum number of places nplaces = 9; else nplaces = (int) strtol( &( format[i + 2] ), NULL, 10 ); i += 2; } else { nplaces = 0; } if ( nplaces == 0 ) { sprintf( DateTime, "%02d", (int) ( sec + 0.5 ) ); } else { sprintf( dynamic_format, "%%0%d.%df", nplaces + 3, nplaces ); sprintf( DateTime, dynamic_format, sec ); if ( format[i] == '.' ) { slen = (int) strlen( DateTime ) - 1; while ( DateTime[slen] == '0' && DateTime[slen - 1] != '.' ) { DateTime[slen] = '\0'; // remove trailing zeros slen--; } } } strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 's' ) { // seconds since 01 Jan 1970 Gregorian secsSince1970 = (int) ( nMJD->time_sec + ( nMJD->base_day - MJD_1970 ) * SecInDay ); sprintf( DateTime, "%d", secsSince1970 ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 't' ) { // tab buf[posn] = '\t'; posn++; if ( posn >= last ) return posn; } else if ( next == 'T' ) { // hour:min:sec, 24 hour clock (00 - 23):(00 - 59):(00 - 59) sprintf( DateTime, "%02d:%02d:%02d", hour, min, second ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'U' ) { // week of year as a number, (00 - 53) start of week is Sunday getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian ); days_in_wk1 = ( nMJD->base_day - doy - 4 ) % 7; w = ( doy + 6 - days_in_wk1 ) / 7; sprintf( DateTime, "%02d", w ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'u' ) { // weekday as a number, 0 = Monday d = 1 + ( nMJD->base_day - 5 ) % 7; sprintf( DateTime, "%01d", d ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'v' ) { // day-MonthName-year day of month ( 1 - 31) - month (Jan ... Dec) - year (yyyy) monthText = getMonth( month ); if ( ysign == 0 ) { if ( day < 10 ) sprintf( DateTime, " %01d-%s-%04d", day, monthText, year ); else sprintf( DateTime, "%02d-%s-%04d", day, monthText, year ); } else { if ( day < 10 ) sprintf( DateTime, " %01d-%s-(-)%04d", day, monthText, year ); else sprintf( DateTime, "%02d-%s-(-)%04d", day, monthText, year ); } strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'V' ) { // week of year as a number, (01 - 53) start of week is Monday and first week has at least 3 days in year getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian ); days_in_wk1 = ( nMJD->base_day - doy - 3 ) % 7; if ( days_in_wk1 <= 3 ) w = ( doy + 6 - days_in_wk1 ) / 7; // ensure first week has at least 3 days in this year else w = 1 + ( doy + 6 - days_in_wk1 ) / 7; if ( w == 0 ) w = 53; sprintf( DateTime, "%02d", w ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'w' ) { // weekday as a number, 0 = Sunday d = ( nMJD->base_day - 4 ) % 7; sprintf( DateTime, "%01d", d ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'W' ) { // week of year as a number, (00 - 53) start of week is Monday getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian ); days_in_wk1 = ( nMJD->base_day - doy - 3 ) % 7; w = ( doy + 6 - days_in_wk1 ) / 7; sprintf( DateTime, "%02d", w ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'x' ) { // date string dayText = getDayOfWeek( nMJD ); monthText = getMonth( month ); if ( ysign == 0 ) sprintf( DateTime, "%s %s %02d, %04d", dayText, monthText, day, year ); else sprintf( DateTime, "%s %s %02d, -%04d", dayText, monthText, day, year ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'X' ) { // time string sprintf( DateTime, "%02d:%02d:%02d", hour, min, second ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'y' ) { // 2 digit year y = year % 100; if ( ysign == 0 ) sprintf( DateTime, "%02d", y ); else sprintf( DateTime, "-%02d", y ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'Y' ) { // 4 digit year if ( ysign == 0 ) sprintf( DateTime, "%04d", year ); else sprintf( DateTime, "-%04d", year ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'Z' ) { // time zone and calendar, always UTC if ( forceJulian ) strncat( &( buf[posn] ), "UTC Julian", last - posn ); else strncat( &( buf[posn] ), "UTC Gregorian", last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == 'z' ) { // time zone, always UTC strncat( &( buf[posn] ), "+0000", last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == '+' ) { // date and time dayText = getDayOfWeek( nMJD ); monthText = getMonth( month ); if ( ysign == 0 ) sprintf( DateTime, "%s %s %02d %02d:%02d:%02d UTC %04d", dayText, monthText, day, hour, min, second, year ); else sprintf( DateTime, "%s %s %02d %02d:%02d:%02d UTC -%04d", dayText, monthText, day, hour, min, second, year ); strncat( &( buf[posn] ), DateTime, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } else if ( next == '.' || isdigit( next ) != 0 ) { // nplaces is number of decimal places ( 0 < nplaces <= 9 ) if ( next == '.' ) // maximum number of places nplaces = 9; else nplaces = (int) strtol( &( format[i] ), NULL, 10 ); // fractional part of seconds to maximum available accuracy sec_fraction = sec - (int) sec; sprintf( dynamic_format, "%%-%d.%df", nplaces + 2, nplaces ); //sprintf(DateTime, "%-11.9f", sec_fraction); sprintf( DateTime, dynamic_format, sec_fraction ); while ( ( ptr = strrchr( &( DateTime[0] ), ' ' ) ) != NULL ) ptr[0] = '\0'; // remove trailing white space if ( next == '.' ) { slen = (int) strlen( DateTime ) - 1; while ( DateTime[slen] == '0' && DateTime[slen - 1] != '.' ) { DateTime[slen] = '\0'; // remove trailing zeros slen--; } } ptr = strchr( DateTime, '.' ); // remove everything in front of the decimal point, and // ignore case (%0) where no decimal point exists at all. if ( ptr != NULL ) strncat( &( buf[posn] ), ptr, last - posn ); posn = strlen( buf ); if ( posn >= last ) return posn; } } else { // regular multi-byte character, pass it on buf[posn] = next; posn++; if ( posn >= last ) return posn; } buf[posn] = '\0'; i++; } return posn; } int geMJDtime_TAI( const MJDtime *number1, const TAI_UTC *number2 ) { // Returns true if number1 >= number2. // N.B. both number1 and number2 must be normalized. if ( number1->base_day > number2->base_day ) { return 1; } else if ( number1->base_day < number2->base_day ) { return 0; } else { return ( number1->time_sec >= number2->time_sec_tai ); } } int geMJDtime_UTC( const MJDtime *number1, const TAI_UTC *number2 ) { // Returns true if number1 >= number2. // N.B. both number1 and number2 must be normalized. if ( number1->base_day > number2->base_day ) { return 1; } else if ( number1->base_day < number2->base_day ) { return 0; } else { return ( number1->time_sec >= number2->time_sec_utc ); } } double leap_second_TAI( const MJDtime *MJD_TAI, int *inleap, int *index ) { // Logic assumes input MJD_TAI is in TAI // *inleap lets the calling routine know whether MJD_TAI corresponds // to an epoch when a positive leap increment is being inserted. MJDtime MJD_value, *MJD = &MJD_value; double leap; int debug = 0; // N.B. geMJDtime_TAI only works for normalized values. *MJD = *MJD_TAI; normalize_MJD( MJD ); // Search for index such that TAI_UTC_lookup_table[*index] <= MJD(TAI) < TAI_UTC_lookup_table[*index+1] 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 ); if ( debug == 2 ) fprintf( stderr, "*index = %d\n", *index ); if ( *index == -1 ) { // MJD is less than first table entry. // Debug: check that condition is met if ( debug && geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index + 1] ) ) { fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index ); exit( EXIT_FAILURE ); } // There is (by assertion) no discontinuity at the start of the table. // Therefore, *inleap cannot be true. *inleap = 0; // Use initial offset for MJD values before first table entry. // Calculate this offset strictly from offset1. The slope term // doesn't enter because offset2 is the same as the UTC of the // first epoch of the table. return -TAI_UTC_lookup_table[*index + 1].offset1; } else if ( *index == number_of_entries_in_tai_utc_table - 1 ) { // MJD is greater than or equal to last table entry. // Debug: check that condition is met if ( debug && !geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index] ) ) { fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index ); exit( EXIT_FAILURE ); } // If beyond end of table, cannot be in middle of leap second insertion. *inleap = 0; // Use final offset for MJD values after last table entry. // The slope term doesn't enter because modern values of the slope // are zero. return -TAI_UTC_lookup_table[*index].offset1; } else if ( *index >= 0 && *index < number_of_entries_in_tai_utc_table ) { // table[*index] <= MJD < table[*index+1]. // Debug: check that condition is met if ( debug && !( geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index] ) && !geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index + 1] ) ) ) { fprintf( stderr, "MJD = {%d, %f}\n", MJD->base_day, MJD->time_sec ); fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index ); exit( EXIT_FAILURE ); } 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 ); // Convert MJD(TAI) to normalized MJD(UTC). MJD->time_sec += leap; normalize_MJD( MJD ); // If MJD(UT) is in the next interval of the corresponding // TAI_UTC_lookup_table, then we are right in the middle of a // leap interval (recently a second but for earlier epochs it could be // less) insertion. Note this logic even works when leap intervals // are taken away from UTC (i.e., leap is positive) since in that // case the UTC *index always corresponds to the TAI *index. *inleap = geMJDtime_UTC( MJD, &TAI_UTC_lookup_table[*index + 1] ); return leap; } else { fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad *index = %d\n", *index ); exit( EXIT_FAILURE ); } } 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 ) { // Configure the transformation between continuous time and broken-down time // that is used for ctimeqsas, btimeqsas, and strfqsas. int forceJulian, ret; MJDtime MJD_value, *MJD = &MJD_value; // Allocate memory for *qsasconfig if that hasn't been done by a // previous call. if ( *qsasconfig == NULL ) { *qsasconfig = (QSASConfig *) malloc( (size_t) sizeof ( QSASConfig ) ); if ( *qsasconfig == NULL ) { fprintf( stderr, "configqsas: out of memory\n" ); exit( EXIT_FAILURE ); } } // Set bhunt_search index to a definite value less than -1 to insure no // initial hunt phase from some random index value is done. ( *qsasconfig )->index = -40; if ( scale != 0. ) { if ( ifbtime_offset ) { if ( ccontrol & 0x1 ) forceJulian = 1; else forceJulian = 0; ret = setFromUT( year, month, day, hour, min, sec, MJD, forceJulian ); if ( ret ) { fprintf( stderr, "configqsas: some problem with broken-down arguments\n" ); exit( EXIT_FAILURE ); } offset1 = (double) MJD->base_day; offset2 = MJD->time_sec / (double) SecInDay; } ( *qsasconfig )->scale = scale; ( *qsasconfig )->offset1 = offset1; ( *qsasconfig )->offset2 = offset2; ( *qsasconfig )->ccontrol = ccontrol; } else { // if scale is 0., then use default values. Currently, that // default is continuous time (stored as a double) is seconds since // 1970-01-01 while broken-down time is Gregorian with no other // additional corrections. ( *qsasconfig )->scale = 1. / (double) SecInDay; ( *qsasconfig )->offset1 = (double) MJD_1970; ( *qsasconfig )->offset2 = 0.; ( *qsasconfig )->ccontrol = 0x0; } } void closeqsas( QSASConfig **qsasconfig ) { // Close library if it has been opened. if ( *qsasconfig != NULL ) { free( (void *) *qsasconfig ); *qsasconfig = NULL; } } int ctimeqsas( int year, int month, int day, int hour, int min, double sec, double * ctime, QSASConfig *qsasconfig ) { MJDtime MJD_value, *MJD = &MJD_value; int forceJulian, ret; if ( qsasconfig == NULL ) { fprintf( stderr, "libqsastime (ctimeqsas) ERROR: configqsas must be called first.\n" ); exit( EXIT_FAILURE ); } if ( qsasconfig->ccontrol & 0x1 ) forceJulian = 1; else forceJulian = 0; ret = setFromUT( year, month, day, hour, min, sec, MJD, forceJulian ); if ( ret ) return ret; *ctime = ( ( (double) ( MJD->base_day ) - qsasconfig->offset1 ) - qsasconfig->offset2 + MJD->time_sec / (double) SecInDay ) / qsasconfig->scale; return 0; } void btimeqsas( int *year, int *month, int *day, int *hour, int *min, double *sec, double ctime, QSASConfig *qsasconfig ) { MJDtime MJD_value, *MJD = &MJD_value; int forceJulian; double integral_offset1, integral_offset2, integral_scaled_ctime; int inleap; if ( qsasconfig == NULL ) { fprintf( stderr, "libqsastime (btimeqsas) ERROR: configqsas must be called first.\n" ); exit( EXIT_FAILURE ); } MJD->time_sec = SecInDay * ( modf( qsasconfig->offset1, &integral_offset1 ) + modf( qsasconfig->offset2, &integral_offset2 ) + modf( ctime * qsasconfig->scale, &integral_scaled_ctime ) ); MJD->base_day = (int) ( integral_offset1 + integral_offset2 + integral_scaled_ctime ); if ( qsasconfig->ccontrol & 0x1 ) forceJulian = 1; else forceJulian = 0; if ( qsasconfig->ccontrol & 0x2 ) MJD->time_sec += leap_second_TAI( MJD, &inleap, &( qsasconfig->index ) ); else inleap = 0; // If in the middle of a positive leap increment insertion, normalize the // broken-down result so that *sec exceeds 60 to mark the insertion // (similar to the way February 29 marks a leap day). if ( inleap ) MJD->time_sec -= 1.; breakDownMJD( year, month, day, hour, min, sec, MJD, forceJulian ); if ( inleap ) *sec += 1.; } size_t strfqsas( char * buf, size_t len, const char *format, double ctime, QSASConfig *qsasconfig ) { MJDtime MJD_value, *MJD = &MJD_value; int forceJulian; double integral_offset1, integral_offset2, integral_scaled_ctime; int inleap; if ( qsasconfig == NULL ) { fprintf( stderr, "libqsastime (strfqsas) ERROR: configqsas must be called first.\n" ); exit( EXIT_FAILURE ); } MJD->time_sec = SecInDay * ( modf( qsasconfig->offset1, &integral_offset1 ) + modf( qsasconfig->offset2, &integral_offset2 ) + modf( ctime * qsasconfig->scale, &integral_scaled_ctime ) ); MJD->base_day = (int) ( integral_offset1 + integral_offset2 + integral_scaled_ctime ); if ( qsasconfig->ccontrol & 0x1 ) forceJulian = 1; else forceJulian = 0; if ( qsasconfig->ccontrol & 0x2 ) MJD->time_sec += leap_second_TAI( MJD, &inleap, &( qsasconfig->index ) ); else inleap = 0; return strfMJD( buf, len, format, MJD, forceJulian, inleap ); } // bhunt_search. Search an ordered table with a binary hunt phase and a // binary search phase. // // On entry *low is used to help the hunt phase speed up the binary // search when consecutive calls to bhunt_search are made with // similar key values. On exit, *low is adjusted such that // base[*low] <= key < base[(*low+1)] with the special cases of //low set to -1 to indicate the key < base[0] and *low set to n-1 // to indicate base[n-1] <= key. The function *ge must return true (1) // if its first argument (the search key) is greater than or equal to // its second argument (a table entry). Otherwise it returns false // (0). Items in the array base must be in ascending order. void bhunt_search( const void *key, const void *base, int n, size_t size, int *low, int ( *ge )( const void *keyval, const void *datum ) ) { const void *indexbase; int mid, high, hunt_inc = 1; // If previous search found below range, then assure one hunt cycle // just in case new key is also below range. if ( *low == -1 ) *low = 0; // Protect against invalid or undefined *low values where hunt // is waste of time. if ( *low < 0 || *low >= n ) { *low = -1; high = n; } else { // binary hunt phase where we are assured 0 <= *low < n indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) ( *low ) ) ); if ( ( *ge )( key, indexbase ) ) { high = ( *low ) + hunt_inc; indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) high ) ); // indexbase is valid if high < n. while ( ( high < n ) && ( ( *ge )( key, indexbase ) ) ) { *low = high; hunt_inc += hunt_inc; high = high + hunt_inc; indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) high ) ); } if ( high >= n ) high = n; // At this point, low is valid and base[low] <= key // and either key < base[high] for valid high or high = n. } else { high = *low; *low = high - hunt_inc; indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) ( *low ) ) ); // indexbase is valid if(*low) >= 0 while ( ( ( *low ) >= 0 ) && !( ( *ge )( key, indexbase ) ) ) { high = *low; hunt_inc += hunt_inc; *low = ( *low ) - hunt_inc; indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) ( *low ) ) ); } if ( ( *low ) < 0 ) *low = -1; // At this point high is valid and key < base[high] // and either base[low] <= key for valid low or low = -1. } } // binary search phase where we are assured base[low] <= key < base[high] // when both low and high are valid with obvious special cases signalled // by low = -1 or high = n. while ( high - *low > 1 ) { mid = *low + ( high - *low ) / 2; indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) mid ) ); if ( ( *ge )( key, indexbase ) ) *low = mid; else high = mid; } }