PLplot  5.10.0
qsastime.c
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines