From: <airwin@us...>  20090208 00:13:56

Revision: 9471 http://plplot.svn.sourceforge.net/plplot/?rev=9471&view=rev Author: airwin Date: 20090208 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 brokendown 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 20090207 17:26:32 UTC (rev 9470) +++ trunk/lib/qsastime/qsastime.c 20090208 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 18581117 (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. 