PLplot  5.10.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
qsastime.c
Go to the documentation of this file.
1 //
2 // This software originally contributed under the LGPL in January 2009 to
3 // PLplot by the
4 // Cluster Science Centre
5 // QSAS team,
6 // Imperial College, London
7 // Copyright (C) 2009 Imperial College, London
8 // Copyright (C) 2009 Alan W. Irwin
9 //
10 // This file is part of PLplot.
11 //
12 // PLplot is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU Library General Public License as published
14 // by the Free Software Foundation; either version 2 of the License, or
15 // (at your option) any later version.
16 //
17 // PLplot is distributed in the hope that it will be useful,
18 // but WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 // GNU Library General Public License for more details.
21 //
22 // You should have received a copy of the GNU Library General Public License
23 // along with PLplot; if not, write to the Free Software
24 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25 //
26 //
27 
28 // MJD measures from the start of 17 Nov 1858
29 
30 // These utilities use the Gregorian calendar after 4 Oct 1582 (Julian) i.e. from 15 Oct 1582 Gregorian
31 // Note C libraries use Gregorian only from 14 Sept 1752
32 // More detailed discussion can be found at http://aa.usno.navy.mil/data/docs/JulianDate.php
33 // These routines have been compared with the results of the US Naval Observatory online converter.
34 // Modified Julian Date (MJD) = Julian Date (JD) - 2400000.5
35 //
36 // In all routines, specifying a day, hour, minute or second field greater than would be valid is
37 // handled with modulo arithmetic and safe.
38 // Thus 2006-12-32 00:62:00.0 will safely, and correctly, be treated as 2007-01-01 01:02:00.0
39 //
40 //
41 #include <ctype.h>
42 #include <math.h>
43 #include "qsastimeP.h"
44 #include "tai-utc.h"
45 
46 // MJD for 0000-01-01 (correctly Jan 01, BCE 1)
47 // Julian proleptic calendar value.
48 #define MJD_0000J -678943
49 // Gregorian proleptic calendar value. (At MJD_0000J the Gregorian proleptic
50 // calendar reads two days behind the Julian proleptic calendar, i.e. - 2 days,
51 // see http://en.wikipedia.org/wiki/Proleptic_Gregorian_calendar,
52 // so MJD_0000G = MJD_0000J+2)
53 #define MJD_0000G -678941
54 // MJD for 0001-01-01 which is 366 days later than previous definitions because
55 // the year 0 is a leap year in both calendars.
56 #define MJD_0001J -678577
57 #define MJD_0001G -678575
58 // MJD for Jan 01, 1970 00:00:00 Gregorian, the Unix epoch.
59 #define MJD_1970 40587
60 
61 static const double SecInDay = 86400; // we ignore leap seconds
62 static const int MonthStartDOY[] = { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 };
63 static const int MonthStartDOY_L[] = { 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 };
64 
65 // Static function declarations.
66 static int geMJDtime_TAI( const MJDtime *number1, const TAI_UTC *number2 );
67 static int geMJDtime_UTC( const MJDtime *number1, const TAI_UTC *number2 );
68 static double leap_second_TAI( const MJDtime *MJD_TAI, int *inleap, int *index );
69 // End of static function declarations.
70 
71 int setFromUT( int year, int month, int day, int hour, int min, double sec, MJDtime *MJD, int forceJulian )
72 {
73  // convert broken-down time to MJD
74  // MJD measures from the start of 17 Nov 1858
75 
76  // If forceJulian is true (non-zero), the Julian proleptic calendar is
77  // used whatever the year. Otherwise, the Gregorian proleptic calendar
78  // is used whatever the year.
79  // Note C libraries use Gregorian only (at least on Linux). In contrast,
80  // the Linux (and Unix?) cal application uses Julian for earlier dates
81  // and Gregorian from 14 Sept 1752 onwards.
82 
83  int leaps, year4, year100, year400;
84  double dbase_day, non_leaps = 365.;
85  //int dbase_day, non_leaps = 365;
86  double time_sec, dextraDays;
87  int extraDays;
88 
89  if ( month < 0 || month > 11 )
90  {
91  fprintf( stderr, "setfromUT: invalid month value\n" );
92  exit( EXIT_FAILURE );
93  }
94  // As year increases, year4/4 increments by 1 at
95  // year = -7, -3, 1, 5, 9, etc.
96  // As year increases, year100/100 increments by 1 at
97  // year = -299, -199, -99, 1, 101, 201, 301, etc.
98  // As year increases, year400/400 increments by 1 at
99  // year = -1199, -799, -399, 1, 401, 801, 1201, etc.
100  if ( year <= 0 )
101  {
102  year4 = year - 4;
103  year100 = year - 100;
104  year400 = year - 400;
105  }
106  else
107  {
108  year4 = year - 1;
109  year100 = year - 1;
110  year400 = year - 1;
111  }
112 
113  if ( forceJulian )
114  {
115  // count leap years on proleptic Julian Calendar starting from MJD_0000J
116  leaps = year4 / 4;
117  if ( year % 4 == 0 )
118  dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000J;
119  else
120  dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000J;
121  }
122  else
123  {
124  // count leap years for proleptic Gregorian Calendar.
125  // Algorithm below for 1858-11-17 (0 MJD) gives
126  // leaps = 450 and hence dbase_day of 678941, so subtract that value
127  // or add MJD_0000G (which is two days different from MJD_0000J, see
128  // above).
129  leaps = year4 / 4 - year100 / 100 + year400 / 400;
130 
131  // left to right associativity means the double value of
132  // non_leaps propagate to make all calculations be
133  // done in double precision without the potential of
134  // integer overflow. The result should be a double which
135  // stores the expected exact integer results of the
136  // calculation with exact representation unless the
137  // result is much larger than the integer overflow limit.
138  if ( ( year % 4 == 0 && year % 100 != 0 ) || ( year % 4 == 0 && year % 400 == 0 ) )
139  dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000G;
140  else
141  dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000G;
142  }
143 
144  time_sec = sec + ( (double) min + (double) hour * 60. ) * 60.;
145 
146  if ( time_sec >= SecInDay )
147  {
148  dextraDays = ( time_sec / SecInDay );
149  // precaution against overflowing extraDays.
150  if ( fabs( dextraDays ) > 2.e9 )
151  {
152  return 3;
153  }
154  extraDays = (int) ( dextraDays );
155  dbase_day += extraDays;
156  time_sec -= extraDays * SecInDay;
157  }
158  // precaution against overflowing MJD->base_day.
159  if ( fabs( dbase_day ) > 2.e9 )
160  {
161  return 4;
162  }
163  else
164  {
165  // The exact integer result should be represented exactly in the
166  // double, dbase_day, and its absolute value should be less than
167  // the integer overflow limit. So the cast to int should be
168  // exact.
169  MJD->base_day = (int) dbase_day;
170  MJD->time_sec = time_sec;
171  return 0;
172  }
173 }
174 
175 void getYAD( int *year, int *ifleapyear, int *doy, const MJDtime *MJD, int forceJulian )
176 {
177  // Get year and day of year from normalized MJD
178 
179  int j, ifcorrect, year4, year100, year400;
180 
181  j = MJD->base_day;
182 
183  if ( forceJulian )
184  {
185  // Shift j epoch to 0000-01-01 for the Julian proleptic calendar.
186  j -= MJD_0000J;
187 
188  // 365.25 is the exact period of the Julian year so year will be correct
189  // if the day offset is set exactly right so that years -4, 0, 4 are
190  // leap years, i.e. years -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 start with
191  // j = -1826 -1461, -1095, -730, -365, 0, 366, 731, 1096, 1461, 1827
192  if ( j >= 366 )
193  {
194  *year = (int) ( (double) ( j ) / 365.25 );
195  year4 = *year - 1;
196  }
197  else
198  {
199  *year = (int) ( (double) ( j - 365 ) / 365.25 );
200  year4 = *year - 4;
201  }
202 
203  *doy = j - *year * 365 - year4 / 4;
204 
205  *ifleapyear = *year % 4 == 0;
206  }
207  else
208  {
209  // Shift j epoch to 0000-01-01 for the Gregorian proleptic calendar.
210  j -= MJD_0000G;
211  // 365.245 is the exact period of the Gregorian year so year will be correct
212  // on average, but because the leap year rule is irregular within
213  // the 400-year Gregorian cycle, the first and last days of the
214  // year may need further adjustment.
215 
216  if ( j >= 366 )
217  {
218  *year = (int) ( (double) ( j ) / 365.2425 );
219  year4 = *year - 1;
220  year100 = *year - 1;
221  year400 = *year - 1;
222  }
223  else
224  {
225  *year = (int) ( (double) ( j - 365 ) / 365.2425 );
226  year4 = *year - 4;
227  year100 = *year - 100;
228  year400 = *year - 400;
229  }
230 
231  *doy = j - *year * 365 - year4 / 4 + year100 / 100 - year400 / 400;
232  *ifleapyear = ( *year % 4 == 0 && *year % 100 != 0 ) || ( *year % 4 == 0 && *year % 400 == 0 );
233 
234  // Rare corrections to above average Gregorian relations.
235  if ( *doy < 1 )
236  {
237  ( *year )--;
238  ifcorrect = 1;
239  }
240  else if ( *doy > 365 && ( !*ifleapyear || *doy > 366 ) )
241  {
242  ( *year )++;
243  ifcorrect = 1;
244  }
245  else
246  {
247  ifcorrect = 0;
248  }
249  if ( ifcorrect )
250  {
251  if ( j >= 366 )
252  {
253  year4 = *year - 1;
254  year100 = *year - 1;
255  year400 = *year - 1;
256  }
257  else
258  {
259  year4 = *year - 4;
260  year100 = *year - 100;
261  year400 = *year - 400;
262  }
263 
264  *doy = j - *year * 365 - year4 / 4 + year100 / 100 - year400 / 400;
265  *ifleapyear = ( *year % 4 == 0 && *year % 100 != 0 ) || ( *year % 4 == 0 && *year % 400 == 0 );
266  }
267  }
268 }
269 
270 void normalize_MJD( MJDtime *MJD )
271 {
272  int extra_days;
273  // Calculate MJDout as normalized version
274  // (i.e., 0. <= MJDout->time_sec < 86400.) of MJDin.
275  if ( MJD->time_sec >= 0 )
276  {
277  extra_days = (int) ( MJD->time_sec / SecInDay );
278  }
279  else
280  {
281  // allow for negative seconds push into previous day even if less than 1 day
282  extra_days = (int) ( MJD->time_sec / SecInDay ) - 1;
283  }
284 
285  MJD->base_day += extra_days;
286  MJD->time_sec -= extra_days * SecInDay;
287 }
288 
289 void breakDownMJD( int *year, int *month, int *day, int *hour, int *min, double *sec, const MJDtime *MJD, int forceJulian )
290 {
291  // Convert MJD struct into date/time elements
292  // Note year 0 CE (AD) [1 BCE (BC)] is a leap year
293 
294  int doy, ifleapyear;
295  MJDtime nMJD_value, *nMJD = &nMJD_value;
296 
297  *nMJD = *MJD;
298  normalize_MJD( nMJD );
299  // Time part
300 
301  *sec = nMJD->time_sec;
302  *hour = (int) ( *sec / 3600. );
303  *sec -= (double) *hour * 3600.;
304  *min = (int) ( *sec / 60. );
305  *sec -= (double) *min * 60.;
306 
307  getYAD( year, &ifleapyear, &doy, nMJD, forceJulian );
308 
309  // calculate month part with doy set to be the day within
310  // the year in the range from 1 to 366
311  *month = -1;
312  if ( ifleapyear )
313  {
314  while ( doy > MonthStartDOY_L[*month + 1] )
315  {
316  ( *month )++;
317  if ( *month == 11 )
318  break;
319  }
320  *day = doy - MonthStartDOY_L[*month];
321  }
322  else
323  {
324  while ( doy > MonthStartDOY[*month + 1] )
325  {
326  ( *month )++;
327  if ( *month == 11 )
328  break;
329  }
330  *day = doy - MonthStartDOY[*month];
331  }
332 }
333 
334 const char * getDayOfWeek( const MJDtime *MJD )
335 {
336  static const char *dow = { "Wed\0Thu\0Fri\0Sat\0Sun\0Mon\0Tue" };
337  int d = MJD->base_day % 7;
338  if ( d < 0 )
339  d += 7;
340  return &( dow[d * 4] );
341 }
342 
343 const char * getLongDayOfWeek( const MJDtime *MJD )
344 {
345  static const char *dow = { "Wednesday\0Thursday\0\0Friday\0\0\0\0Saturday\0\0Sunday\0\0\0\0Monday\0\0\0\0Tuesday" };
346  int d = MJD->base_day % 7;
347  if ( d < 0 )
348  d += 7;
349  return &( dow[d * 10] );
350 }
351 
352 const char * getMonth( int m )
353 {
354  static const char *months = { "Jan\0Feb\0Mar\0Apr\0May\0Jun\0Jul\0Aug\0Sep\0Oct\0Nov\0Dec" };
355  return &( months[( m ) * 4] );
356 }
357 
358 const char * getLongMonth( int m )
359 {
360  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" };
361  return &( months[( m ) * 10] );
362 }
363 
364 
365 size_t strfMJD( char * buf, size_t len, const char *format, const MJDtime *MJD, int forceJulian, int inleap )
366 {
367  // Format a text string according to the format string.
368  // Uses the same syntax as strftime() but does not use current locale.
369  // The null terminator is included in len for safety.
370 
371  // if inleap is true (non-zero) then renormalize so that (int) sec
372  // is 60 to mark results as a flag that a leap increment (recently
373  // always a second, but historically it was sometimes smaller than
374  // that) was in the process of being inserted for this particular
375  // epoch just prior to a positive discontinuity in TAI-UTC.
376 
377  int year, month, day, hour, min, ysign, second, d, y;
378  int y1, ifleapyear;
379  int i, secsSince1970;
380  int nplaces, fmtlen, slen;
381  int resolution;
382  double shiftPlaces;
383  char * ptr;
384  double sec, sec_fraction;
385  int w, doy, days_in_wk1;
386  const char *dayText;
387  const char *monthText;
388  char DateTime[80];
389  size_t posn = 0;
390  size_t last = len - 1;
391  MJDtime nMJD_value, *nMJD = &nMJD_value;
392  char dynamic_format[10];
393 
394  // Find required resolution
395  resolution = 0;
396  fmtlen = (int) strlen( format );
397  i = 0;
398  while ( i < fmtlen )
399  {
400  char next = format[i];
401  if ( next == '%' )
402  {
403  // find seconds format if used
404  i++;
405  next = format[i];
406  if ( isdigit( next ) != 0 )
407  {
408  nplaces = (int) strtol( &( format[i] ), NULL, 10 );
409  if ( nplaces > resolution )
410  resolution = nplaces;
411  }
412  else if ( next == '.' )
413  {
414  resolution = 9; // maximum resolution allowed
415  }
416  }
417  i++;
418  }
419 
420  // ensure rounding is done before breakdown
421  shiftPlaces = pow( 10, (double) resolution );
422  *nMJD = *MJD;
423  nMJD->time_sec += 0.5 / shiftPlaces;
424 
425  buf[last] = '\0';
426  buf[0] = '\0'; // force overwrite of old buffer since strnctat() used hereafter
427 
428  if ( inleap )
429  nMJD->time_sec -= 1.;
430 
431  breakDownMJD( &year, &month, &day, &hour, &min, &sec, nMJD, forceJulian );
432  if ( inleap )
433  sec += 1.;
434 
435  if ( year < 0 )
436  {
437  ysign = 1;
438  year = -year;
439  }
440  else
441  ysign = 0;
442 
443  //truncate seconds to resolution to stop formatting rounding up
444  sec = floor( sec * shiftPlaces ) / shiftPlaces;
445  second = (int) sec;
446 
447  // Read format string, character at a time
448  i = 0;
449  while ( i < fmtlen )
450  {
451  char next = format[i];
452  if ( next == '%' )
453  {
454  // format character or escape
455  i++;
456  next = format[i];
457  if ( next == '%' )
458  {
459  // escaped %, pass it on
460  buf[posn] = next;
461  posn++;
462  if ( posn >= last )
463  return posn;
464  }
465  else if ( next == 'a' )
466  {
467  // short day name
468  dayText = getDayOfWeek( nMJD );
469  strncat( &( buf[posn] ), dayText, last - posn );
470  posn = strlen( buf );
471  if ( posn >= last )
472  return posn;
473  }
474  else if ( next == 'A' )
475  {
476  // long day name
477  dayText = getLongDayOfWeek( nMJD );
478  strncat( &( buf[posn] ), dayText, last - posn );
479  posn = strlen( buf );
480  if ( posn >= last )
481  return posn;
482  }
483  else if ( next == 'b' || next == 'h' )
484  {
485  // short month name
486  monthText = getMonth( month );
487  strncat( &( buf[posn] ), monthText, last - posn );
488  posn = strlen( buf );
489  if ( posn >= last )
490  return posn;
491  }
492  else if ( next == 'B' )
493  {
494  // long month name
495  monthText = getLongMonth( month );
496  strncat( &( buf[posn] ), monthText, last - posn );
497  posn = strlen( buf );
498  if ( posn >= last )
499  return posn;
500  }
501  else if ( next == 'c' )
502  {
503  // Date and Time with day of week
504  dayText = getDayOfWeek( nMJD );
505  monthText = getMonth( month );
506  if ( ysign == 0 )
507  sprintf( DateTime, "%s %s %02d %02d:%02d:%02d %04d", dayText, monthText, day, hour, min, second, year );
508  else
509  sprintf( DateTime, "%s %s %02d %02d:%02d:%02d -%04d", dayText, monthText, day, hour, min, second, year );
510 
511  strncat( &( buf[posn] ), DateTime, last - posn );
512  posn = strlen( buf );
513  if ( posn >= last )
514  return posn;
515  }
516  else if ( next == 'C' )
517  {
518  // year / 100 so, e.g. 1989 is 20th century but comes out as 19
519  int century = year / 100;
520  if ( ysign == 0 )
521  sprintf( DateTime, "%02d", century );
522  else
523  sprintf( DateTime, "-%02d", century + 1 );
524 
525  strncat( &( buf[posn] ), DateTime, last - posn );
526  posn = strlen( buf );
527  if ( posn >= last )
528  return posn;
529  }
530  else if ( next == 'd' )
531  {
532  // day of month (01 - 31)
533  sprintf( DateTime, "%02d", day );
534 
535  strncat( &( buf[posn] ), DateTime, last - posn );
536  posn = strlen( buf );
537  if ( posn >= last )
538  return posn;
539  }
540  else if ( next == 'D' )
541  {
542  // month/day/year
543  y = year % 100;
544  if ( ysign == 0 )
545  sprintf( DateTime, "%02d/%02d/%02d", month + 1, day, y );
546  else
547  sprintf( DateTime, "%02d/%02d/-%02d", month + 1, day, y );
548 
549  strncat( &( buf[posn] ), DateTime, last - posn );
550  posn = strlen( buf );
551  if ( posn >= last )
552  return posn;
553  }
554  else if ( next == 'e' )
555  {
556  // day of month ( 1 - 31)
557  if ( day < 10 )
558  sprintf( DateTime, " %01d", day );
559  else
560  sprintf( DateTime, "%02d", day );
561 
562  strncat( &( buf[posn] ), DateTime, last - posn );
563  posn = strlen( buf );
564  if ( posn >= last )
565  return posn;
566  }
567  else if ( next == 'F' )
568  {
569  // year-month-day
570  if ( ysign == 0 )
571  sprintf( DateTime, "%04d-%02d-%02d", year, month + 1, day );
572  else
573  sprintf( DateTime, "-%04d-%02d-%02d", year, month + 1, day );
574 
575  strncat( &( buf[posn] ), DateTime, last - posn );
576  posn = strlen( buf );
577  if ( posn >= last )
578  return posn;
579  }
580  else if ( next == 'H' )
581  {
582  // hour, 24 hour clock (00 - 23)
583  sprintf( DateTime, "%02d", hour );
584 
585  strncat( &( buf[posn] ), DateTime, last - posn );
586  posn = strlen( buf );
587  if ( posn >= last )
588  return posn;
589  }
590  else if ( next == 'I' )
591  {
592  // hour, 12 hour clock (01 - 12)
593  if ( hour == 0 )
594  sprintf( DateTime, "%02d", hour + 12 );
595  else if ( hour > 12 )
596  sprintf( DateTime, "%02d", hour - 12 );
597  else
598  sprintf( DateTime, "%02d", hour );
599 
600  strncat( &( buf[posn] ), DateTime, last - posn );
601  posn = strlen( buf );
602  if ( posn >= last )
603  return posn;
604  }
605  else if ( next == 'j' )
606  {
607  // day of year
608  getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian );
609  sprintf( DateTime, "%03d", doy );
610 
611  strncat( &( buf[posn] ), DateTime, last - posn );
612  posn = strlen( buf );
613  if ( posn >= last )
614  return posn;
615  }
616  else if ( next == 'k' )
617  {
618  // hour, 24 hour clock ( 0 - 23)
619  if ( hour < 10 )
620  sprintf( DateTime, " %01d", hour );
621  else
622  sprintf( DateTime, "%02d", hour );
623 
624  strncat( &( buf[posn] ), DateTime, last - posn );
625  posn = strlen( buf );
626  if ( posn >= last )
627  return posn;
628  }
629  else if ( next == 'l' )
630  {
631  // hour, 12 hour clock ( 1 - 12)
632  if ( hour == 0 )
633  sprintf( DateTime, "%02d", hour + 12 );
634  else if ( hour < 10 )
635  sprintf( DateTime, " %01d", hour );
636  else if ( hour <= 12 )
637  sprintf( DateTime, "%02d", hour );
638  else if ( hour < 22 )
639  sprintf( DateTime, " %01d", hour - 12 );
640  else
641  sprintf( DateTime, "%02d", hour - 12 );
642 
643  strncat( &( buf[posn] ), DateTime, last - posn );
644  posn = strlen( buf );
645  if ( posn >= last )
646  return posn;
647  }
648  else if ( next == 'm' )
649  {
650  // month (01 - 12)
651  sprintf( DateTime, "%02d", month + 1 );
652 
653  strncat( &( buf[posn] ), DateTime, last - posn );
654  posn = strlen( buf );
655  if ( posn >= last )
656  return posn;
657  }
658  else if ( next == 'M' )
659  {
660  // minute (00 - 59)
661  sprintf( DateTime, "%02d", min );
662 
663  strncat( &( buf[posn] ), DateTime, last - posn );
664  posn = strlen( buf );
665  if ( posn >= last )
666  return posn;
667  }
668  else if ( next == 'n' )
669  {
670  // newline
671  buf[posn] = '\n';
672  posn++;
673  if ( posn >= last )
674  return posn;
675  }
676  else if ( next == 'p' )
677  {
678  // am/pm on12 hour clock
679  if ( hour < 0 )
680  sprintf( DateTime, "AM" );
681  else
682  sprintf( DateTime, "PM" );
683 
684  strncat( &( buf[posn] ), DateTime, last - posn );
685  posn = strlen( buf );
686  if ( posn >= last )
687  return posn;
688  }
689  else if ( next == 'r' )
690  {
691  // hour:min:sec AM, 12 hour clock (01 - 12):(00 - 59):(00 - 59) (AM - PM)
692  if ( hour == 0 )
693  sprintf( DateTime, "%02d:%02d:%02d AM", hour + 12, min, second );
694  else if ( hour > 12 )
695  sprintf( DateTime, "%02d:%02d:%02d PM", hour - 12, min, second );
696  else if ( hour == 12 )
697  sprintf( DateTime, "%02d:%02d:%02d PM", hour, min, second );
698  else
699  sprintf( DateTime, "%02d:%02d:%02d AM", hour, min, second );
700 
701  strncat( &( buf[posn] ), DateTime, last - posn );
702  posn = strlen( buf );
703  if ( posn >= last )
704  return posn;
705  }
706  else if ( next == 'R' )
707  {
708  // hour:min, 24 hour clock (00 - 23):(00 - 59)
709  sprintf( DateTime, "%02d:%02d", hour, min );
710 
711  strncat( &( buf[posn] ), DateTime, last - posn );
712  posn = strlen( buf );
713  if ( posn >= last )
714  return posn;
715  }
716  else if ( next == 'S' )
717  {
718  // second (00 - 59 with optional decimal point and numbers after the decimal point.)
719  if ( i + 2 < fmtlen && format[i + 1] == '%' && ( format[i + 2] == '.' || isdigit( format[i + 2] ) != 0 ) )
720  {
721  // nplaces is number of decimal places ( 0 < nplaces <= 9 )
722  if ( format[i + 2] == '.' )
723  // maximum number of places
724  nplaces = 9;
725  else
726  nplaces = (int) strtol( &( format[i + 2] ), NULL, 10 );
727  i += 2;
728  }
729  else
730  {
731  nplaces = 0;
732  }
733 
734  if ( nplaces == 0 )
735  {
736  sprintf( DateTime, "%02d", (int) ( sec + 0.5 ) );
737  }
738  else
739  {
740  sprintf( dynamic_format, "%%0%d.%df", nplaces + 3, nplaces );
741  sprintf( DateTime, dynamic_format, sec );
742  if ( format[i] == '.' )
743  {
744  slen = (int) strlen( DateTime ) - 1;
745  while ( DateTime[slen] == '0' && DateTime[slen - 1] != '.' )
746  {
747  DateTime[slen] = '\0'; // remove trailing zeros
748  slen--;
749  }
750  }
751  }
752 
753  strncat( &( buf[posn] ), DateTime, last - posn );
754  posn = strlen( buf );
755  if ( posn >= last )
756  return posn;
757  }
758  else if ( next == 's' )
759  {
760  // seconds since 01 Jan 1970 Gregorian
761  secsSince1970 = (int) ( nMJD->time_sec + ( nMJD->base_day - MJD_1970 ) * SecInDay );
762  sprintf( DateTime, "%d", secsSince1970 );
763 
764  strncat( &( buf[posn] ), DateTime, last - posn );
765  posn = strlen( buf );
766  if ( posn >= last )
767  return posn;
768  }
769  else if ( next == 't' )
770  {
771  // tab
772  buf[posn] = '\t';
773  posn++;
774  if ( posn >= last )
775  return posn;
776  }
777  else if ( next == 'T' )
778  {
779  // hour:min:sec, 24 hour clock (00 - 23):(00 - 59):(00 - 59)
780  sprintf( DateTime, "%02d:%02d:%02d", hour, min, second );
781 
782  strncat( &( buf[posn] ), DateTime, last - posn );
783  posn = strlen( buf );
784  if ( posn >= last )
785  return posn;
786  }
787  else if ( next == 'U' )
788  {
789  // week of year as a number, (00 - 53) start of week is Sunday
790  getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian );
791  days_in_wk1 = ( nMJD->base_day - doy - 4 ) % 7;
792 
793  w = ( doy + 6 - days_in_wk1 ) / 7;
794 
795  sprintf( DateTime, "%02d", w );
796 
797  strncat( &( buf[posn] ), DateTime, last - posn );
798  posn = strlen( buf );
799  if ( posn >= last )
800  return posn;
801  }
802  else if ( next == 'u' )
803  {
804  // weekday as a number, 0 = Monday
805  d = 1 + ( nMJD->base_day - 5 ) % 7;
806 
807  sprintf( DateTime, "%01d", d );
808 
809  strncat( &( buf[posn] ), DateTime, last - posn );
810  posn = strlen( buf );
811  if ( posn >= last )
812  return posn;
813  }
814  else if ( next == 'v' )
815  {
816  // day-MonthName-year day of month ( 1 - 31) - month (Jan ... Dec) - year (yyyy)
817 
818  monthText = getMonth( month );
819 
820  if ( ysign == 0 )
821  {
822  if ( day < 10 )
823  sprintf( DateTime, " %01d-%s-%04d", day, monthText, year );
824  else
825  sprintf( DateTime, "%02d-%s-%04d", day, monthText, year );
826  }
827  else
828  {
829  if ( day < 10 )
830  sprintf( DateTime, " %01d-%s-(-)%04d", day, monthText, year );
831  else
832  sprintf( DateTime, "%02d-%s-(-)%04d", day, monthText, year );
833  }
834 
835  strncat( &( buf[posn] ), DateTime, last - posn );
836  posn = strlen( buf );
837  if ( posn >= last )
838  return posn;
839  }
840  else if ( next == 'V' )
841  {
842  // week of year as a number, (01 - 53) start of week is Monday and first week has at least 3 days in year
843  getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian );
844  days_in_wk1 = ( nMJD->base_day - doy - 3 ) % 7;
845 
846  if ( days_in_wk1 <= 3 )
847  w = ( doy + 6 - days_in_wk1 ) / 7; // ensure first week has at least 3 days in this year
848  else
849  w = 1 + ( doy + 6 - days_in_wk1 ) / 7;
850 
851  if ( w == 0 )
852  w = 53;
853  sprintf( DateTime, "%02d", w );
854 
855  strncat( &( buf[posn] ), DateTime, last - posn );
856  posn = strlen( buf );
857  if ( posn >= last )
858  return posn;
859  }
860  else if ( next == 'w' )
861  {
862  // weekday as a number, 0 = Sunday
863  d = ( nMJD->base_day - 4 ) % 7;
864 
865  sprintf( DateTime, "%01d", d );
866 
867  strncat( &( buf[posn] ), DateTime, last - posn );
868  posn = strlen( buf );
869  if ( posn >= last )
870  return posn;
871  }
872  else if ( next == 'W' )
873  {
874  // week of year as a number, (00 - 53) start of week is Monday
875  getYAD( &y1, &ifleapyear, &doy, nMJD, forceJulian );
876  days_in_wk1 = ( nMJD->base_day - doy - 3 ) % 7;
877 
878  w = ( doy + 6 - days_in_wk1 ) / 7;
879 
880  sprintf( DateTime, "%02d", w );
881 
882  strncat( &( buf[posn] ), DateTime, last - posn );
883  posn = strlen( buf );
884  if ( posn >= last )
885  return posn;
886  }
887  else if ( next == 'x' )
888  {
889  // date string
890  dayText = getDayOfWeek( nMJD );
891  monthText = getMonth( month );
892  if ( ysign == 0 )
893  sprintf( DateTime, "%s %s %02d, %04d", dayText, monthText, day, year );
894  else
895  sprintf( DateTime, "%s %s %02d, -%04d", dayText, monthText, day, year );
896 
897  strncat( &( buf[posn] ), DateTime, last - posn );
898  posn = strlen( buf );
899  if ( posn >= last )
900  return posn;
901  }
902  else if ( next == 'X' )
903  {
904  // time string
905  sprintf( DateTime, "%02d:%02d:%02d", hour, min, second );
906 
907  strncat( &( buf[posn] ), DateTime, last - posn );
908  posn = strlen( buf );
909  if ( posn >= last )
910  return posn;
911  }
912  else if ( next == 'y' )
913  {
914  // 2 digit year
915  y = year % 100;
916 
917  if ( ysign == 0 )
918  sprintf( DateTime, "%02d", y );
919  else
920  sprintf( DateTime, "-%02d", y );
921 
922  strncat( &( buf[posn] ), DateTime, last - posn );
923  posn = strlen( buf );
924  if ( posn >= last )
925  return posn;
926  }
927  else if ( next == 'Y' )
928  {
929  // 4 digit year
930  if ( ysign == 0 )
931  sprintf( DateTime, "%04d", year );
932  else
933  sprintf( DateTime, "-%04d", year );
934 
935  strncat( &( buf[posn] ), DateTime, last - posn );
936  posn = strlen( buf );
937  if ( posn >= last )
938  return posn;
939  }
940  else if ( next == 'Z' )
941  {
942  // time zone and calendar, always UTC
943  if ( forceJulian )
944  strncat( &( buf[posn] ), "UTC Julian", last - posn );
945  else
946  strncat( &( buf[posn] ), "UTC Gregorian", last - posn );
947 
948  posn = strlen( buf );
949  if ( posn >= last )
950  return posn;
951  }
952  else if ( next == 'z' )
953  {
954  // time zone, always UTC
955  strncat( &( buf[posn] ), "+0000", last - posn );
956  posn = strlen( buf );
957  if ( posn >= last )
958  return posn;
959  }
960  else if ( next == '+' )
961  {
962  // date and time
963  dayText = getDayOfWeek( nMJD );
964  monthText = getMonth( month );
965  if ( ysign == 0 )
966  sprintf( DateTime, "%s %s %02d %02d:%02d:%02d UTC %04d", dayText, monthText, day, hour, min, second, year );
967  else
968  sprintf( DateTime, "%s %s %02d %02d:%02d:%02d UTC -%04d", dayText, monthText, day, hour, min, second, year );
969 
970  strncat( &( buf[posn] ), DateTime, last - posn );
971  posn = strlen( buf );
972  if ( posn >= last )
973  return posn;
974  }
975  else if ( next == '.' || isdigit( next ) != 0 )
976  {
977  // nplaces is number of decimal places ( 0 < nplaces <= 9 )
978  if ( next == '.' )
979  // maximum number of places
980  nplaces = 9;
981  else
982  nplaces = (int) strtol( &( format[i] ), NULL, 10 );
983 
984  // fractional part of seconds to maximum available accuracy
985  sec_fraction = sec - (int) sec;
986  sprintf( dynamic_format, "%%-%d.%df", nplaces + 2, nplaces );
987  //sprintf(DateTime, "%-11.9f", sec_fraction);
988  sprintf( DateTime, dynamic_format, sec_fraction );
989  while ( ( ptr = strrchr( &( DateTime[0] ), ' ' ) ) != NULL )
990  ptr[0] = '\0'; // remove trailing white space
991 
992  if ( next == '.' )
993  {
994  slen = (int) strlen( DateTime ) - 1;
995  while ( DateTime[slen] == '0' && DateTime[slen - 1] != '.' )
996  {
997  DateTime[slen] = '\0'; // remove trailing zeros
998  slen--;
999  }
1000  }
1001 
1002  ptr = strchr( DateTime, '.' );
1003  // remove everything in front of the decimal point, and
1004  // ignore case (%0) where no decimal point exists at all.
1005  if ( ptr != NULL )
1006  strncat( &( buf[posn] ), ptr, last - posn );
1007  posn = strlen( buf );
1008  if ( posn >= last )
1009  return posn;
1010  }
1011  }
1012  else
1013  {
1014  // regular multi-byte character, pass it on
1015  buf[posn] = next;
1016  posn++;
1017  if ( posn >= last )
1018  return posn;
1019  }
1020  buf[posn] = '\0';
1021  i++;
1022  }
1023  return posn;
1024 }
1025 
1026 int geMJDtime_TAI( const MJDtime *number1, const TAI_UTC *number2 )
1027 {
1028  // Returns true if number1 >= number2.
1029  // N.B. both number1 and number2 must be normalized.
1030  if ( number1->base_day > number2->base_day )
1031  {
1032  return 1;
1033  }
1034  else if ( number1->base_day < number2->base_day )
1035  {
1036  return 0;
1037  }
1038  else
1039  {
1040  return ( number1->time_sec >= number2->time_sec_tai );
1041  }
1042 }
1043 
1044 int geMJDtime_UTC( const MJDtime *number1, const TAI_UTC *number2 )
1045 {
1046  // Returns true if number1 >= number2.
1047  // N.B. both number1 and number2 must be normalized.
1048  if ( number1->base_day > number2->base_day )
1049  {
1050  return 1;
1051  }
1052  else if ( number1->base_day < number2->base_day )
1053  {
1054  return 0;
1055  }
1056  else
1057  {
1058  return ( number1->time_sec >= number2->time_sec_utc );
1059  }
1060 }
1061 
1062 double leap_second_TAI( const MJDtime *MJD_TAI, int *inleap, int *index )
1063 {
1064  // Logic assumes input MJD_TAI is in TAI
1065  // *inleap lets the calling routine know whether MJD_TAI corresponds
1066  // to an epoch when a positive leap increment is being inserted.
1067 
1068  MJDtime MJD_value, *MJD = &MJD_value;
1069  double leap;
1070  int debug = 0;
1071  // N.B. geMJDtime_TAI only works for normalized values.
1072  *MJD = *MJD_TAI;
1073  normalize_MJD( MJD );
1074  // Search for index such that TAI_UTC_lookup_table[*index] <= MJD(TAI) < TAI_UTC_lookup_table[*index+1]
1075  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 );
1076  if ( debug == 2 )
1077  fprintf( stderr, "*index = %d\n", *index );
1078  if ( *index == -1 )
1079  {
1080  // MJD is less than first table entry.
1081  // Debug: check that condition is met
1082  if ( debug && geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index + 1] ) )
1083  {
1084  fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index );
1085  exit( EXIT_FAILURE );
1086  }
1087  // There is (by assertion) no discontinuity at the start of the table.
1088  // Therefore, *inleap cannot be true.
1089  *inleap = 0;
1090  // Use initial offset for MJD values before first table entry.
1091  // Calculate this offset strictly from offset1. The slope term
1092  // doesn't enter because offset2 is the same as the UTC of the
1093  // first epoch of the table.
1094  return -TAI_UTC_lookup_table[*index + 1].offset1;
1095  }
1096  else if ( *index == number_of_entries_in_tai_utc_table - 1 )
1097  {
1098  // MJD is greater than or equal to last table entry.
1099  // Debug: check that condition is met
1100  if ( debug && !geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index] ) )
1101  {
1102  fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index );
1103  exit( EXIT_FAILURE );
1104  }
1105  // If beyond end of table, cannot be in middle of leap second insertion.
1106  *inleap = 0;
1107  // Use final offset for MJD values after last table entry.
1108  // The slope term doesn't enter because modern values of the slope
1109  // are zero.
1110  return -TAI_UTC_lookup_table[*index].offset1;
1111  }
1112  else if ( *index >= 0 && *index < number_of_entries_in_tai_utc_table )
1113  {
1114  // table[*index] <= MJD < table[*index+1].
1115  // Debug: check that condition is met
1116  if ( debug && !( geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index] ) && !geMJDtime_TAI( MJD, &TAI_UTC_lookup_table[*index + 1] ) ) )
1117  {
1118  fprintf( stderr, "MJD = {%d, %f}\n", MJD->base_day, MJD->time_sec );
1119  fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index );
1120  exit( EXIT_FAILURE );
1121  }
1122  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 );
1123  // Convert MJD(TAI) to normalized MJD(UTC).
1124  MJD->time_sec += leap;
1125  normalize_MJD( MJD );
1126  // If MJD(UT) is in the next interval of the corresponding
1127  // TAI_UTC_lookup_table, then we are right in the middle of a
1128  // leap interval (recently a second but for earlier epochs it could be
1129  // less) insertion. Note this logic even works when leap intervals
1130  // are taken away from UTC (i.e., leap is positive) since in that
1131  // case the UTC *index always corresponds to the TAI *index.
1132  *inleap = geMJDtime_UTC( MJD, &TAI_UTC_lookup_table[*index + 1] );
1133  return leap;
1134  }
1135  else
1136  {
1137  fprintf( stderr, "libqsastime (leap_second_TAI) logic ERROR: bad *index = %d\n", *index );
1138  exit( EXIT_FAILURE );
1139  }
1140 }
1141 
1142 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 )
1143 {
1144  // Configure the transformation between continuous time and broken-down time
1145  // that is used for ctimeqsas, btimeqsas, and strfqsas.
1146  int forceJulian, ret;
1147  MJDtime MJD_value, *MJD = &MJD_value;
1148 
1149  // Allocate memory for *qsasconfig if that hasn't been done by a
1150  // previous call.
1151  if ( *qsasconfig == NULL )
1152  {
1153  *qsasconfig = (QSASConfig *) malloc( (size_t) sizeof ( QSASConfig ) );
1154  if ( *qsasconfig == NULL )
1155  {
1156  fprintf( stderr, "configqsas: out of memory\n" );
1157  exit( EXIT_FAILURE );
1158  }
1159  }
1160 
1161  // Set bhunt_search index to a definite value less than -1 to insure no
1162  // initial hunt phase from some random index value is done.
1163  ( *qsasconfig )->index = -40;
1164 
1165  if ( scale != 0. )
1166  {
1167  if ( ifbtime_offset )
1168  {
1169  if ( ccontrol & 0x1 )
1170  forceJulian = 1;
1171  else
1172  forceJulian = 0;
1173  ret = setFromUT( year, month, day, hour, min, sec, MJD, forceJulian );
1174  if ( ret )
1175  {
1176  fprintf( stderr, "configqsas: some problem with broken-down arguments\n" );
1177  exit( EXIT_FAILURE );
1178  }
1179  offset1 = (double) MJD->base_day;
1180  offset2 = MJD->time_sec / (double) SecInDay;
1181  }
1182  ( *qsasconfig )->scale = scale;
1183  ( *qsasconfig )->offset1 = offset1;
1184  ( *qsasconfig )->offset2 = offset2;
1185  ( *qsasconfig )->ccontrol = ccontrol;
1186  }
1187  else
1188  {
1189  // if scale is 0., then use default values. Currently, that
1190  // default is continuous time (stored as a double) is seconds since
1191  // 1970-01-01 while broken-down time is Gregorian with no other
1192  // additional corrections.
1193  ( *qsasconfig )->scale = 1. / (double) SecInDay;
1194  ( *qsasconfig )->offset1 = (double) MJD_1970;
1195  ( *qsasconfig )->offset2 = 0.;
1196  ( *qsasconfig )->ccontrol = 0x0;
1197  }
1198 }
1199 
1200 void closeqsas( QSASConfig **qsasconfig )
1201 {
1202  // Close library if it has been opened.
1203  if ( *qsasconfig != NULL )
1204  {
1205  free( (void *) *qsasconfig );
1206  *qsasconfig = NULL;
1207  }
1208 }
1209 
1210 int ctimeqsas( int year, int month, int day, int hour, int min, double sec, double * ctime, QSASConfig *qsasconfig )
1211 {
1212  MJDtime MJD_value, *MJD = &MJD_value;
1213  int forceJulian, ret;
1214 
1215  if ( qsasconfig == NULL )
1216  {
1217  fprintf( stderr, "libqsastime (ctimeqsas) ERROR: configqsas must be called first.\n" );
1218  exit( EXIT_FAILURE );
1219  }
1220 
1221  if ( qsasconfig->ccontrol & 0x1 )
1222  forceJulian = 1;
1223  else
1224  forceJulian = 0;
1225 
1226  ret = setFromUT( year, month, day, hour, min, sec, MJD, forceJulian );
1227  if ( ret )
1228  return ret;
1229  *ctime = ( ( (double) ( MJD->base_day ) - qsasconfig->offset1 ) - qsasconfig->offset2 + MJD->time_sec / (double) SecInDay ) / qsasconfig->scale;
1230  return 0;
1231 }
1232 
1233 void btimeqsas( int *year, int *month, int *day, int *hour, int *min, double *sec, double ctime, QSASConfig *qsasconfig )
1234 {
1235  MJDtime MJD_value, *MJD = &MJD_value;
1236  int forceJulian;
1237  double integral_offset1, integral_offset2, integral_scaled_ctime;
1238  int inleap;
1239 
1240  if ( qsasconfig == NULL )
1241  {
1242  fprintf( stderr, "libqsastime (btimeqsas) ERROR: configqsas must be called first.\n" );
1243  exit( EXIT_FAILURE );
1244  }
1245 
1246  MJD->time_sec = SecInDay * ( modf( qsasconfig->offset1, &integral_offset1 ) + modf( qsasconfig->offset2, &integral_offset2 ) + modf( ctime * qsasconfig->scale, &integral_scaled_ctime ) );
1247  MJD->base_day = (int) ( integral_offset1 + integral_offset2 + integral_scaled_ctime );
1248 
1249  if ( qsasconfig->ccontrol & 0x1 )
1250  forceJulian = 1;
1251  else
1252  forceJulian = 0;
1253 
1254  if ( qsasconfig->ccontrol & 0x2 )
1255  MJD->time_sec += leap_second_TAI( MJD, &inleap, &( qsasconfig->index ) );
1256  else
1257  inleap = 0;
1258 
1259  // If in the middle of a positive leap increment insertion, normalize the
1260  // broken-down result so that *sec exceeds 60 to mark the insertion
1261  // (similar to the way February 29 marks a leap day).
1262 
1263  if ( inleap )
1264  MJD->time_sec -= 1.;
1265 
1266  breakDownMJD( year, month, day, hour, min, sec, MJD, forceJulian );
1267  if ( inleap )
1268  *sec += 1.;
1269 }
1270 
1271 size_t strfqsas( char * buf, size_t len, const char *format, double ctime, QSASConfig *qsasconfig )
1272 {
1273  MJDtime MJD_value, *MJD = &MJD_value;
1274  int forceJulian;
1275  double integral_offset1, integral_offset2, integral_scaled_ctime;
1276  int inleap;
1277 
1278  if ( qsasconfig == NULL )
1279  {
1280  fprintf( stderr, "libqsastime (strfqsas) ERROR: configqsas must be called first.\n" );
1281  exit( EXIT_FAILURE );
1282  }
1283  MJD->time_sec = SecInDay * ( modf( qsasconfig->offset1, &integral_offset1 ) + modf( qsasconfig->offset2, &integral_offset2 ) + modf( ctime * qsasconfig->scale, &integral_scaled_ctime ) );
1284  MJD->base_day = (int) ( integral_offset1 + integral_offset2 + integral_scaled_ctime );
1285 
1286  if ( qsasconfig->ccontrol & 0x1 )
1287  forceJulian = 1;
1288  else
1289  forceJulian = 0;
1290 
1291  if ( qsasconfig->ccontrol & 0x2 )
1292  MJD->time_sec += leap_second_TAI( MJD, &inleap, &( qsasconfig->index ) );
1293  else
1294  inleap = 0;
1295 
1296  return strfMJD( buf, len, format, MJD, forceJulian, inleap );
1297 }
1298 
1299 // bhunt_search. Search an ordered table with a binary hunt phase and a
1300 // binary search phase.
1301 //
1302 // On entry *low is used to help the hunt phase speed up the binary
1303 // search when consecutive calls to bhunt_search are made with
1304 // similar key values. On exit, *low is adjusted such that
1305 // base[*low] <= key < base[(*low+1)] with the special cases of
1306 //low set to -1 to indicate the key < base[0] and *low set to n-1
1307 // to indicate base[n-1] <= key. The function *ge must return true (1)
1308 // if its first argument (the search key) is greater than or equal to
1309 // its second argument (a table entry). Otherwise it returns false
1310 // (0). Items in the array base must be in ascending order.
1311 
1312 void bhunt_search( const void *key, const void *base, int n, size_t size, int *low, int ( *ge )( const void *keyval, const void *datum ) )
1313 {
1314  const void *indexbase;
1315  int mid, high, hunt_inc = 1;
1316  // If previous search found below range, then assure one hunt cycle
1317  // just in case new key is also below range.
1318  if ( *low == -1 )
1319  *low = 0;
1320  // Protect against invalid or undefined *low values where hunt
1321  // is waste of time.
1322  if ( *low < 0 || *low >= n )
1323  {
1324  *low = -1;
1325  high = n;
1326  }
1327  else
1328  {
1329  // binary hunt phase where we are assured 0 <= *low < n
1330  indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) ( *low ) ) );
1331  if ( ( *ge )( key, indexbase ) )
1332  {
1333  high = ( *low ) + hunt_inc;
1334  indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) high ) );
1335  // indexbase is valid if high < n.
1336  while ( ( high < n ) && ( ( *ge )( key, indexbase ) ) )
1337  {
1338  *low = high;
1339  hunt_inc += hunt_inc;
1340  high = high + hunt_inc;
1341  indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) high ) );
1342  }
1343  if ( high >= n )
1344  high = n;
1345  // At this point, low is valid and base[low] <= key
1346  // and either key < base[high] for valid high or high = n.
1347  }
1348  else
1349  {
1350  high = *low;
1351  *low = high - hunt_inc;
1352  indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) ( *low ) ) );
1353  // indexbase is valid if(*low) >= 0
1354  while ( ( ( *low ) >= 0 ) && !( ( *ge )( key, indexbase ) ) )
1355  {
1356  high = *low;
1357  hunt_inc += hunt_inc;
1358  *low = ( *low ) - hunt_inc;
1359  indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) ( *low ) ) );
1360  }
1361  if ( ( *low ) < 0 )
1362  *low = -1;
1363  // At this point high is valid and key < base[high]
1364  // and either base[low] <= key for valid low or low = -1.
1365  }
1366  }
1367  // binary search phase where we are assured base[low] <= key < base[high]
1368  // when both low and high are valid with obvious special cases signalled
1369  // by low = -1 or high = n.
1370  while ( high - *low > 1 )
1371  {
1372  mid = *low + ( high - *low ) / 2;
1373  indexbase = (const void *) ( ( (const char *) base ) + ( size * (size_t) mid ) );
1374  if ( ( *ge )( key, indexbase ) )
1375  *low = mid;
1376  else
1377  high = mid;
1378  }
1379 }