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
setFromUT(int year,int month,int day,int hour,int min,double sec,MJDtime * MJD,int forceJulian)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
getYAD(int * year,int * ifleapyear,int * doy,const MJDtime * MJD,int forceJulian)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
normalize_MJD(MJDtime * MJD)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
breakDownMJD(int * year,int * month,int * day,int * hour,int * min,double * sec,const MJDtime * MJD,int forceJulian)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
getDayOfWeek(const MJDtime * MJD)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
getLongDayOfWeek(const MJDtime * MJD)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
getMonth(int m)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
getLongMonth(int m)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
strfMJD(char * buf,size_t len,const char * format,const MJDtime * MJD,int forceJulian,int inleap)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
geMJDtime_TAI(const MJDtime * number1,const TAI_UTC * number2)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
geMJDtime_UTC(const MJDtime * number1,const TAI_UTC * number2)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
leap_second_TAI(const MJDtime * MJD_TAI,int * inleap,int * index)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
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)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
closeqsas(QSASConfig ** qsasconfig)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
ctimeqsas(int year,int month,int day,int hour,int min,double sec,double * ctime,QSASConfig * qsasconfig)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
btimeqsas(int * year,int * month,int * day,int * hour,int * min,double * sec,double ctime,QSASConfig * qsasconfig)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
strfqsas(char * buf,size_t len,const char * format,double ctime,QSASConfig * qsasconfig)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
bhunt_search(const void * key,const void * base,int n,size_t size,int * low,int (* ge)(const void * keyval,const void * datum))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 }
1380