From: <ai...@us...> - 2009-02-06 01:18:00
|
Revision: 9462 http://plplot.svn.sourceforge.net/plplot/?rev=9462&view=rev Author: airwin Date: 2009-02-06 01:17:56 +0000 (Fri, 06 Feb 2009) Log Message: ----------- Protect against bad value of month index input to setFromUT. Within setFromUT, do much improved protection against integer overflows. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-06 00:03:55 UTC (rev 9461) +++ trunk/lib/qsastime/qsastime.c 2009-02-06 01:17:56 UTC (rev 9462) @@ -55,15 +55,11 @@ /* Note C libraries use Gregorian only from 14 Sept 1752 onwards */ int leaps, lastyear; + double dbase_day, time_sec, dextraDays; + int extraDays; - /* Approximate precaution to avoid overflowing integer portion of - MJD. MJD epoch is 1858-11-17. The months, days, etc., portion - of this calculation are only included to guard against - extremely large values being used for some/all of them. */ - if(abs(365.25*(year-1858) + 12.*(month-10) + (day-17) + hour/24. - + min/1440. + sec/86400.) > 2.e9) + if(month < 1 || month > 12) return 1; - if(year <= 0) { /* count leap years on Julian Calendar */ @@ -72,9 +68,16 @@ leaps = year / 4 - 1 ; /* (note leaps is negative here and year 0 (1 BCE) was a leap year */ if(year%4 == 0) - MJD->base_day = year * 365 + leaps + MonthStartDOY_L[month-1] + day - 678943; + /* left to right associativity means the double constant + value of 365. propagates to make all calculations be + done in double precision without the potential of + integer overflow. The result should be a double which + stores the expected exact integer results of the + calculation with exact representation unless the + result is much larger than the integer overflow limit. */ + dbase_day = year * 365. + leaps + MonthStartDOY_L[month-1] + day - 678943; else - MJD->base_day = year * 365 + leaps + MonthStartDOY[month-1] + day - 678943; + dbase_day = year * 365. + leaps + MonthStartDOY[month-1] + day - 678943; } else if(year < 1582 || (year == 1582 && month < 10) || (year == 1582 && month == 10 && day < 15) || forceJulian) @@ -84,35 +87,50 @@ leaps = (year -1 ) / 4; if(year%4 == 0) - MJD->base_day = year * 365 + leaps + MonthStartDOY_L[month-1] + day - 678943; + dbase_day = year * 365. + leaps + MonthStartDOY_L[month-1] + day - 678943; else - MJD->base_day = year * 365 + leaps + MonthStartDOY[month-1] + day - 678943; + dbase_day = year * 365. + leaps + MonthStartDOY[month-1] + day - 678943; } else { /* count leap years Gregorian Calendar - modern dates */ /* Algorithm below for 17 Nov 1858 (0 MJD) gives */ - /* leaps = 450 and hence base_day of 678941, so subtract it to give MJD day */ + /* leaps = 450 and hence dbase_day of 678941, so subtract it to give MJD day */ lastyear = year - 1; leaps = lastyear / 4 - lastyear / 100 + lastyear / 400; if( (year%4 == 0 && year%100 != 0) || (year%4 == 0 && year%400 == 0) ) - MJD->base_day = year * 365 + leaps + MonthStartDOY_L[month-1] + day - 678941; + dbase_day = year * 365. + leaps + MonthStartDOY_L[month-1] + day - 678941; else - MJD->base_day = year * 365 + leaps + MonthStartDOY[month-1] + day - 678941; + dbase_day = year * 365. + leaps + MonthStartDOY[month-1] + day - 678941; } - MJD->time_sec = sec + ( (double) min + (double) hour * 60. ) * 60.; + time_sec = sec + ( (double) min + (double) hour * 60. ) * 60.; - if(MJD->time_sec >= SecInDay) + if(time_sec >= SecInDay) { - int extraDays = (int) (MJD->time_sec / SecInDay); - MJD->base_day += extraDays; - MJD->time_sec -= extraDays * SecInDay; + dextraDays = (time_sec / SecInDay); + /* precaution against overflowing extraDays. */ + if(abs(dextraDays) > 2.e9) { + return 2; + } + extraDays = (int) (dextraDays); + dbase_day += extraDays; + time_sec -= extraDays * SecInDay; } - - return 0; + /* precaution against overflowing MJD->base_day. */ + if(abs(dbase_day) > 2.e9){ + return 3; + } else { + /* The exact integer result should be represented exactly in the + double, dbase_day, and its absolute value should be less than + the integer overflow limit. So the cast to int should be + exact. */ + MJD->base_day = (int) dbase_day; + MJD->time_sec = time_sec; + return 0; + } } const char * getDayOfWeek( const MJDtime *MJD) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-07 04:39:05
|
Revision: 9468 http://plplot.svn.sourceforge.net/plplot/?rev=9468&view=rev Author: airwin Date: 2009-02-07 04:39:01 +0000 (Sat, 07 Feb 2009) Log Message: ----------- Implement forceJulian == -1 logic which is interpreted as a force Gregorian mode to follow what is done for the Linux C library where the Gregorian calendar is assumed for all time for 64-bit systems. For 32-bit systems, time_t starts overflowing for dates before 1902, so worries about the Julian calendar (the initial transition from Julian to Gregorian occurred in 1582) are irrelevant for those systems. I assume that is why the Linux C library date interpretation is still Gregorian for all the much earlier dates that become accessible with 64-bit systems. Personally, I think sticking with Gregorian for all dates represented on a computer is an excellent idea from the KISS perspective. Long may it last!) Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-07 04:27:09 UTC (rev 9467) +++ trunk/lib/qsastime/qsastime.c 2009-02-07 04:39:01 UTC (rev 9468) @@ -26,7 +26,7 @@ /* MJD measures from the start of 17 Nov 1858 */ -/* These utilities use the Gregorian calendar after 4 Oct 1582 (Julian) i.e. from 15 Oct 1582 Gregoriam +/* These utilities use the Gregorian calendar after 4 Oct 1582 (Julian) i.e. from 15 Oct 1582 Gregorian Note C libraries use Gregorian only from 14 Sept 1752 More detailed discussion can be found at http://aa.usno.navy.mil/data/docs/JulianDate.php These routines have been compared with the results of the US Naval Observatory online converter. @@ -60,7 +60,9 @@ if(month < 0 || month > 11) return 1; - if(year <= 0) + if(forceJulian < -1 || forceJulian > 1) + return 2; + if(year <= 0 && forceJulian != -1) { /* count leap years on Julian Calendar */ /* MJD for Jan 1 0000 (correctly Jan 01, BCE 1) is - 678943, count from there */ @@ -80,7 +82,7 @@ dbase_day = year * 365. + leaps + MonthStartDOY[month] + day - 678943; } - else if(year < 1582 || (year == 1582 && month < 9) || (year == 1582 && month == 9 && day < 15) || forceJulian) + else if((forceJulian == 0 && (year < 1582 || (year == 1582 && month < 9) || (year == 1582 && month == 9 && day < 15))) || forceJulian == 1) { /* count leap years on Julian Calendar */ /* MJD for Jan 1 0000 (correctly Jan 01, BCE 1) is - 678943, count from there */ @@ -113,7 +115,7 @@ dextraDays = (time_sec / SecInDay); /* precaution against overflowing extraDays. */ if(abs(dextraDays) > 2.e9) { - return 2; + return 3; } extraDays = (int) (dextraDays); dbase_day += extraDays; @@ -121,7 +123,7 @@ } /* precaution against overflowing MJD->base_day. */ if(abs(dbase_day) > 2.e9){ - return 3; + return 4; } else { /* The exact integer result should be represented exactly in the double, dbase_day, and its absolute value should be less than This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-08 00:13:56
|
Revision: 9471 http://plplot.svn.sourceforge.net/plplot/?rev=9471&view=rev Author: airwin Date: 2009-02-08 00:13:51 +0000 (Sun, 08 Feb 2009) Log Message: ----------- Simplify and correct leap year code. Also improve documentation of that code. Result is correct near JD epoch (forceJulian=1), near year epoch (forceJulian=1 and forceJulian=-1), and near MJD epoch (forceJulian=-1). Conversion from broken-down to continuous time with setFromUT also agrees exactly with timegm results from the Linux C library for extremely wide range of proleptic Gregorian calendar, (-5000000 to + 5000000 years) with coarse sampling of 1000000 years. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-07 17:26:32 UTC (rev 9470) +++ trunk/lib/qsastime/qsastime.c 2009-02-08 00:13:51 UTC (rev 9471) @@ -54,57 +54,65 @@ /* default is to use Gregorian after 4 Oct 1582 (Julian) i.e. from 15 Oct 1582 Gregorian */ /* Note C libraries use Gregorian only from 14 Sept 1752 onwards */ - int leaps, lastyear; - double dbase_day, time_sec, dextraDays; + int leaps, year4, year100, year400; + double dbase_day, non_leaps = 365.; + //int dbase_day, non_leaps = 365; + double time_sec, dextraDays; int extraDays; if(month < 0 || month > 11) return 1; if(forceJulian < -1 || forceJulian > 1) return 2; - if(year <= 0 && forceJulian != -1) + /* As year increases, year4/4 increments by 1 at + year = -7, -3, 1, 5, 9, etc. */ + /* As year increases, year100/100 increments by 1 at + year = -299, -199, -99, 1, 101, 201, 301, etc. */ + /* As year increases, year400/400 increments by 1 at + year = -1199, -799, -399, 1, 401, 801, 1201, etc. */ + if(year <=0) { + year4 = year - 4; + year100 = year - 100; + year400 = year - 400; + } else { + year4 = year - 1; + year100 = year - 1; + year400 = year - 1; + } + if((forceJulian == 0 && (year < 1582 || (year == 1582 && month < 9) || (year == 1582 && month == 9 && day < 15))) || forceJulian == 1) { - /* count leap years on Julian Calendar */ + /* count leap years on proleptic Julian Calendar */ /* MJD for Jan 1 0000 (correctly Jan 01, BCE 1) is - 678943, count from there */ - /* negative CE (AD) years convert to BCE (BC) as BCE = 1 - CE, e.g. 2 BCE = -1 CE */ - - leaps = year / 4 - 1 ; /* (note leaps is negative here and year 0 (1 BCE) was a leap year */ - if(year%4 == 0) - /* left to right associativity means the double constant - value of 365. propagates to make all calculations be - done in double precision without the potential of - integer overflow. The result should be a double which - stores the expected exact integer results of the - calculation with exact representation unless the - result is much larger than the integer overflow limit. */ - dbase_day = year * 365. + leaps + MonthStartDOY_L[month] + day - 678943; - else - dbase_day = year * 365. + leaps + MonthStartDOY[month] + day - 678943; - - } - else if((forceJulian == 0 && (year < 1582 || (year == 1582 && month < 9) || (year == 1582 && month == 9 && day < 15))) || forceJulian == 1) - { - /* count leap years on Julian Calendar */ - /* MJD for Jan 1 0000 (correctly Jan 01, BCE 1) is - 678943, count from there */ - leaps = (year -1 ) / 4; + leaps = year4 / 4; if(year%4 == 0) - dbase_day = year * 365. + leaps + MonthStartDOY_L[month] + day - 678943; + dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day - 678943; else - dbase_day = year * 365. + leaps + MonthStartDOY[month] + day - 678943; + dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day - 678943; } else { - /* count leap years Gregorian Calendar - modern dates */ - /* Algorithm below for 17 Nov 1858 (0 MJD) gives */ - /* leaps = 450 and hence dbase_day of 678941, so subtract it to give MJD day */ + /* count leap years for proleptic Gregorian Calendar. */ + /* Algorithm below for 1858-11-17 (0 MJD) gives */ + /* leaps = 450 and hence dbase_day of 678941, so subtract that value */ + /* to give MJD day. Note this zero point is two days different from */ + /* Julian zero point above, because proleptic Gregorian is behind */ + /* Julian by two days at year = 1, see */ + /* http://en.wikipedia.org/wiki/Proleptic_Gregorian_calendar, */ - lastyear = year - 1; - leaps = lastyear / 4 - lastyear / 100 + lastyear / 400; + leaps = year4 / 4 - year100 / 100 + year400 / 400; + + /* left to right associativity means the double value of + non_leaps propagate to make all calculations be + done in double precision without the potential of + integer overflow. The result should be a double which + stores the expected exact integer results of the + calculation with exact representation unless the + result is much larger than the integer overflow limit. */ if( (year%4 == 0 && year%100 != 0) || (year%4 == 0 && year%400 == 0) ) - dbase_day = year * 365. + leaps + MonthStartDOY_L[month] + day - 678941; + dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day - 678941; else - dbase_day = year * 365. + leaps + MonthStartDOY[month] + day - 678941; + dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day - 678941; } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-10 00:27:14
|
Revision: 9485 http://plplot.svn.sourceforge.net/plplot/?rev=9485&view=rev Author: airwin Date: 2009-02-10 00:27:11 +0000 (Tue, 10 Feb 2009) Log Message: ----------- Indentation + tweak two comments. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-09 22:56:38 UTC (rev 9484) +++ trunk/lib/qsastime/qsastime.c 2009-02-10 00:27:11 UTC (rev 9485) @@ -193,36 +193,36 @@ j = MJD->base_day + extra_days; - if( j < -678943) { + if( j < -678943) { - /* Negative CE dates */ + /* Negative CE dates */ - j += 678943; + j += 678943; - /* negative years */ - year = (int) ((float)(j-365) / 365.25); - doy = j +1 - year * 365.25; + /* negative years */ + year = (int) ((float)(j-365) / 365.25); + doy = j +1 - year * 365.25; - } - else if( j < -678577) { + } + else if( j < -678577) { - /* CE = 0, BCE = 1 dates */ + /* CE = 0, BCE = 1 dates */ - j += 678943; + j += 678943; - /* negative years */ - year = 0; - doy = j +1; + /* negative years */ + year = 0; + doy = j +1; - } + } else if( j < -100840 || forceJulian == 1) { /* Julian Dates */ j += 678943; year = (int) ((float)j / 365.25); - lastyear = year - 1; - doy = j - year * 365 - lastyear / 4; + lastyear = year - 1; + doy = j - year * 365 - lastyear / 4; } @@ -267,17 +267,17 @@ /* BCE dates */ j += 678943; - if( j>=0) - { - *year = 0; - j++; - } - else - { - /* negative years */ - *year = (int) ((float)(j-365) / 365.25); - j = j +1 - *year * 365.25; - } + if( j>=0) + { + *year = 0; + j++; + } + else + { + /* negative years */ + *year = (int) ((float)(j-365) / 365.25); + j = j +1 - *year * 365.25; + } /* j is now always positive */ *month = -1; @@ -303,13 +303,13 @@ } else if( j < -100840 || forceJulian == 1) { - /* Julian Dates */ + /* Julian proleptic dates */ j += 678943; *year = (int) ((float)j / 365.25); - lastyear = *year - 1; - j = j - *year * 365 - lastyear / 4; + lastyear = *year - 1; + j = j - *year * 365 - lastyear / 4; *month = -1; if(*year%4 != 0) @@ -333,7 +333,7 @@ } else { - /* Gregorian Dates */ + /* Gregorian proleptic Dates */ j += 678941; *year = (int) ((float)j / 365.2425); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-10 04:05:45
|
Revision: 9486 http://plplot.svn.sourceforge.net/plplot/?rev=9486&view=rev Author: airwin Date: 2009-02-10 04:05:40 +0000 (Tue, 10 Feb 2009) Log Message: ----------- Use #defines for MJD_0000J, MJD_0000G, MJD_0001J, and MJD_0001G to replace numerical constants scattered through the code. Some curly brackets (control flow) style changes. Improve commentary. Initial attempt at forceJulian = -1 (Gregorian proleptic calendar) logic for getDOY (untested so far) and breakDownMJD (the first test shows some work still to do). Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-10 00:27:11 UTC (rev 9485) +++ trunk/lib/qsastime/qsastime.c 2009-02-10 04:05:40 UTC (rev 9486) @@ -39,6 +39,18 @@ */ #include <ctype.h> #include "qsastime.h" +/* MJD for 0000-01-01 (correctly Jan 01, BCE 1) */ +/* Julian proleptic calendar value. */ +#define MJD_0000J -678943 +/* Gregorian proleptic calendar value. (At MJD_0000J the Gregorian proleptic + calendar reads two days behind the Julian proleptic calendar, i.e. - 2 days, + see http://en.wikipedia.org/wiki/Proleptic_Gregorian_calendar, + so MJD_0000G = MJD_0000J+2) */ +#define MJD_0000G -678941 +/* MJD for 0001-01-01 which is 366 days later than previous definitions because + the year 0 is a leap year in both calendars.*/ +#define MJD_0001J -678577 +#define MJD_0001G -678575 static double SecInDay = 86400; /* we ignore leap seconds */ static int MJD_1970 = 40587; /* MJD for Jan 01, 1970 00:00:00 */ @@ -60,10 +72,14 @@ double time_sec, dextraDays; int extraDays; - if(month < 0 || month > 11) - return 1; - if(forceJulian < -1 || forceJulian > 1) - return 2; + if(month < 0 || month > 11) { + fprintf(stderr, "setfromUT: invalid month value\n"); + exit(EXIT_FAILURE); + } + if(forceJulian < -1 || forceJulian > 1) { + fprintf(stderr, "setfromUT: invalid forceJulian value\n"); + exit(EXIT_FAILURE); + } /* As year increases, year4/4 increments by 1 at year = -7, -3, 1, 5, 9, etc. */ /* As year increases, year100/100 increments by 1 at @@ -82,25 +98,20 @@ if((forceJulian == 0 && (year < 1582 || (year == 1582 && month < 9) || (year == 1582 && month == 9 && day < 15))) || forceJulian == 1) { - /* count leap years on proleptic Julian Calendar */ - /* MJD for Jan 1 0000 (correctly Jan 01, BCE 1) is - 678943, count from there */ - + /* count leap years on proleptic Julian Calendar starting from MJD_0000J */ leaps = year4 / 4; if(year%4 == 0) - dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day - 678943; + dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000J; else - dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day - 678943; + dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000J; } else { /* count leap years for proleptic Gregorian Calendar. */ - /* Algorithm below for 1858-11-17 (0 MJD) gives */ - /* leaps = 450 and hence dbase_day of 678941, so subtract that value */ - /* to give MJD day. Note this zero point is two days different from */ - /* Julian zero point above, because proleptic Gregorian is behind */ - /* Julian by two days at year = 1, see */ - /* http://en.wikipedia.org/wiki/Proleptic_Gregorian_calendar, */ - + /* Algorithm below for 1858-11-17 (0 MJD) gives + leaps = 450 and hence dbase_day of 678941, so subtract that value + or add MJD_0000G (which is two days different from MJD_0000J, see + above). */ leaps = year4 / 4 - year100 / 100 + year400 / 400; /* left to right associativity means the double value of @@ -111,9 +122,9 @@ calculation with exact representation unless the result is much larger than the integer overflow limit. */ if( (year%4 == 0 && year%100 != 0) || (year%4 == 0 && year%400 == 0) ) - dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day - 678941; + dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000G; else - dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day - 678941; + dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000G; } @@ -180,6 +191,10 @@ int extra_days,j,lastyear; + if(forceJulian < -1 || forceJulian > 1) { + fprintf(stderr, "getDOY: invalid forceJulian value\n"); + exit(EXIT_FAILURE); + } if(MJD->time_sec >= 0) { extra_days = (int) (MJD->time_sec / SecInDay); @@ -193,50 +208,43 @@ j = MJD->base_day + extra_days; - if( j < -678943) { + if(forceJulian != -1 && j < MJD_0000J) { - /* Negative CE dates */ + /* Change epoch so j is measured in days after 0000-01-01 on the + Julian proleptic calendar. */ - j += 678943; + j -= MJD_0000J; - /* negative years */ - year = (int) ((float)(j-365) / 365.25); + /* j must be strictly negative from above logic. Therefore, year must be + strictly negative as well. */ + year = (int) ((double)(j-365) / 365.25); doy = j +1 - year * 365.25; - } - else if( j < -678577) { + } else if(forceJulian != -1 && j < MJD_0001J ) { + /* in year 0 which was a leap year. */ - /* CE = 0, BCE = 1 dates */ - - j += 678943; + j -= MJD_0000J; - /* negative years */ year = 0; doy = j +1; - } - else if( j < -100840 || forceJulian == 1) - { - /* Julian Dates */ - j += 678943; - - year = (int) ((float)j / 365.25); - lastyear = year - 1; - doy = j - year * 365 - lastyear / 4; + } else if( forceJulian != -1 && (j < -100840 || forceJulian == 1)) { + j -= MJD_0000J; + year = (int) ((double)j / 365.25); + lastyear = year - 1; + doy = j - year * 365 - lastyear / 4; + } else { + /* forceJulian == -1 || (j >= -100840 && forceJulian == 0) */ + /* Change epoch so j is measured in days after 0000-01-01 on the + Gregorian proleptic calendar. */ + j -= MJD_0000G; - - } - else - { - /* Gregorian Dates */ - j += 678941; - - year = (int) ((float)j / 365.2425); - lastyear = year - 1; - doy = j - year * 365 - lastyear / 4 + lastyear / 100 - lastyear / 400; - - } + year = (int) ((float)j / 365.2425); + lastyear = year - 1; + doy = j - year * 365 - lastyear / 4 + lastyear / 100 - lastyear / 400; + } + return doy; } @@ -244,123 +252,101 @@ { /* Convert MJD struct into date/time elements */ /* Note year 0 CE (AD) [1 BCE (BC)] is a leap year */ - /* There are 678943 days from year 0 to MJD(0) */ int extra_days,j,lastyear; double seconds; - if(MJD->time_sec >= 0) - { + if(forceJulian < -1 || forceJulian > 1) { + fprintf(stderr, "breakDownMJD: invalid forceJulian value\n"); + exit(EXIT_FAILURE); + } + if(MJD->time_sec >= 0) { extra_days = (int) (MJD->time_sec / SecInDay); - } - else - { - /* allow for negative seconds push into previous day even if less than 1 day */ - extra_days = (int) (MJD->time_sec / SecInDay) - 1 ; - } - - + } else { + /* allow for negative seconds push into previous day even if less than 1 day */ + extra_days = (int) (MJD->time_sec / SecInDay) - 1 ; + } + j = MJD->base_day + extra_days; - if( j <= -678577) { + if( forceJulian != -1 && j <= MJD_0001J) { /* BCE dates */ - j += 678943; - if( j>=0) - { + j -= MJD_0000J; + if( j>=0) { *year = 0; j++; - } - else - { - /* negative years */ - *year = (int) ((float)(j-365) / 365.25); - j = j +1 - *year * 365.25; - } + } else { + /* negative years */ + *year = (int) ((float)(j-365) / 365.25); + j = j +1 - *year * 365.25; + } /* j is now always positive */ *month = -1; - if(*year%4 != 0) - { - while(j > MonthStartDOY[*month +1]) - { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY[*month]; + if(*year%4 != 0) { + while(j > MonthStartDOY[*month +1]) { + (*month)++; + if(*month == 11) break; } - else - { - while(j > MonthStartDOY_L[*month +1]) - { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY_L[*month]; + *day = j - MonthStartDOY[*month]; + } else { + while(j > MonthStartDOY_L[*month +1]) { + (*month)++; + if(*month == 11) break; } - } - else if( j < -100840 || forceJulian == 1) - { - /* Julian proleptic dates */ - j += 678943; - - *year = (int) ((float)j / 365.25); - - lastyear = *year - 1; - j = j - *year * 365 - lastyear / 4; - - *month = -1; - if(*year%4 != 0) - { - while(j > MonthStartDOY[*month +1]) - { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY[*month]; - } - else - { - while(j > MonthStartDOY_L[*month + 1]) - { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY_L[*month]; - } + *day = j - MonthStartDOY_L[*month]; } - else - { - /* Gregorian proleptic Dates */ - j += 678941; - - *year = (int) ((float)j / 365.2425); - lastyear = *year - 1; - j = j - *year * 365 - lastyear / 4 + lastyear / 100 - lastyear / 400; - - *month = -1; - if((*year%4 == 0 && *year%100 != 0) || (*year%4 == 0 && *year%400 == 0) ) - { - while(j > MonthStartDOY_L[*month +1]) - { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY_L[*month]; - } - else - { - while(j > MonthStartDOY[*month +1]) - { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY[*month]; - } - + } else if( forceJulian !=-1 && (j < -100840 || forceJulian == 1)) { + /* Julian proleptic dates */ + j -= MJD_0000J; + + *year = (int) ((float)j / 365.25); + + lastyear = *year - 1; + j = j - *year * 365 - lastyear / 4; + + *month = -1; + if(*year%4 != 0) { + while(j > MonthStartDOY[*month +1]) { + (*month)++; + if(*month == 11) break; + } + *day = j - MonthStartDOY[*month]; + } else { + while(j > MonthStartDOY_L[*month + 1]) { + (*month)++; + if(*month == 11) break; + } + *day = j - MonthStartDOY_L[*month]; } + } else { + /* forceJulian == -1 || (j >= -100840 && forceJulian == 0) */ + /* Gregorian proleptic Dates */ + j -= MJD_0000G; + + *year = (int) ((float)j / 365.2425); + lastyear = *year - 1; + j = j - *year * 365 - lastyear / 4 + lastyear / 100 - lastyear / 400; + + *month = -1; + if((*year%4 == 0 && *year%100 != 0) || (*year%4 == 0 && *year%400 == 0) ) { + while(j > MonthStartDOY_L[*month +1]) { + (*month)++; + if(*month == 11) break; + } + *day = j - MonthStartDOY_L[*month]; + } else { + while(j > MonthStartDOY[*month +1]) { + (*month)++; + if(*month == 11) break; + } + *day = j - MonthStartDOY[*month]; + } + + } /* Time part */ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-10 07:16:46
|
Revision: 9489 http://plplot.svn.sourceforge.net/plplot/?rev=9489&view=rev Author: airwin Date: 2009-02-10 07:16:44 +0000 (Tue, 10 Feb 2009) Log Message: ----------- Fix breakDownMJD forceJulian=-1 (forced Gregorian proleptic calendar) issues. Results pass test 1 of qsastime_testlib (which tests setFromUT, breakDownMJD, and strfMJD for both forceJulian=-1 and forceJulian=1 near the epoch of Julian date). Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-10 07:12:08 UTC (rev 9488) +++ trunk/lib/qsastime/qsastime.c 2009-02-10 07:16:44 UTC (rev 9489) @@ -253,7 +253,7 @@ /* Convert MJD struct into date/time elements */ /* Note year 0 CE (AD) [1 BCE (BC)] is a leap year */ - int extra_days,j,lastyear; + int extra_days,j, year4, year100, year400; double seconds; if(forceJulian < -1 || forceJulian > 1) { @@ -300,13 +300,23 @@ *day = j - MonthStartDOY_L[*month]; } } else if( forceJulian !=-1 && (j < -100840 || forceJulian == 1)) { - /* Julian proleptic dates */ + /* Shift j epoch to 0000-01-01 for the Julian proleptic calendar.*/ j -= MJD_0000J; - *year = (int) ((float)j / 365.25); + /* 365.25 is exact period of Julian year so year will be correct if + day offset is set exactly right so that years -4, 0, 4 are leap years, + i.e. years -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 start with j = + -1826 -1461, -1095, -730, -365, 0, 366, 731, 1096, 1461, 1827 */ + if(j >= 0) { + *year = (int) ((float)j / 365.25); + year4 = year-1; + } else { + *year = (int) ((float)j/ 365.25) - 1; + year4 = year-4; + } - lastyear = *year - 1; - j = j - *year * 365 - lastyear / 4; + /* shift j epoch to day of year */ + j = j - *year * 365 - year4 / 4; *month = -1; if(*year%4 != 0) { @@ -324,12 +334,24 @@ } } else { /* forceJulian == -1 || (j >= -100840 && forceJulian == 0) */ - /* Gregorian proleptic Dates */ + /* Shift j epoch to 0000-01-01 for the Gregorian proleptic calendar.*/ j -= MJD_0000G; + /* 365.245 is exact period of Gregorian year so year will be correct + if exactly correct j offset is chosen (see above). */ - *year = (int) ((float)j / 365.2425); - lastyear = *year - 1; - j = j - *year * 365 - lastyear / 4 + lastyear / 100 - lastyear / 400; + if(j >=0) { + *year = (int) ((float)j / 365.2425); + year4 = *year - 1; + year100 = *year - 1; + year400 = *year - 1; + } else { + *year = (int) ((float)j / 365.2425) - 1; + year4 = *year - 4; + year100 = *year - 100; + year400 = *year - 400; + } + + j = j - *year * 365 - year4 / 4 + year100 / 100 - year400 / 400; *month = -1; if((*year%4 == 0 && *year%100 != 0) || (*year%4 == 0 && *year%400 == 0) ) { This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-10 18:47:47
|
Revision: 9498 http://plplot.svn.sourceforge.net/plplot/?rev=9498&view=rev Author: airwin Date: 2009-02-10 18:47:43 +0000 (Tue, 10 Feb 2009) Log Message: ----------- Put functions in more logical order. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-10 12:50:17 UTC (rev 9497) +++ trunk/lib/qsastime/qsastime.c 2009-02-10 18:47:43 UTC (rev 9498) @@ -155,35 +155,6 @@ } } -const char * getDayOfWeek( const MJDtime *MJD) -{ - static char *dow = {"Wed\0Thu\0Fri\0Sat\0Sun\0Mon\0Tue"}; - int d = MJD->base_day % 7; - if(d < 0) d += 7; - return &(dow[d*4]); -} - -const char * getLongDayOfWeek( const MJDtime *MJD) -{ - static char *dow = {"Wednesday\0Thursday\0\0Friday\0\0\0\0Saturday\0\0Sunday\0\0\0\0Monday\0\0\0\0Tuesday"}; - int d = MJD->base_day % 7; - if(d < 0) d += 7; - return &(dow[d*10]); -} - -const char * getMonth( int m) -{ - static char *months = {"Jan\0Feb\0Mar\0Apr\0May\0Jun\0Jul\0Aug\0Sep\0Oct\0Nov\0Dec"}; - return &(months[(m)*4]); -} - -const char * getLongMonth( int m) -{ - static char *months = {"January\0\0\0February\0\0March\0\0\0\0\0April\0\0\0\0\0May\0\0\0\0\0\0\0June\0\0\0\0\0\0July\0\0\0\0\0\0August\0\0\0\0September\0October\0\0\0November\0\0December"}; - return &(months[(m)*10]); -} - - int getDOY(const MJDtime *MJD, int forceJulian) { /* Get from Day Of Year */ @@ -242,7 +213,6 @@ year = (int) ((float)j / 365.2425); lastyear = year - 1; doy = j - year * 365 - lastyear / 4 + lastyear / 100 - lastyear / 400; - } return doy; @@ -379,6 +349,35 @@ *sec = seconds - (double) *min * 60.; } +const char * getDayOfWeek( const MJDtime *MJD) +{ + static char *dow = {"Wed\0Thu\0Fri\0Sat\0Sun\0Mon\0Tue"}; + int d = MJD->base_day % 7; + if(d < 0) d += 7; + return &(dow[d*4]); +} + +const char * getLongDayOfWeek( const MJDtime *MJD) +{ + static char *dow = {"Wednesday\0Thursday\0\0Friday\0\0\0\0Saturday\0\0Sunday\0\0\0\0Monday\0\0\0\0Tuesday"}; + int d = MJD->base_day % 7; + if(d < 0) d += 7; + return &(dow[d*10]); +} + +const char * getMonth( int m) +{ + static char *months = {"Jan\0Feb\0Mar\0Apr\0May\0Jun\0Jul\0Aug\0Sep\0Oct\0Nov\0Dec"}; + return &(months[(m)*4]); +} + +const char * getLongMonth( int m) +{ + static char *months = {"January\0\0\0February\0\0March\0\0\0\0\0April\0\0\0\0\0May\0\0\0\0\0\0\0June\0\0\0\0\0\0July\0\0\0\0\0\0August\0\0\0\0September\0October\0\0\0November\0\0December"}; + return &(months[(m)*10]); +} + + size_t strfMJD(char * buf, size_t len, const char *format, const MJDtime *MJD, int forceJulian) { /* Format a text string according to the format string. This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-10 23:22:09
|
Revision: 9500 http://plplot.svn.sourceforge.net/plplot/?rev=9500&view=rev Author: airwin Date: 2009-02-10 22:05:30 +0000 (Tue, 10 Feb 2009) Log Message: ----------- Logic changes for breakDownMJD. Use combined month logic to simplify code. For Julian case bypass negative j logic (this bypassed logic to be removed later after more testing) and use generalized logic instead for both negative and positive j. The result passes test1 (standard tests for a small range of epochs around the Julian day epoch for both the Julian and Gregorian proleptic calendar cases), but needs more comprehensive testing now. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-10 21:56:29 UTC (rev 9499) +++ trunk/lib/qsastime/qsastime.c 2009-02-10 22:05:30 UTC (rev 9500) @@ -210,7 +210,7 @@ Gregorian proleptic calendar. */ j -= MJD_0000G; - year = (int) ((float)j / 365.2425); + year = (int) ((double)j / 365.2425); lastyear = year - 1; doy = j - year * 365 - lastyear / 4 + lastyear / 100 - lastyear / 400; } @@ -223,7 +223,7 @@ /* Convert MJD struct into date/time elements */ /* Note year 0 CE (AD) [1 BCE (BC)] is a leap year */ - int extra_days,j, year4, year100, year400; + int extra_days,j, year4, year100, year400, ifleapyear; double seconds; if(forceJulian < -1 || forceJulian > 1) { @@ -239,36 +239,20 @@ j = MJD->base_day + extra_days; - if( forceJulian != -1 && j <= MJD_0001J) { + if( 0 && forceJulian != -1 && j < MJD_0001J) { - /* BCE dates */ - j -= MJD_0000J; if( j>=0) { *year = 0; j++; } else { /* negative years */ - *year = (int) ((float)(j-365) / 365.25); - j = j +1 - *year * 365.25; + *year = (int) ((double)(j-365) / 365.25); + j = j + 1 - *year * 365.25; } - + ifleapyear = *year%4 == 0; /* j is now always positive */ - *month = -1; - - if(*year%4 != 0) { - while(j > MonthStartDOY[*month +1]) { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY[*month]; - } else { - while(j > MonthStartDOY_L[*month +1]) { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY_L[*month]; - } + } else if( forceJulian !=-1 && (j < -100840 || forceJulian == 1)) { /* Shift j epoch to 0000-01-01 for the Julian proleptic calendar.*/ j -= MJD_0000J; @@ -278,30 +262,18 @@ i.e. years -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 start with j = -1826 -1461, -1095, -730, -365, 0, 366, 731, 1096, 1461, 1827 */ if(j >= 0) { - *year = (int) ((float)j / 365.25); - year4 = year-1; + *year = (int) ((double)(j) / 365.25); + year4 = *year-1; } else { - *year = (int) ((float)j/ 365.25) - 1; - year4 = year-4; + *year = (int) ((double)(j-365)/ 365.25); + year4 = *year-4; } /* shift j epoch to day of year */ j = j - *year * 365 - year4 / 4; - *month = -1; - if(*year%4 != 0) { - while(j > MonthStartDOY[*month +1]) { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY[*month]; - } else { - while(j > MonthStartDOY_L[*month + 1]) { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY_L[*month]; - } + ifleapyear = *year%4 == 0; + } else { /* forceJulian == -1 || (j >= -100840 && forceJulian == 0) */ /* Shift j epoch to 0000-01-01 for the Gregorian proleptic calendar.*/ @@ -310,12 +282,12 @@ if exactly correct j offset is chosen (see above). */ if(j >=0) { - *year = (int) ((float)j / 365.2425); + *year = (int) ((double)(j) / 365.2425); year4 = *year - 1; year100 = *year - 1; year400 = *year - 1; } else { - *year = (int) ((float)j / 365.2425) - 1; + *year = (int) ((double)(j-365) / 365.2425); year4 = *year - 4; year100 = *year - 100; year400 = *year - 400; @@ -323,23 +295,25 @@ j = j - *year * 365 - year4 / 4 + year100 / 100 - year400 / 400; - *month = -1; - if((*year%4 == 0 && *year%100 != 0) || (*year%4 == 0 && *year%400 == 0) ) { - while(j > MonthStartDOY_L[*month +1]) { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY_L[*month]; - } else { - while(j > MonthStartDOY[*month +1]) { - (*month)++; - if(*month == 11) break; - } - *day = j - MonthStartDOY[*month]; - } - + ifleapyear = (*year%4 == 0 && *year%100 != 0) || (*year%4 == 0 && *year%400 == 0); } + /* calculate month part with j already set to be the day number within + the year in the range from 1 to 366 */ + *month = -1; + if(ifleapyear) { + while(j > MonthStartDOY_L[*month +1]) { + (*month)++; + if(*month == 11) break; + } + *day = j - MonthStartDOY_L[*month]; + } else { + while(j > MonthStartDOY[*month +1]) { + (*month)++; + if(*month == 11) break; + } + *day = j - MonthStartDOY[*month]; + } /* Time part */ seconds = MJD->time_sec - extra_days * SecInDay; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-12 05:21:18
|
Revision: 9513 http://plplot.svn.sourceforge.net/plplot/?rev=9513&view=rev Author: airwin Date: 2009-02-12 05:21:07 +0000 (Thu, 12 Feb 2009) Log Message: ----------- Solve 'year 0' problem in breakDownMJD as revealed by test02 of qsastime_testlib. Subsequently I strengthed the test suite and found a "December 31st" problem, so there is more work to do. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-11 21:07:10 UTC (rev 9512) +++ trunk/lib/qsastime/qsastime.c 2009-02-12 05:21:07 UTC (rev 9513) @@ -261,7 +261,7 @@ day offset is set exactly right so that years -4, 0, 4 are leap years, i.e. years -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 start with j = -1826 -1461, -1095, -730, -365, 0, 366, 731, 1096, 1461, 1827 */ - if(j >= 0) { + if(j >= 366) { *year = (int) ((double)(j) / 365.25); year4 = *year-1; } else { @@ -281,7 +281,7 @@ /* 365.245 is exact period of Gregorian year so year will be correct if exactly correct j offset is chosen (see above). */ - if(j >=0) { + if(j >=366) { *year = (int) ((double)(j) / 365.2425); year4 = *year - 1; year100 = *year - 1; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-13 04:06:17
|
Revision: 9519 http://plplot.svn.sourceforge.net/plplot/?rev=9519&view=rev Author: airwin Date: 2009-02-13 04:06:13 +0000 (Fri, 13 Feb 2009) Log Message: ----------- Apply some end-of-year corrections to the Gregorian calculations for breakDownMJD. This result passes some thorough testing which consists of tests 1 (verbose near JD epoch), 2 (verbose near year epoch), 3 (verbose near MJD epoch), and 4 (non-verbose for every year between -5 million and + 5 million which is very close to the complete possible range since going to +/- 6 million years overflows the double MJD->base_day. For each of these tests, the dates covered for every year are at least Jan 1, Feb 28, Feb 29, Mar 1, and Dec 31 for both the Julian and Gregorian cases. These tests complete all my planned tests for the case where broken-down time is the argument. More tests are planned where a time_t is the argument, although the above tests are already fairly complete. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-13 03:59:29 UTC (rev 9518) +++ trunk/lib/qsastime/qsastime.c 2009-02-13 04:06:13 UTC (rev 9519) @@ -223,7 +223,7 @@ /* Convert MJD struct into date/time elements */ /* Note year 0 CE (AD) [1 BCE (BC)] is a leap year */ - int extra_days,j, year4, year100, year400, ifleapyear; + int extra_days, j, doy, ifcorrect, year4, year100, year400, ifleapyear; double seconds; if(forceJulian < -1 || forceJulian > 1) { @@ -257,10 +257,10 @@ /* Shift j epoch to 0000-01-01 for the Julian proleptic calendar.*/ j -= MJD_0000J; - /* 365.25 is exact period of Julian year so year will be correct if - day offset is set exactly right so that years -4, 0, 4 are leap years, - i.e. years -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 start with j = - -1826 -1461, -1095, -730, -365, 0, 366, 731, 1096, 1461, 1827 */ + /* 365.25 is the exact period of the Julian year so year will be correct + if the day offset is set exactly right so that years -4, 0, 4 are + leap years, i.e. years -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 start with + j = -1826 -1461, -1095, -730, -365, 0, 366, 731, 1096, 1461, 1827 */ if(j >= 366) { *year = (int) ((double)(j) / 365.25); year4 = *year-1; @@ -269,8 +269,7 @@ year4 = *year-4; } - /* shift j epoch to day of year */ - j = j - *year * 365 - year4 / 4; + doy = j - *year * 365 - year4 / 4; ifleapyear = *year%4 == 0; @@ -278,8 +277,10 @@ /* forceJulian == -1 || (j >= -100840 && forceJulian == 0) */ /* Shift j epoch to 0000-01-01 for the Gregorian proleptic calendar.*/ j -= MJD_0000G; - /* 365.245 is exact period of Gregorian year so year will be correct - if exactly correct j offset is chosen (see above). */ + /* 365.245 is the exact period of the Gregorian year so year will be correct + on average, but because the leap year rule is irregular within + the 400-year Gregorian cycle, the first and last days of the + year may need further adjustment. */ if(j >=366) { *year = (int) ((double)(j) / 365.2425); @@ -293,26 +294,50 @@ year400 = *year - 400; } - j = j - *year * 365 - year4 / 4 + year100 / 100 - year400 / 400; - + doy = j - *year * 365 - year4 / 4 + year100 / 100 - year400 / 400; ifleapyear = (*year%4 == 0 && *year%100 != 0) || (*year%4 == 0 && *year%400 == 0); + + /* Rare corrections to above average Gregorian relations. */ + if(doy < 1) { + (*year)--; + ifcorrect = 1; + } else if(doy > 365 && (!ifleapyear || doy > 366)) { + (*year)++; + ifcorrect = 1; + } else { + ifcorrect = 0; + } + if(ifcorrect) { + if(j >=366) { + year4 = *year - 1; + year100 = *year - 1; + year400 = *year - 1; + } else { + year4 = *year - 4; + year100 = *year - 100; + year400 = *year - 400; + } + + doy = j - *year * 365 - year4 / 4 + year100 / 100 - year400 / 400; + ifleapyear = (*year%4 == 0 && *year%100 != 0) || (*year%4 == 0 && *year%400 == 0); + } } - /* calculate month part with j already set to be the day number within + /* calculate month part with doy set to be the day within the year in the range from 1 to 366 */ *month = -1; if(ifleapyear) { - while(j > MonthStartDOY_L[*month +1]) { + while(doy > MonthStartDOY_L[*month +1]) { (*month)++; if(*month == 11) break; } - *day = j - MonthStartDOY_L[*month]; + *day = doy - MonthStartDOY_L[*month]; } else { - while(j > MonthStartDOY[*month +1]) { + while(doy > MonthStartDOY[*month +1]) { (*month)++; if(*month == 11) break; } - *day = j - MonthStartDOY[*month]; + *day = doy - MonthStartDOY[*month]; } /* Time part */ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <sm...@us...> - 2009-02-13 23:22:50
|
Revision: 9527 http://plplot.svn.sourceforge.net/plplot/?rev=9527&view=rev Author: smekal Date: 2009-02-13 22:10:27 +0000 (Fri, 13 Feb 2009) Log Message: ----------- ANSI C only allows to declare a variable at the beginning of a block. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-13 22:06:13 UTC (rev 9526) +++ trunk/lib/qsastime/qsastime.c 2009-02-13 22:10:27 UTC (rev 9527) @@ -684,9 +684,10 @@ } else if(next == 'V') { + int days_in_wk1; /* week of year as a number, (01 - 53) start of week is Monday and first week has at least 3 days in year */ getYAD(&y1, &ifleapyear, &doy, pnMJD, forceJulian); - int days_in_wk1 = (pnMJD->base_day - doy - 3) % 7; + days_in_wk1 = (pnMJD->base_day - doy - 3) % 7; if(days_in_wk1 <= 3) w = (doy +6 - days_in_wk1) / 7; /* ensure first week has at least 3 days in this year */ else w = 1 + (doy + 6 - days_in_wk1) / 7; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-18 06:29:14
|
Revision: 9543 http://plplot.svn.sourceforge.net/plplot/?rev=9543&view=rev Author: airwin Date: 2009-02-18 06:29:13 +0000 (Wed, 18 Feb 2009) Log Message: ----------- Implement public API which consists of configqsas, closeqsas, ctimeqsas, btimeqsas, and strfqsas. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-18 06:26:32 UTC (rev 9542) +++ trunk/lib/qsastime/qsastime.c 2009-02-18 06:29:13 UTC (rev 9543) @@ -38,6 +38,7 @@ */ #include <ctype.h> +#include <math.h> #include "qsastime.h" #include "qsastimeP.h" /* MJD for 0000-01-01 (correctly Jan 01, BCE 1) */ @@ -58,6 +59,111 @@ static const int MonthStartDOY[] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334}; static const int MonthStartDOY_L[] = {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}; +void configqsas(double scale, double offset1, double offset2, int ccontrol) +{ + /* Configure the transformation between continuous time and broken-down time + that is used for ctimeqsas, btimeqsas, and strfqsas. */ + + /* Allocate memory for qsasconfig if that hasn't been done by a + previous call. */ + if(qsasconfig == NULL) { + qsasconfig = (QSASConfig *) malloc((size_t) sizeof(QSASConfig)); + if (qsasconfig == NULL) { + fprintf(stderr, "configqsas: out of memory\n"); + exit(EXIT_FAILURE); + } + } + + if(scale != 0.) { + qsasconfig->scale = scale; + qsasconfig->offset1 = offset1; + qsasconfig->offset2 = offset2; + qsasconfig->ccontrol = ccontrol; + } else { + /* if scale is 0., then use default values. Currently, that + default is continuous time (stored as a double) is seconds since + 1970-01-01 while broken-down time is Gregorian with no other + additional corrections. */ + qsasconfig->scale = 1./86400.; + qsasconfig->offset1 = (double) MJD_1970; + qsasconfig->offset2 = 0.; + qsasconfig->ccontrol = 0x0; + } +} + +void closeqsas(void) +{ + /* Close library if it has been opened. */ + if(qsasconfig != NULL) { + free((void *) qsasconfig); + qsasconfig = NULL; + } +} + +int ctimeqsas(int year, int month, int day, int hour, int min, double sec, double * ctime){ + MJDtime MJD_value, *MJD=&MJD_value; + int forceJulian, ret; + double integral_offset1, integral_offset2, integral_scaled_ctime; + + if(qsasconfig == NULL) { + fprintf(stderr, "libqsastime (ctimeqsas) ERROR: configqsas must be called first.\n"); + exit(EXIT_FAILURE); + } + + if(qsasconfig->ccontrol & 0x1) + forceJulian = 1; + else + forceJulian = -1; + + ret = setFromUT(year, month, day, hour, min, sec, MJD, forceJulian); + if(ret) + return ret; + *ctime = (((double)(MJD->base_day) - qsasconfig->offset1) - qsasconfig->offset2 + MJD->time_sec/(double) SecInDay)/qsasconfig->scale; + return 0; + +} + +void btimeqsas(int *year, int *month, int *day, int *hour, int *min, double *sec, double ctime){ + MJDtime MJD_value, *MJD=&MJD_value; + int forceJulian; + double integral_offset1, integral_offset2, integral_scaled_ctime; + + if(qsasconfig == NULL) { + fprintf(stderr, "libqsastime (btimeqsas) ERROR: configqsas must be called first.\n"); + exit(EXIT_FAILURE); + } + + MJD->time_sec = SecInDay*(modf(qsasconfig->offset1, &integral_offset1) + modf(qsasconfig->offset2, &integral_offset2) + modf(ctime*qsasconfig->scale, &integral_scaled_ctime)); + MJD->base_day = (int) (integral_offset1+integral_offset2+integral_scaled_ctime); + + if(qsasconfig->ccontrol & 0x1) + forceJulian = 1; + else + forceJulian = -1; + + breakDownMJD(year, month, day, hour, min, sec, MJD, forceJulian); +} + +size_t strfqsas(char * buf, size_t len, const char *format, double ctime){ + MJDtime MJD_value, *MJD=&MJD_value; + int forceJulian; + double integral_offset1, integral_offset2, integral_scaled_ctime; + + if(qsasconfig == NULL) { + fprintf(stderr, "libqsastime (strfqsas) ERROR: configqsas must be called first.\n"); + exit(EXIT_FAILURE); + } + MJD->time_sec = SecInDay*(modf(qsasconfig->offset1, &integral_offset1) + modf(qsasconfig->offset2, &integral_offset2) + modf(ctime*qsasconfig->scale, &integral_scaled_ctime)); + MJD->base_day = (int) (integral_offset1+integral_offset2+integral_scaled_ctime); + + if(qsasconfig->ccontrol & 0x1) + forceJulian = 1; + else + forceJulian = -1; + + return strfMJD(buf, len, format, MJD, forceJulian); +} + int setFromUT(int year, int month, int day, int hour, int min, double sec, MJDtime *MJD, int forceJulian) { /* convert Gregorian date plus time to MJD */ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <and...@us...> - 2009-02-18 09:36:09
|
Revision: 9548 http://plplot.svn.sourceforge.net/plplot/?rev=9548&view=rev Author: andrewross Date: 2009-02-18 09:36:04 +0000 (Wed, 18 Feb 2009) Log Message: ----------- Make global variables in qsastime.c const as they are constants. This makes it clear that they are thread-safe. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-18 09:35:46 UTC (rev 9547) +++ trunk/lib/qsastime/qsastime.c 2009-02-18 09:36:04 UTC (rev 9548) @@ -53,8 +53,8 @@ #define MJD_0001J -678577 #define MJD_0001G -678575 -static double SecInDay = 86400; /* we ignore leap seconds */ -static int MJD_1970 = 40587; /* MJD for Jan 01, 1970 00:00:00 */ +static const double SecInDay = 86400; /* we ignore leap seconds */ +static const int MJD_1970 = 40587; /* MJD for Jan 01, 1970 00:00:00 */ static const int MonthStartDOY[] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334}; static const int MonthStartDOY_L[] = {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-23 03:08:35
|
Revision: 9582 http://plplot.svn.sourceforge.net/plplot/?rev=9582&view=rev Author: airwin Date: 2009-02-23 03:08:32 +0000 (Mon, 23 Feb 2009) Log Message: ----------- Initial untested programming attempt for bhunt_search, a binary method of hunting and searching an ordered table. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-22 17:24:44 UTC (rev 9581) +++ trunk/lib/qsastime/qsastime.c 2009-02-23 03:08:32 UTC (rev 9582) @@ -5,6 +5,7 @@ QSAS team, Imperial College, London Copyright (C) 2009 Imperial College, London + Copyright (C) 2009 Alan W. Irwin This file is part of PLplot. @@ -60,6 +61,8 @@ static const int MonthStartDOY[] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334}; static const int MonthStartDOY_L[] = {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}; +int bhunt_search(const void *key, const void *base, size_t n, size_t size, int *index, int (*cmp)(const void *keyval, const void *datum)); + int setFromUT(int year, int month, int day, int hour, int min, double sec, MJDtime *MJD, int forceJulian) { /* convert Gregorian date plus time to MJD */ @@ -963,3 +966,62 @@ return strfMJD(buf, len, format, MJD, forceJulian); } + +/* bhunt_search. Search an ordered table with a binary hunt phase and a + binary search phase. + + On entry *low is used to help the hunt phase speed up the binary + search when consecutive calls to bhunt_search are made with + similar key values. On exit, *low is adjusted such that + base[*(low)-1] < key <= base[*low] with the special cases of + *low set to 0 to indicate the key <= base[0] and *low set to n + to indicate base[n-1] < key. The function *gt must return true (1) + if its first argument (the search key) is greater than + its second argument (a table entry). Otherwise it returns false + (0). Items in the array base must be in ascending order. */ + +int bhunt_search(const void *key, const void *base, size_t n, size_t size, int *low, int (*gt)(const void *keyval, const void *datum)) { + const void *indexbase; + size_t mid, high, hunt_inc=1; + /* Protect against badly defined or undefined *low values. */ + if(*low <= 0 || *low > n) { + *low = 0; + high = n+1; + } else { + /* hunt phase */ + indexbase = (void *) (((const char *) base) + (size*(*low))); + if((*gt)(key, indexbase)) { + high = (*low) + hunt_inc; + indexbase = (void *) (((const char *) base) + (size*high)); + while( (high < n) && ((*gt)(key, indexbase)) ) { + *low = high; + hunt_inc+=hunt_inc; + high = high + hunt_inc; + indexbase = (void *) (((const char *) base) + (size*high)); + } + if(high >= n) + high = n+1; + } else { + high = *low; + *low = high - hunt_inc; + indexbase = (void *) (((const char *) base) + (size*(*low))); + while( ((*low) >= 0) && ((*gt)(key, indexbase)) ) { + high = *low; + hunt_inc+=hunt_inc; + *low = (*low) - hunt_inc; + indexbase = (void *) (((const char *) base) + (size*(*low))); + } + if(*low < 0) + *low = 0; + } + } + /* search phase */ + while(high - *low > 1) { + mid = *low + (high-*low)/2; + indexbase = (void *) (((const char *) base) + (size*mid)); + if((*gt)(key, indexbase)) + *low = mid; + else + high = mid; + } +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-23 16:53:05
|
Revision: 9588 http://plplot.svn.sourceforge.net/plplot/?rev=9588&view=rev Author: airwin Date: 2009-02-23 16:52:54 +0000 (Mon, 23 Feb 2009) Log Message: ----------- Second pass at bhunt_search with errors fixed that were found by visual inspection. The result builds. However, the implementation has not yet been tested to show that it finds the correct interval in the table. Also, the efficiency of the algorithm has not been tested yet (although that should be pretty good if there have been no logic screwups in either the hunt phase or the binary search phase). Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-23 16:00:30 UTC (rev 9587) +++ trunk/lib/qsastime/qsastime.c 2009-02-23 16:52:54 UTC (rev 9588) @@ -61,7 +61,7 @@ static const int MonthStartDOY[] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334}; static const int MonthStartDOY_L[] = {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}; -int bhunt_search(const void *key, const void *base, size_t n, size_t size, int *index, int (*cmp)(const void *keyval, const void *datum)); +int bhunt_search(const void *key, const void *base, size_t n, size_t size, int *low, int (*cmp)(const void *keyval, const void *datum)); int setFromUT(int year, int month, int day, int hour, int min, double sec, MJDtime *MJD, int forceJulian) { @@ -973,53 +973,66 @@ On entry *low is used to help the hunt phase speed up the binary search when consecutive calls to bhunt_search are made with similar key values. On exit, *low is adjusted such that - base[*(low)-1] < key <= base[*low] with the special cases of - *low set to 0 to indicate the key <= base[0] and *low set to n - to indicate base[n-1] < key. The function *gt must return true (1) - if its first argument (the search key) is greater than + base[*low] <= key < base[(*low+1)] with the special cases of + *low set to -1 to indicate the key < base[0] and *low set to n-1 + to indicate base[n-1] <= key. The function *ge must return true (1) + if its first argument (the search key) is greater than or equal to its second argument (a table entry). Otherwise it returns false (0). Items in the array base must be in ascending order. */ -int bhunt_search(const void *key, const void *base, size_t n, size_t size, int *low, int (*gt)(const void *keyval, const void *datum)) { +int bhunt_search(const void *key, const void *base, size_t n, size_t size, int *low, int (*ge)(const void *keyval, const void *datum)) { const void *indexbase; - size_t mid, high, hunt_inc=1; - /* Protect against badly defined or undefined *low values. */ - if(*low <= 0 || *low > n) { + int mid, high, hunt_inc=1; + /* If previous search found below range, then assure one hunt cycle + just in case new key is also below range. */ + if(*low == -1) *low = 0; - high = n+1; + /* Protect against invalid or undefined *low values where hunt + is waste of time. */ + if(*low < 0 || *low >= n) { + *low = -1; + high = n; } else { - /* hunt phase */ + /* binary hunt phase where we are assured 0 <= *low < n */ indexbase = (void *) (((const char *) base) + (size*(*low))); - if((*gt)(key, indexbase)) { - high = (*low) + hunt_inc; + if((*ge)(key, indexbase)) { + high = (*low) + hunt_inc; indexbase = (void *) (((const char *) base) + (size*high)); - while( (high < n) && ((*gt)(key, indexbase)) ) { + /* indexbase is valid if high < n. */ + while( (high < n) && ((*ge)(key, indexbase)) ) { *low = high; hunt_inc+=hunt_inc; high = high + hunt_inc; indexbase = (void *) (((const char *) base) + (size*high)); } if(high >= n) - high = n+1; + high = n; + /* At this point, low is valid and base[low] <= key + and either key < base[high] for valid high or high = n. */ } else { high = *low; *low = high - hunt_inc; indexbase = (void *) (((const char *) base) + (size*(*low))); - while( ((*low) >= 0) && ((*gt)(key, indexbase)) ) { + /* indexbase is valid if(*low) >= 0 */ + while( ((*low) >= 0) && !((*ge)(key, indexbase)) ) { high = *low; hunt_inc+=hunt_inc; *low = (*low) - hunt_inc; indexbase = (void *) (((const char *) base) + (size*(*low))); } - if(*low < 0) - *low = 0; + if((*low) < 0) + *low = -1; + /* At this point high is valid and key < base[high] + and either base[low] <= key for valid low or low = -1. */ } } - /* search phase */ + /* binary search phase where we are assured base[low] <= key < base[high] + when both low and high are valid with obvious special cases signalled + by low = -1 or high = n. */ while(high - *low > 1) { mid = *low + (high-*low)/2; indexbase = (void *) (((const char *) base) + (size*mid)); - if((*gt)(key, indexbase)) + if((*ge)(key, indexbase)) *low = mid; else high = mid; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-24 19:18:47
|
Revision: 9603 http://plplot.svn.sourceforge.net/plplot/?rev=9603&view=rev Author: airwin Date: 2009-02-24 19:18:43 +0000 (Tue, 24 Feb 2009) Log Message: ----------- Tweak name of callback function in formal argument list for bhunt_search declaration. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-24 13:08:40 UTC (rev 9602) +++ trunk/lib/qsastime/qsastime.c 2009-02-24 19:18:43 UTC (rev 9603) @@ -61,7 +61,7 @@ static const int MonthStartDOY[] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334}; static const int MonthStartDOY_L[] = {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}; -int bhunt_search(const void *key, const void *base, size_t n, size_t size, int *low, int (*cmp)(const void *keyval, const void *datum)); +int bhunt_search(const void *key, const void *base, size_t n, size_t size, int *low, int (*ge)(const void *keyval, const void *datum)); int setFromUT(int year, int month, int day, int hour, int min, double sec, MJDtime *MJD, int forceJulian) { This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-26 21:33:49
|
Revision: 9625 http://plplot.svn.sourceforge.net/plplot/?rev=9625&view=rev Author: airwin Date: 2009-02-26 21:33:45 +0000 (Thu, 26 Feb 2009) Log Message: ----------- Fix bugs in strfMJD handling of %. and %N (N= 0-9) format converters. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-26 20:13:40 UTC (rev 9624) +++ trunk/lib/qsastime/qsastime.c 2009-02-26 21:33:45 UTC (rev 9625) @@ -321,12 +321,12 @@ Uses the same syntax as strftime() but does not use current locale. The null terminator is included in len for safety. */ - int year, month, day, hour, min, ysign, sec1, second, d, y; + int year, month, day, hour, min, ysign, second, d, y; int y1, ifleapyear; int i, count,secsSince1970; int nplaces,fmtlen,slen; char * ptr; - double sec,sec_fraction; + double sec, sec_fraction; int w,doy,days_in_wk1; const char *dayText; const char *monthText; @@ -349,8 +349,6 @@ else ysign = 0; second = (int) sec; - sec1 = (int)sec/10; - sec -= (double) sec1*10; /* Read format string, character at a time */ fmtlen = strlen(format); @@ -794,39 +792,31 @@ posn = strlen(buf); if(posn >= last) return posn; } - else if( isdigit(next) != 0 ) + else if( next == '.' || isdigit(next) != 0) { - nplaces = strtol(&(format[i]), NULL, 10); - /* numeric value is number of decimal places ( > 0 ) */ - sec_fraction = sec - (double) second; + /* nplaces is number of decimal places ( 0 < nplaces <= 9 ) */ + if( next == '.' ) + /* maximum number of places*/ + nplaces = 9; + else + nplaces = strtol(&(format[i]), NULL, 10); - for(count=0; count<nplaces; count++) sec_fraction *= 10; - sprintf(DateTime, ".%d", (int) sec_fraction); - - /* append 0 to pad to length */ - slen = strlen(DateTime); - while(slen < nplaces+1) - { - DateTime[slen] = '0'; - slen++; - DateTime[slen] = '\0'; - } - strncat(&(buf[posn]), DateTime, last - posn); - posn = strlen(buf); - if(posn >= last) return posn; - } - else if( next == '.' ) - { /* fractional part of seconds to maximum available accuracy */ - sec_fraction = sec - (double) second; - sprintf(DateTime, "%-11.9f", sec_fraction); + sec_fraction = sec - (int) sec; + char dynamic_format[10]; + sprintf(dynamic_format, "%%-%d.%df", nplaces+2, nplaces); + /*sprintf(DateTime, "%-11.9f", sec_fraction);*/ + sprintf(DateTime, dynamic_format, sec_fraction); while( ( ptr = strrchr(&(DateTime[0]), ' ')) != NULL) ptr[0] ='\0'; /* remove trailing white space */ - slen = strlen(DateTime) -1; - while( DateTime[slen] == '0' && DateTime[slen-1] != '.') { - DateTime[slen] ='\0'; /* remove trailing zeros */ - slen --; + + if( next == '.' ) { + slen = strlen(DateTime) -1; + while( DateTime[slen] == '0' && DateTime[slen-1] != '.') { + DateTime[slen] ='\0'; /* remove trailing zeros */ + slen --; + } } - + ptr = strchr(DateTime, '.'); /* remove possible lead 0 */ strncat(&(buf[posn]), ptr, last - posn); posn = strlen(buf); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <sm...@us...> - 2009-02-27 21:28:17
|
Revision: 9634 http://plplot.svn.sourceforge.net/plplot/?rev=9634&view=rev Author: smekal Date: 2009-02-27 21:28:16 +0000 (Fri, 27 Feb 2009) Log Message: ----------- ANSI C only allows to declare variables at the top of a block. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-27 21:25:39 UTC (rev 9633) +++ trunk/lib/qsastime/qsastime.c 2009-02-27 21:28:16 UTC (rev 9634) @@ -794,6 +794,8 @@ } else if( next == '.' || isdigit(next) != 0) { + char dynamic_format[10]; + /* nplaces is number of decimal places ( 0 < nplaces <= 9 ) */ if( next == '.' ) /* maximum number of places*/ @@ -803,7 +805,6 @@ /* fractional part of seconds to maximum available accuracy */ sec_fraction = sec - (int) sec; - char dynamic_format[10]; sprintf(dynamic_format, "%%-%d.%df", nplaces+2, nplaces); /*sprintf(DateTime, "%-11.9f", sec_fraction);*/ sprintf(DateTime, dynamic_format, sec_fraction); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-28 03:11:56
|
Revision: 9641 http://plplot.svn.sourceforge.net/plplot/?rev=9641&view=rev Author: airwin Date: 2009-02-28 03:11:55 +0000 (Sat, 28 Feb 2009) Log Message: ----------- Apply Tony Allen's patch to properly normalize strfMJD results using pre-rounding. The implementation uses the most lightly rounded seconds value in the format string to set the degree of rounding used. This change should yield proper normalization of the strfMJD formatted results in all cases. To give an example, that is tested with Test 05 of qsastime_testlib, if H:M:S = 0:59:59.99999999000, internally, that will format as a properly normalized 1:00:00.0000 if sufficient rounding (e.g., %S%4 which corresponds to rounding errors of 50 ms) is used, and will format as a properly normalized 0:59:59.999999990 if light rounding (e.g., %S%9 which corresponds to rounding errors of 0.5 ns) is used. In the latter case, just like in any floating-point calculation, reducing the rounding means you will see more details, but those details potentially include numerical errors. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-28 02:38:37 UTC (rev 9640) +++ trunk/lib/qsastime/qsastime.c 2009-02-28 03:11:55 UTC (rev 9641) @@ -325,6 +325,8 @@ int y1, ifleapyear; int i, count,secsSince1970; int nplaces,fmtlen,slen; + int resolution; + double shiftPlaces; char * ptr; double sec, sec_fraction; int w,doy,days_in_wk1; @@ -334,10 +336,41 @@ size_t posn = 0; size_t last = len -1; MJDtime nMJD, *pnMJD=&nMJD; + MJDtime roundMJD; char dynamic_format[10]; - normalize_MJD(pnMJD, MJD); + /* Find required resolution */ + resolution = 0; + fmtlen = strlen(format); + i=0; + while(i<fmtlen) + { + char next = format[i]; + if( next == '%') + { + /* find seconds format if used */ + i++; + next = format[i]; + if( isdigit(next) != 0 ) + { + nplaces = strtol(&(format[i]), NULL, 10); + if(nplaces > resolution) resolution = nplaces; + } + else if( next == '.' ) + { + resolution = 9; /* maximum resolution allowed */ + } + } + i++; + } + /* ensure rounding is done before breakdown */ + shiftPlaces = pow(10,(double)resolution); + roundMJD = *MJD; + roundMJD.time_sec += 0.5/shiftPlaces; + + normalize_MJD(pnMJD, &roundMJD); + buf[last] = '\0'; buf[0] = '\0'; /* force overwrite of old buffer since strnctat() used hereafter */ @@ -349,10 +382,11 @@ } else ysign = 0; + /*truncate seconds to resolution to stop formatting rounding up */ + sec = floor(sec *shiftPlaces) / shiftPlaces; second = (int) sec; /* Read format string, character at a time */ - fmtlen = strlen(format); i=0; while(i<fmtlen) { This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-02-28 03:22:09
|
Revision: 9638 http://plplot.svn.sourceforge.net/plplot/?rev=9638&view=rev Author: airwin Date: 2009-02-28 02:06:59 +0000 (Sat, 28 Feb 2009) Log Message: ----------- Sort out a number of formatting issues for seconds. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-02-27 22:17:00 UTC (rev 9637) +++ trunk/lib/qsastime/qsastime.c 2009-02-28 02:06:59 UTC (rev 9638) @@ -334,6 +334,7 @@ size_t posn = 0; size_t last = len -1; MJDtime nMJD, *pnMJD=&nMJD; + char dynamic_format[10]; normalize_MJD(pnMJD, MJD); @@ -586,8 +587,32 @@ } else if(next == 'S') { - /* second (00 - 59) */ - sprintf(DateTime, "%02d", second); + /* second (00 - 59 with optional decimal point and numbers after the decimal point.) */ + if(i+2 < fmtlen && format[i+1] == '%' && (format[i+2] == '.' || isdigit(format[i+2]) != 0)) { + /* nplaces is number of decimal places ( 0 < nplaces <= 9 ) */ + if( format[i+2] == '.' ) + /* maximum number of places*/ + nplaces = 9; + else + nplaces = strtol(&(format[i+2]), NULL, 10); + i+=2; + } else { + nplaces = 0; + } + + if(nplaces == 0) { + sprintf(DateTime, "%02d", (int) (sec + 0.5)); + }else { + sprintf(dynamic_format, "%%0%d.%df", nplaces+3, nplaces); + sprintf(DateTime, dynamic_format, sec); + if( format[i] == '.' ) { + slen = strlen(DateTime) -1; + while( DateTime[slen] == '0' && DateTime[slen-1] != '.') { + DateTime[slen] ='\0'; /* remove trailing zeros */ + slen --; + } + } + } strncat(&(buf[posn]), DateTime, last - posn); posn = strlen(buf); @@ -794,8 +819,6 @@ } else if( next == '.' || isdigit(next) != 0) { - char dynamic_format[10]; - /* nplaces is number of decimal places ( 0 < nplaces <= 9 ) */ if( next == '.' ) /* maximum number of places*/ @@ -818,8 +841,11 @@ } } - ptr = strchr(DateTime, '.'); /* remove possible lead 0 */ - strncat(&(buf[posn]), ptr, last - posn); + ptr = strchr(DateTime, '.'); + /* remove everything in front of the decimal point, and + ignore case (%0) where no decimal point exists at all. */ + if(ptr != NULL) + strncat(&(buf[posn]), ptr, last - posn); posn = strlen(buf); if(posn >= last) return posn; } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-03-01 07:04:19
|
Revision: 9647 http://plplot.svn.sourceforge.net/plplot/?rev=9647&view=rev Author: airwin Date: 2009-03-01 07:04:05 +0000 (Sun, 01 Mar 2009) Log Message: ----------- Implement leap-second option for btimeqsas. Details for last second before discontinuity need more effort to give broken-down time with (int) sec = 60 to flag such epochs where an extra amount of time is being inserted. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-03-01 06:55:15 UTC (rev 9646) +++ trunk/lib/qsastime/qsastime.c 2009-03-01 07:04:05 UTC (rev 9647) @@ -41,6 +41,7 @@ #include <ctype.h> #include <math.h> #include "qsastimeP.h" +#include "tai-utc.h" /* MJD for 0000-01-01 (correctly Jan 01, BCE 1) */ /* Julian proleptic calendar value. */ @@ -61,6 +62,7 @@ static const int MonthStartDOY[] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334}; static const int MonthStartDOY_L[] = {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}; +static int geMJDtimeTAI_UTC(const MJDtime *number1, const TAI_UTC *number2); int bhunt_search(const void *key, const void *base, size_t n, size_t size, int *low, int (*ge)(const void *keyval, const void *datum)); int setFromUT(int year, int month, int day, int hour, int min, double sec, MJDtime *MJD, int forceJulian) @@ -897,6 +899,18 @@ return posn; } +int geMJDtimeTAI_UTC(const MJDtime *number1, const TAI_UTC *number2) { + /* Returns true if number1 >= number2. */ + /* N.B. both number1 and number2 must be normalized. */ + if(number1->base_day > number2->base_day) { + return 1; + } else if (number1->base_day < number2->base_day) { + return 0; + } else { + return (number1->time_sec >= number2->time_sec); + } +} + void configqsas(double scale, double offset1, double offset2, int ccontrol, int ifbtime_offset, int year, int month, int day, int hour, int min, double sec, QSASConfig **qsasconfig) { /* Configure the transformation between continuous time and broken-down time @@ -977,9 +991,10 @@ } void btimeqsas(int *year, int *month, int *day, int *hour, int *min, double *sec, double ctime, const QSASConfig *qsasconfig){ - MJDtime MJD_value, *MJD=&MJD_value; + MJDtime MJD_value, *MJD=&MJD_value, MJD1_value; int forceJulian; double integral_offset1, integral_offset2, integral_scaled_ctime; + int index; if(qsasconfig == NULL) { fprintf(stderr, "libqsastime (btimeqsas) ERROR: configqsas must be called first.\n"); @@ -994,6 +1009,51 @@ else forceJulian = 0; + if(qsasconfig->ccontrol & 0x2) { + /* leap-second logic assuming the above offsets have converted time argument to MJD(TAI). */ + /* N.B. geMJDtimeTAI_UTC only works for normalized values. */ + MJD1_value = MJD_value; + normalize_MJD(MJD, &MJD1_value); + /* Search for index such that TAI_TO_UTC_lookup_table[index] <= MJD(TAI) < TAI_TO_UTC_lookup_table[index+1] */ + bhunt_search(MJD, TAI_TO_UTC_lookup_table, number_of_entries_in_tai_utc_table, sizeof(TAI_UTC), &index, (int (*)(const void *, const void *)) geMJDtimeTAI_UTC); + if(index == -1) { + /* MJD is less than first table entry. */ + /* Debug: check that condition is met */ + if(geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index+1])) { + fprintf(stderr, "libqsastime (btimeqsas) logic ERROR: bad condition for index = %d\n", index); + exit(EXIT_FAILURE); + } + /* Use initial offset for MJD values before first table entry. */ + /* Calculate this offset strictly from offset1. The slope term + doesn't enter because offset2 is the same as the UTC of the + first epoch of the table. */ + MJD->time_sec -= TAI_TO_UTC_lookup_table[index+1].offset1; + } else if(index == number_of_entries_in_tai_utc_table-1) { + /* MJD is greater than or equal to last table entry. */ + /* Debug: check that condition is met */ + if(!geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index])) { + fprintf(stderr, "libqsastime (btimeqsas) logic ERROR: bad condition for index = %d\n", index); + exit(EXIT_FAILURE); + } + /* Use final offset for MJD values after last table entry. + The slope term doesn't enter because modern values of the slope + are zero.*/ + MJD->time_sec -= TAI_TO_UTC_lookup_table[index].offset1; + } else if(index >= 0 && index < number_of_entries_in_tai_utc_table) { + /* table[index] <= MJD < table[index+1]. */ + /* Debug: check that condition is met */ + if(!(geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index]) && !geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index+1]))) { + fprintf(stderr, "MJD = {%d, %f}\n", MJD->base_day, MJD->time_sec); + fprintf(stderr, "libqsastime (btimeqsas) logic ERROR: bad condition for index = %d\n", index); + exit(EXIT_FAILURE); + } + MJD->time_sec -= (TAI_TO_UTC_lookup_table[index].offset1 + ((MJD->base_day-TAI_TO_UTC_lookup_table[index].offset2) + MJD->time_sec/SecInDay)*TAI_TO_UTC_lookup_table[index].slope)/(1. + TAI_TO_UTC_lookup_table[index].slope/SecInDay); + } else { + fprintf(stderr, "libqsastime (btimeqsas) logic ERROR: bad index = %d\n", index); + exit(EXIT_FAILURE); + } + } + breakDownMJD(year, month, day, hour, min, sec, MJD, forceJulian); } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-03-01 23:36:46
|
Revision: 9650 http://plplot.svn.sourceforge.net/plplot/?rev=9650&view=rev Author: airwin Date: 2009-03-01 23:36:43 +0000 (Sun, 01 Mar 2009) Log Message: ----------- Split off leap-second calculations for MJD(TAI) into its own little static function, leap_second_TAI. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-03-01 19:51:59 UTC (rev 9649) +++ trunk/lib/qsastime/qsastime.c 2009-03-01 23:36:43 UTC (rev 9650) @@ -63,6 +63,7 @@ static const int MonthStartDOY_L[] = {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}; static int geMJDtimeTAI_UTC(const MJDtime *number1, const TAI_UTC *number2); +static double leap_second_TAI(const MJDtime *MJD_TAI); int bhunt_search(const void *key, const void *base, size_t n, size_t size, int *low, int (*ge)(const void *keyval, const void *datum)); int setFromUT(int year, int month, int day, int hour, int min, double sec, MJDtime *MJD, int forceJulian) @@ -911,6 +912,54 @@ } } +double leap_second_TAI(const MJDtime *MJD_TAI) { + /* Logic assumes input MJD_TAI is in TAI */ + + MJDtime MJD_value, *MJD = &MJD_value; + int index; + + /* N.B. geMJDtimeTAI_UTC only works for normalized values. */ + normalize_MJD(MJD, MJD_TAI); + /* Search for index such that TAI_TO_UTC_lookup_table[index] <= MJD(TAI) < TAI_TO_UTC_lookup_table[index+1] */ + bhunt_search(MJD, TAI_TO_UTC_lookup_table, number_of_entries_in_tai_utc_table, sizeof(TAI_UTC), &index, (int (*)(const void *, const void *)) geMJDtimeTAI_UTC); + if(index == -1) { + /* MJD is less than first table entry. */ + /* Debug: check that condition is met */ + if(geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index+1])) { + fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for index = %d\n", index); + exit(EXIT_FAILURE); + } + /* Use initial offset for MJD values before first table entry. */ + /* Calculate this offset strictly from offset1. The slope term + doesn't enter because offset2 is the same as the UTC of the + first epoch of the table. */ + return -TAI_TO_UTC_lookup_table[index+1].offset1; + } else if(index == number_of_entries_in_tai_utc_table-1) { + /* MJD is greater than or equal to last table entry. */ + /* Debug: check that condition is met */ + if(!geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index])) { + fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for index = %d\n", index); + exit(EXIT_FAILURE); + } + /* Use final offset for MJD values after last table entry. + The slope term doesn't enter because modern values of the slope + are zero.*/ + return -TAI_TO_UTC_lookup_table[index].offset1; + } else if(index >= 0 && index < number_of_entries_in_tai_utc_table) { + /* table[index] <= MJD < table[index+1]. */ + /* Debug: check that condition is met */ + if(!(geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index]) && !geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index+1]))) { + fprintf(stderr, "MJD = {%d, %f}\n", MJD->base_day, MJD->time_sec); + fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for index = %d\n", index); + exit(EXIT_FAILURE); + } + return -(TAI_TO_UTC_lookup_table[index].offset1 + ((MJD->base_day-TAI_TO_UTC_lookup_table[index].offset2) + MJD->time_sec/SecInDay)*TAI_TO_UTC_lookup_table[index].slope)/(1. + TAI_TO_UTC_lookup_table[index].slope/SecInDay); + } else { + fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad index = %d\n", index); + exit(EXIT_FAILURE); + } +} + void configqsas(double scale, double offset1, double offset2, int ccontrol, int ifbtime_offset, int year, int month, int day, int hour, int min, double sec, QSASConfig **qsasconfig) { /* Configure the transformation between continuous time and broken-down time @@ -991,10 +1040,9 @@ } void btimeqsas(int *year, int *month, int *day, int *hour, int *min, double *sec, double ctime, const QSASConfig *qsasconfig){ - MJDtime MJD_value, *MJD=&MJD_value, MJD1_value; + MJDtime MJD_value, *MJD=&MJD_value; int forceJulian; double integral_offset1, integral_offset2, integral_scaled_ctime; - int index; if(qsasconfig == NULL) { fprintf(stderr, "libqsastime (btimeqsas) ERROR: configqsas must be called first.\n"); @@ -1010,48 +1058,7 @@ forceJulian = 0; if(qsasconfig->ccontrol & 0x2) { - /* leap-second logic assuming the above offsets have converted time argument to MJD(TAI). */ - /* N.B. geMJDtimeTAI_UTC only works for normalized values. */ - MJD1_value = MJD_value; - normalize_MJD(MJD, &MJD1_value); - /* Search for index such that TAI_TO_UTC_lookup_table[index] <= MJD(TAI) < TAI_TO_UTC_lookup_table[index+1] */ - bhunt_search(MJD, TAI_TO_UTC_lookup_table, number_of_entries_in_tai_utc_table, sizeof(TAI_UTC), &index, (int (*)(const void *, const void *)) geMJDtimeTAI_UTC); - if(index == -1) { - /* MJD is less than first table entry. */ - /* Debug: check that condition is met */ - if(geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index+1])) { - fprintf(stderr, "libqsastime (btimeqsas) logic ERROR: bad condition for index = %d\n", index); - exit(EXIT_FAILURE); - } - /* Use initial offset for MJD values before first table entry. */ - /* Calculate this offset strictly from offset1. The slope term - doesn't enter because offset2 is the same as the UTC of the - first epoch of the table. */ - MJD->time_sec -= TAI_TO_UTC_lookup_table[index+1].offset1; - } else if(index == number_of_entries_in_tai_utc_table-1) { - /* MJD is greater than or equal to last table entry. */ - /* Debug: check that condition is met */ - if(!geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index])) { - fprintf(stderr, "libqsastime (btimeqsas) logic ERROR: bad condition for index = %d\n", index); - exit(EXIT_FAILURE); - } - /* Use final offset for MJD values after last table entry. - The slope term doesn't enter because modern values of the slope - are zero.*/ - MJD->time_sec -= TAI_TO_UTC_lookup_table[index].offset1; - } else if(index >= 0 && index < number_of_entries_in_tai_utc_table) { - /* table[index] <= MJD < table[index+1]. */ - /* Debug: check that condition is met */ - if(!(geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index]) && !geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[index+1]))) { - fprintf(stderr, "MJD = {%d, %f}\n", MJD->base_day, MJD->time_sec); - fprintf(stderr, "libqsastime (btimeqsas) logic ERROR: bad condition for index = %d\n", index); - exit(EXIT_FAILURE); - } - MJD->time_sec -= (TAI_TO_UTC_lookup_table[index].offset1 + ((MJD->base_day-TAI_TO_UTC_lookup_table[index].offset2) + MJD->time_sec/SecInDay)*TAI_TO_UTC_lookup_table[index].slope)/(1. + TAI_TO_UTC_lookup_table[index].slope/SecInDay); - } else { - fprintf(stderr, "libqsastime (btimeqsas) logic ERROR: bad index = %d\n", index); - exit(EXIT_FAILURE); - } + MJD->time_sec += leap_second_TAI(MJD); } breakDownMJD(year, month, day, hour, min, sec, MJD, forceJulian); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-03-02 23:22:59
|
Revision: 9660 http://plplot.svn.sourceforge.net/plplot/?rev=9660&view=rev Author: airwin Date: 2009-03-02 23:22:47 +0000 (Mon, 02 Mar 2009) Log Message: ----------- Add int if60secformat to strfMJD argument list. This forces renormalization of results so that the integer part of the second is 60 to signal that the epoch is in the middle of leap-interval insertion just prior to a positive discontinuity in TAI-UTC. The implementation of the actual renormalization is still in progress, but I am making the commit now because this builds, and all the other changes in this source file have been thoroughly tested and do work. Return inleap value from leap_second_TAI to signal when the epoch is in the middle of a leap-interval insertion just prior to a positive discontinuity in TAI-UTC. This inleap logic should be complete, but needs testing when the if60secformat implementation is complete for strfMJD. bhunt_search declaration moved to qsastimeP.h. Simplify normalize_MJD argument list. Drop "p" form of MJDtime variable name. Instead, use _value prefix on value form of variable and raw name for pointer to that value. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-03-02 23:02:37 UTC (rev 9659) +++ trunk/lib/qsastime/qsastime.c 2009-03-02 23:22:47 UTC (rev 9660) @@ -62,18 +62,22 @@ static const int MonthStartDOY[] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334}; static const int MonthStartDOY_L[] = {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}; +/* Static function declarations. */ static int geMJDtimeTAI_UTC(const MJDtime *number1, const TAI_UTC *number2); -static double leap_second_TAI(const MJDtime *MJD_TAI); -int bhunt_search(const void *key, const void *base, size_t n, size_t size, int *low, int (*ge)(const void *keyval, const void *datum)); +static double leap_second_TAI(const MJDtime *MJD_TAI, int * inleap); +/* End of static function declarations. */ int setFromUT(int year, int month, int day, int hour, int min, double sec, MJDtime *MJD, int forceJulian) { - /* convert Gregorian date plus time to MJD */ + /* convert broken-down time to MJD */ /* MJD measures from the start of 17 Nov 1858 */ - /* the int flag Julian forces use of Julian calendar whatever the year */ - /* default is to use Gregorian after 4 Oct 1582 (Julian) i.e. from 15 Oct 1582 Gregorian */ - /* Note C libraries use Gregorian only from 14 Sept 1752 onwards */ + /* If forceJulian is true (non-zero), the Julian proleptic calendar is + used whatever the year. Otherwise, the Gregorian proleptic calendar + is used whatever the year. */ + /* Note C libraries use Gregorian only (at least on Linux). In contrast, + the Linux (and Unix?) cal application uses Julian for earlier dates + and Gregorian from 14 Sept 1752 onwards. */ int leaps, year4, year100, year400; double dbase_day, non_leaps = 365.; @@ -233,23 +237,22 @@ *ifleapyear = (*year%4 == 0 && *year%100 != 0) || (*year%4 == 0 && *year%400 == 0); } } - } -void normalize_MJD(MJDtime *MJDout, const MJDtime *MJDin) +void normalize_MJD(MJDtime *MJD) { int extra_days; /* Calculate MJDout as normalized version (i.e., 0. <= MJDout->time_sec < 86400.) of MJDin. */ - if(MJDin->time_sec >= 0) { - extra_days = (int) (MJDin->time_sec / SecInDay); + if(MJD->time_sec >= 0) { + extra_days = (int) (MJD->time_sec / SecInDay); } else { /* allow for negative seconds push into previous day even if less than 1 day */ - extra_days = (int) (MJDin->time_sec / SecInDay) - 1 ; + extra_days = (int) (MJD->time_sec / SecInDay) - 1 ; } - MJDout->base_day = MJDin->base_day + extra_days; - MJDout->time_sec = MJDin->time_sec - extra_days * SecInDay; + MJD->base_day += extra_days; + MJD->time_sec -= extra_days * SecInDay; } void breakDownMJD(int *year, int *month, int *day, int *hour, int *min, double *sec, const MJDtime *MJD, int forceJulian) @@ -258,18 +261,19 @@ /* Note year 0 CE (AD) [1 BCE (BC)] is a leap year */ int doy, ifleapyear; - MJDtime nMJD, *pnMJD=&nMJD; + MJDtime nMJD_value, *nMJD=&nMJD_value; - normalize_MJD(pnMJD, MJD); + *nMJD = *MJD; + normalize_MJD(nMJD); /* Time part */ - *sec = pnMJD->time_sec; + *sec = nMJD->time_sec; *hour = (int)( *sec / 3600.); *sec -= (double) *hour * 3600.; *min = (int) ( *sec / 60.); *sec -= (double) *min * 60.; - getYAD(year, &ifleapyear, &doy, pnMJD, forceJulian); + getYAD(year, &ifleapyear, &doy, nMJD, forceJulian); /* calculate month part with doy set to be the day within the year in the range from 1 to 366 */ @@ -318,11 +322,19 @@ } -size_t strfMJD(char * buf, size_t len, const char *format, const MJDtime *MJD, int forceJulian) +size_t strfMJD(char * buf, size_t len, const char *format, const MJDtime *MJD, int forceJulian, int if60secformat) { /* Format a text string according to the format string. Uses the same syntax as strftime() but does not use current locale. The null terminator is included in len for safety. */ + + /* if if60secformat is true (non-zero) AND the time is somewhere in the + first second of a year, then renormalize to the previous year so that the + date/time comes out beyond the normal last second of the year with the + seconds greater than 60. This form of normalization is used as a flag + that a leap increment (recently always a second, but historically it + was sometimes smaller than that) was in the process of being inserted + for this particular epoch just prior to a discontinuity in TAI-UTC. */ int year, month, day, hour, min, ysign, second, d, y; int y1, ifleapyear; @@ -338,8 +350,7 @@ char DateTime[80]; size_t posn = 0; size_t last = len -1; - MJDtime nMJD, *pnMJD=&nMJD; - MJDtime roundMJD; + MJDtime nMJD_value, *nMJD=&nMJD_value; char dynamic_format[10]; /* Find required resolution */ @@ -369,15 +380,23 @@ /* ensure rounding is done before breakdown */ shiftPlaces = pow(10,(double)resolution); - roundMJD = *MJD; - roundMJD.time_sec += 0.5/shiftPlaces; + *nMJD = *MJD; + nMJD->time_sec += 0.5/shiftPlaces; - normalize_MJD(pnMJD, &roundMJD); + normalize_MJD(nMJD); buf[last] = '\0'; buf[0] = '\0'; /* force overwrite of old buffer since strnctat() used hereafter */ - breakDownMJD(&year, &month, &day, &hour, &min, &sec, pnMJD, forceJulian); + breakDownMJD(&year, &month, &day, &hour, &min, &sec, nMJD, forceJulian); + if(if60secformat && month == 0 && day == 1 && hour == 0 && min == 0 && (int) sec == 0) { + year = year-1; + month = 11; + day = 31; + hour = 23; + min = 59; + sec+= 60.; + } if(year < 0) { ysign = 1; @@ -409,7 +428,7 @@ else if(next == 'a') { /* short day name */ - dayText = getDayOfWeek(pnMJD); + dayText = getDayOfWeek(nMJD); strncat(&(buf[posn]), dayText, last - posn); posn = strlen(buf); if(posn >= last) return posn; @@ -417,7 +436,7 @@ else if(next == 'A') { /* long day name */ - dayText = getLongDayOfWeek(pnMJD); + dayText = getLongDayOfWeek(nMJD); strncat(&(buf[posn]), dayText, last - posn); posn = strlen(buf); if(posn >= last) return posn; @@ -441,7 +460,7 @@ else if(next == 'c') { /* Date and Time with day of week */ - dayText = getDayOfWeek(pnMJD); + dayText = getDayOfWeek(nMJD); monthText = getMonth(month); if(ysign == 0) sprintf(DateTime, "%s %s %02d %02d:%02d:%02d %04d", dayText, monthText, day, hour, min, second, year ); @@ -534,7 +553,7 @@ else if(next == 'j') { /* day of year */ - getYAD(&y1, &ifleapyear, &doy, pnMJD, forceJulian); + getYAD(&y1, &ifleapyear, &doy, nMJD, forceJulian); sprintf(DateTime, "%03d", doy); strncat(&(buf[posn]), DateTime, last - posn); @@ -658,7 +677,7 @@ else if(next == 's') { /* seconds since 01 Jan 1970 Gregorian */ - secsSince1970 = (int)(pnMJD->time_sec + (pnMJD->base_day - MJD_1970) * SecInDay); + secsSince1970 = (int)(nMJD->time_sec + (nMJD->base_day - MJD_1970) * SecInDay); sprintf(DateTime, "%d", secsSince1970); strncat(&(buf[posn]), DateTime, last - posn); @@ -684,8 +703,8 @@ else if(next == 'U') { /* week of year as a number, (00 - 53) start of week is Sunday */ - getYAD(&y1, &ifleapyear, &doy, pnMJD, forceJulian); - days_in_wk1 = (pnMJD->base_day - doy - 4) % 7; + getYAD(&y1, &ifleapyear, &doy, nMJD, forceJulian); + days_in_wk1 = (nMJD->base_day - doy - 4) % 7; w = (doy + 6 - days_in_wk1) / 7; @@ -698,7 +717,7 @@ else if(next == 'u') { /* weekday as a number, 0 = Monday */ - d = 1 + (pnMJD->base_day - 5) % 7; + d = 1 + (nMJD->base_day - 5) % 7; sprintf(DateTime, "%01d", d); @@ -735,8 +754,8 @@ { int days_in_wk1; /* week of year as a number, (01 - 53) start of week is Monday and first week has at least 3 days in year */ - getYAD(&y1, &ifleapyear, &doy, pnMJD, forceJulian); - days_in_wk1 = (pnMJD->base_day - doy - 3) % 7; + getYAD(&y1, &ifleapyear, &doy, nMJD, forceJulian); + days_in_wk1 = (nMJD->base_day - doy - 3) % 7; if(days_in_wk1 <= 3) w = (doy +6 - days_in_wk1) / 7; /* ensure first week has at least 3 days in this year */ else w = 1 + (doy + 6 - days_in_wk1) / 7; @@ -751,7 +770,7 @@ else if(next == 'w') { /* weekday as a number, 0 = Sunday */ - d = (pnMJD->base_day - 4) % 7; + d = (nMJD->base_day - 4) % 7; sprintf(DateTime, "%01d", d); @@ -762,8 +781,8 @@ else if(next == 'W') { /* week of year as a number, (00 - 53) start of week is Monday */ - getYAD(&y1, &ifleapyear, &doy, pnMJD, forceJulian); - days_in_wk1 = (pnMJD->base_day - doy - 3) % 7; + getYAD(&y1, &ifleapyear, &doy, nMJD, forceJulian); + days_in_wk1 = (nMJD->base_day - doy - 3) % 7; w = (doy +6 - days_in_wk1) / 7; @@ -776,7 +795,7 @@ else if(next == 'x') { /* date string */ - dayText = getDayOfWeek(pnMJD); + dayText = getDayOfWeek(nMJD); monthText = getMonth(month); if(ysign == 0) sprintf(DateTime, "%s %s %02d, %04d", dayText, monthText, day, year ); @@ -843,7 +862,7 @@ else if(next == '+') { /* date and time */ - dayText = getDayOfWeek(pnMJD); + dayText = getDayOfWeek(nMJD); monthText = getMonth(month); if(ysign == 0) sprintf(DateTime, "%s %s %02d %02d:%02d:%02d UTC %04d", dayText, monthText, day, hour, min, second, year ); @@ -912,14 +931,17 @@ } } -double leap_second_TAI(const MJDtime *MJD_TAI) { +double leap_second_TAI(const MJDtime *MJD_TAI, int * inleap) { /* Logic assumes input MJD_TAI is in TAI */ + /* *inleap lets the calling routine know whether MJD_TAI corresponds + to an epoch when a positive leap increment is being inserted. */ MJDtime MJD_value, *MJD = &MJD_value; int index; - + double leap; /* N.B. geMJDtimeTAI_UTC only works for normalized values. */ - normalize_MJD(MJD, MJD_TAI); + *MJD = *MJD_TAI; + normalize_MJD(MJD); /* Search for index such that TAI_TO_UTC_lookup_table[index] <= MJD(TAI) < TAI_TO_UTC_lookup_table[index+1] */ bhunt_search(MJD, TAI_TO_UTC_lookup_table, number_of_entries_in_tai_utc_table, sizeof(TAI_UTC), &index, (int (*)(const void *, const void *)) geMJDtimeTAI_UTC); if(index == -1) { @@ -929,6 +951,9 @@ fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for index = %d\n", index); exit(EXIT_FAILURE); } + /* There is (by assertion) no discontinuity at the start of the table. + Therefore, *inleap cannot be true. */ + *inleap = 0; /* Use initial offset for MJD values before first table entry. */ /* Calculate this offset strictly from offset1. The slope term doesn't enter because offset2 is the same as the UTC of the @@ -941,6 +966,8 @@ fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for index = %d\n", index); exit(EXIT_FAILURE); } + /* If beyond end of table, cannot be in middle of leap second insertion. */ + *inleap = 0; /* Use final offset for MJD values after last table entry. The slope term doesn't enter because modern values of the slope are zero.*/ @@ -953,7 +980,18 @@ fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for index = %d\n", index); exit(EXIT_FAILURE); } - return -(TAI_TO_UTC_lookup_table[index].offset1 + ((MJD->base_day-TAI_TO_UTC_lookup_table[index].offset2) + MJD->time_sec/SecInDay)*TAI_TO_UTC_lookup_table[index].slope)/(1. + TAI_TO_UTC_lookup_table[index].slope/SecInDay); + leap = -(TAI_TO_UTC_lookup_table[index].offset1 + ((MJD->base_day-TAI_TO_UTC_lookup_table[index].offset2) + MJD->time_sec/SecInDay)*TAI_TO_UTC_lookup_table[index].slope)/(1. + TAI_TO_UTC_lookup_table[index].slope/SecInDay); + /* Convert MJD(TAI) to normalized MJD(UTC). */ + MJD->time_sec += leap; + normalize_MJD(MJD); + /* If MJD(UT) is in the next interval of the corresponding + UTC_TO_TAI_lookup_table, then we are right in the middle of a + leap interval (recently a second but for earlier epochs it could be + less) insertion. Note this logic even works when leap intervals + are taken away from UTC (i.e., leap is positive) since in that + case the UTC index always corresponds to the TAI index. */ + *inleap = geMJDtimeTAI_UTC(MJD, &UTC_TO_TAI_lookup_table[index+1]); + return leap; } else { fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad index = %d\n", index); exit(EXIT_FAILURE); @@ -1043,6 +1081,7 @@ MJDtime MJD_value, *MJD=&MJD_value; int forceJulian; double integral_offset1, integral_offset2, integral_scaled_ctime; + int inleap; if(qsasconfig == NULL) { fprintf(stderr, "libqsastime (btimeqsas) ERROR: configqsas must be called first.\n"); @@ -1058,7 +1097,7 @@ forceJulian = 0; if(qsasconfig->ccontrol & 0x2) { - MJD->time_sec += leap_second_TAI(MJD); + MJD->time_sec += leap_second_TAI(MJD, &inleap); } breakDownMJD(year, month, day, hour, min, sec, MJD, forceJulian); @@ -1069,6 +1108,7 @@ MJDtime MJD_value, *MJD=&MJD_value; int forceJulian; double integral_offset1, integral_offset2, integral_scaled_ctime; + int inleap; if(qsasconfig == NULL) { fprintf(stderr, "libqsastime (strfqsas) ERROR: configqsas must be called first.\n"); @@ -1082,7 +1122,11 @@ else forceJulian = 0; - return strfMJD(buf, len, format, MJD, forceJulian); + if(qsasconfig->ccontrol & 0x2) { + MJD->time_sec += leap_second_TAI(MJD, &inleap); + } + + return strfMJD(buf, len, format, MJD, forceJulian, inleap); } /* bhunt_search. Search an ordered table with a binary hunt phase and a This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-03-03 00:38:53
|
Revision: 9661 http://plplot.svn.sourceforge.net/plplot/?rev=9661&view=rev Author: airwin Date: 2009-03-03 00:38:52 +0000 (Tue, 03 Mar 2009) Log Message: ----------- Complete (I hope) all leap interval changes to btimeqsas and strfqsas to optionally calculate and format broken-down UTC values as a function of continous TAI (or any time scale that is a linear transform of TAI). The results appear to be correct for all current tests (which include detailed plots near epochs of a leap interval decrement and both small and large leap second increments) in examples/python xw29.py. The equivalent leap-interval changes to ctimeqsas to optionally calculate continuous time in TAI (or any time scale that is a linear transformation of TAI) as a function of broken-down UTC still need to be implemented. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-03-02 23:22:47 UTC (rev 9660) +++ trunk/lib/qsastime/qsastime.c 2009-03-03 00:38:52 UTC (rev 9661) @@ -322,19 +322,17 @@ } -size_t strfMJD(char * buf, size_t len, const char *format, const MJDtime *MJD, int forceJulian, int if60secformat) +size_t strfMJD(char * buf, size_t len, const char *format, const MJDtime *MJD, int forceJulian, int inleap) { /* Format a text string according to the format string. Uses the same syntax as strftime() but does not use current locale. The null terminator is included in len for safety. */ - /* if if60secformat is true (non-zero) AND the time is somewhere in the - first second of a year, then renormalize to the previous year so that the - date/time comes out beyond the normal last second of the year with the - seconds greater than 60. This form of normalization is used as a flag - that a leap increment (recently always a second, but historically it - was sometimes smaller than that) was in the process of being inserted - for this particular epoch just prior to a discontinuity in TAI-UTC. */ + /* if inleap is true (non-zero) then renormalize so that (int) sec + is 60 to mark results as a flag that a leap increment (recently + always a second, but historically it was sometimes smaller than + that) was in the process of being inserted for this particular + epoch just prior to a positive discontinuity in TAI-UTC. */ int year, month, day, hour, min, ysign, second, d, y; int y1, ifleapyear; @@ -383,20 +381,16 @@ *nMJD = *MJD; nMJD->time_sec += 0.5/shiftPlaces; - normalize_MJD(nMJD); - buf[last] = '\0'; buf[0] = '\0'; /* force overwrite of old buffer since strnctat() used hereafter */ + if(inleap) + nMJD->time_sec -= 1.; + breakDownMJD(&year, &month, &day, &hour, &min, &sec, nMJD, forceJulian); - if(if60secformat && month == 0 && day == 1 && hour == 0 && min == 0 && (int) sec == 0) { - year = year-1; - month = 11; - day = 31; - hour = 23; - min = 59; - sec+= 60.; - } + if(inleap) + sec += 1.; + if(year < 0) { ysign = 1; @@ -1096,11 +1090,17 @@ else forceJulian = 0; - if(qsasconfig->ccontrol & 0x2) { + if(qsasconfig->ccontrol & 0x2) MJD->time_sec += leap_second_TAI(MJD, &inleap); - } + else + inleap = 0; + if(inleap) + MJD->time_sec -= 1.; + breakDownMJD(year, month, day, hour, min, sec, MJD, forceJulian); + if(inleap) + *sec += 1.; } size_t strfqsas(char * buf, size_t len, const char *format, double ctime, const QSASConfig *qsasconfig) @@ -1122,9 +1122,10 @@ else forceJulian = 0; - if(qsasconfig->ccontrol & 0x2) { + if(qsasconfig->ccontrol & 0x2) MJD->time_sec += leap_second_TAI(MJD, &inleap); - } + else + inleap = 0; return strfMJD(buf, len, format, MJD, forceJulian, inleap); } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ai...@us...> - 2009-04-28 18:37:40
|
Revision: 9861 http://plplot.svn.sourceforge.net/plplot/?rev=9861&view=rev Author: airwin Date: 2009-04-28 18:37:39 +0000 (Tue, 28 Apr 2009) Log Message: ----------- Standardize to emacs indentation. Modified Paths: -------------- trunk/lib/qsastime/qsastime.c Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-04-28 17:46:35 UTC (rev 9860) +++ trunk/lib/qsastime/qsastime.c 2009-04-28 18:37:39 UTC (rev 9861) @@ -107,11 +107,11 @@ if(forceJulian) { /* count leap years on proleptic Julian Calendar starting from MJD_0000J */ - leaps = year4 / 4; - if(year%4 == 0) - dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000J; - else - dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000J; + leaps = year4 / 4; + if(year%4 == 0) + dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000J; + else + dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000J; } else { /* count leap years for proleptic Gregorian Calendar. */ /* Algorithm below for 1858-11-17 (0 MJD) gives @@ -128,11 +128,11 @@ calculation with exact representation unless the result is much larger than the integer overflow limit. */ if( (year%4 == 0 && year%100 != 0) || (year%4 == 0 && year%400 == 0) ) - dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000G; - else - dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000G; + dbase_day = year * non_leaps + leaps + MonthStartDOY_L[month] + day + MJD_0000G; + else + dbase_day = year * non_leaps + leaps + MonthStartDOY[month] + day + MJD_0000G; - } + } time_sec = sec + ( (double) min + (double) hour * 60. ) * 60.; @@ -356,25 +356,25 @@ fmtlen = strlen(format); i=0; while(i<fmtlen) - { - char next = format[i]; - if( next == '%') { - /* find seconds format if used */ + char next = format[i]; + if( next == '%') + { + /* find seconds format if used */ + i++; + next = format[i]; + if( isdigit(next) != 0 ) + { + nplaces = strtol(&(format[i]), NULL, 10); + if(nplaces > resolution) resolution = nplaces; + } + else if( next == '.' ) + { + resolution = 9; /* maximum resolution allowed */ + } + } i++; - next = format[i]; - if( isdigit(next) != 0 ) - { - nplaces = strtol(&(format[i]), NULL, 10); - if(nplaces > resolution) resolution = nplaces; - } - else if( next == '.' ) - { - resolution = 9; /* maximum resolution allowed */ - } } - i++; - } /* ensure rounding is done before breakdown */ shiftPlaces = pow(10,(double)resolution); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |