/* New Calendar and Clock 
   Copyright (c) 2008 by Libor Spacek

   Note: time; returns Posix time based on UTC (Universal Temps Coordinat) or 
   TAI (Temps Atomique International) rather than the daylight saving time */
   
using system;   // imports printf, time, ctime, gmtime, gettimeofday, strftime  

namespace date;

private mdayposix jdayposix cycledays secsinday trueyear myyear gregyear 
   lunarmonth fullmoon venussyn venusinf;
   
private infixl 2200 dmod;

public mjday mday2jday mday2pday jday2mday jday2pday pday2jday pday2mday 
   secs2days days2secs days2hours hours2days dhms2secs secs2dhms opdhms 
   hms2secs secs2hms htds2secs secs2htds days2mayan mayan2days addmayan 
   greg2jdays greg2pdays greg2mdays greg2psecs jdays2greg pdays2greg mdays2greg
   gmparse phase completion;
   
public nonfix wday mon day utc year;
   
// constants in whole day units
const mdayposix  = 1856305;// Mayan day for the posix epoch 1 Jan 1970
const jdayposix  = 2440588;// Julian day (since 1 Jan 4713 BC),posix epoch
const cycledays  = 1872000;// end of cycle: total days in 13 Baktuns

// constants in whole seconds units
const secsinday  = 86400;  // number of seconds in a day
const trueyear   = 31556941;// current true year (divisible by 13) 
const myyear     = 31556943;// div by 2277, secsinday compatible = 365.2424 days
const gregyear   = 31556952;// div by 40824, mean gregorian year = 365.2425 days
const lunarmonth = 2551443; // duration of the lunar synodic month
const fullmoon   = 1213810200;//18th June 2008,17:30 full moon in posix seconds
const venussyn   = 50450688; //duration of the Venus synodic cycle
const venusinf   = 1187409600;//18th August 2007,4am Venus inferior conjunction

/* timex macro: prints cputime and realtime in miliseconds 
   (however, cputime on many machines only counts in 10 milisecond units).
   timex can be used safely as a wrapper anywhere in the code */
   
def timex f = f , (void (printf "%s: %dms cpu %8.1fms realtime\n" 
  ((str (quote (f))), (int ((clock-c0)/CLOCKS_PER_SEC*1000+0.5)),
   (1000.0*(gettimeofday-d0))))) when d0 = gettimeofday; c0 = clock end;

/* dmod operator to work on double first argument, 
   so that int, bigint and double times can be conveniently used */
x::double dmod y::int | x::double dmod y::bigint = 
   (x - intx) + (intx mod y) when intx = (int x) end; // mod of a double
x dmod y = x mod y; // polymorphic catchall for ints and bigints

/******************************************************************************/
// either mayan or julian posix epoch (plus posix seconds), gives a double mjday
// to get current pday, use simply: (secs2days time) or (secs2days gettimeofday)
mjday epoch::int secs::int | 
mjday epoch::int secs::bigint |
mjday epoch::int secs::double = epoch+secs/secsinday;// use time or gettimeofday

// all conversions between Julian (j) Mayan (m) and Posix (p), in days
// jday mday pday are numbers of days since their relevant origins (epochs)
jday2mday d::int | jday2mday d::double =  d - jdayposix + mdayposix;
mday2jday d::int | mday2jday d::double =  d - mdayposix + jdayposix; 
jday2pday d::int | jday2pday d::double =  d - jdayposix;
mday2pday d::int | mday2pday d::double =  d - mdayposix;
pday2jday d::int | pday2jday d::double =  d + jdayposix;
pday2mday d::int | pday2mday d::double =  d + mdayposix;

// inner units conversions for convenience and readability   
secs2days s::int | secs2days s::bigint | secs2days s::double = s / secsinday;
days2secs d::int | days2secs d::bigint | days2secs d::double = secsinday * d;
days2hours d::int| days2hours d::bigint| days2hours d::double= 24*d;
hours2days h::int = h / 24;

/* conversions from/to days:hours:minutes:seconds format
   seconds can be int, bigint or double. d,h,m are ints */
dhms2secs (d::int:h::int:m::int:s::int) |
dhms2secs (d::int:h::int:m::int:s::bigint) |
dhms2secs (d::int:h::int:m::int:s::double) = 60*(60*(24*d+h)+m)+s;
secs2dhms secs::int | secs2dhms secs::bigint | secs2dhms secs::double = 
   d:(h mod 24):(m mod 60):(secs dmod 60)
   when m::int = int (secs / 60); h::int = m div 60; d::int = h div 24 end;

/* conversions from/to hours:minutes:seconds format for displaying time of day.
   hours may be more than 24 but use d:h:m:s above for longer periods of time */
hms2secs (h::int:m::int:s::int) |
hms2secs (h::int:m::int:s::bigint) |
hms2secs (h::int:m::int:s::double) = 60*(60*h+m)+s;
secs2hms secs::int | secs2hms secs::bigint | secs2hms secs::double = 
   h:(m mod 60):(secs dmod 60)
   when m::int = int (secs / 60); h::int = m div 60 end;

// any binary operator op applied to two args in (days:hours:minutes:seconds)
opdhms op (d1::int:h1::int:m1::int:s1)(d2::int:h2::int:m2::int:s2) =
   secs2dhms (op (dhms2secs (d1:h1:m1:s1)) (dhms2secs (d2:h2:m2:s2)));

/* New Time Format: hours:3mins:10secs:secs = hours:tres:dicis:secs = h:t:d:s
   the normal seconds are now just a single digit 0-9
   dicis:secs are easy to read: 6:0 means 60 seconds, 12:5 125 seconds etc.
   tres - multiply by three to get traditional babylonian minutes.
   hours remain as usual (24 hour clock) */
htds2secs (h::int:t::int:d::int:s::int)| 
htds2secs (h::int:t::int:d::int:s::bigint)| 
htds2secs (h::int:t::int:d::int:s::double) = 10*(18*(20*h+t)+d)+s;
secs2htds secs::int | secs2htds secs::bigint | secs2htds secs::double =
   h:(t mod 20):(d mod 18):(secs dmod 10)
   when d::int = int (secs / 10); t::int = d div 18; h::int = t div 20 end;

// Mayan 'long count' calendar presentation format   
days2mayan d::int = baktun:(katun mod 20):(tun mod 20):(vinal mod 18):(d mod 20)
   when 
      vinal=d div 20; tun=vinal div 18; katun=tun div 20; baktun=katun div 20
   end;
mayan2days (baktun::int:katun::int:tun::int:vinal::int:kin::int) = 
   20*(18*(20*(20*baktun+katun)+tun)+vinal)+kin;

/* Calculations in Mayan long count format, 
	e.g. addmayan day1 day2 probably not needed, as it is the same as: days2mayan day1+day2; */
addmayan (baktun1::int:katun1::int:tun1::int:vinal1::int:kin1::int)
   	   (baktun2::int:katun2::int:tun2::int:vinal2::int:kin2::int) =
   baktun:(katun mod 20):(tun mod 20):(vinal mod 18):(kin mod 20) 
   when 
   kin = kin1+kin2; vinal = vinal1+vinal2+(kin div 20); 
   tun = tun1+tun2+(vinal div 18); katun = katun1+katun2+(tun div 20);
   baktun = baktun1+baktun2+(katun div 20) 
   end;
 
/* Gregorian calendar presentation format: (d,m,y) 
   from Dr Albert Graef's Q lib and others before.
   unlike the Mayan long count, these dates are historically correct only after
   the introduction of Gregorian calendar in 1582 (much later in non-catholic
   countries). Ten days had been deleted by pope Gregory. However, due to
   ignoring 'pagan' advice, the corrected drift now builds up over the periods 
   of 4,100,200,400 years. This buildup is currently 2.15 days between 1900 and
   2100! On top of that, an uncorrected drift still remains, estimated as the
   minimum of 8 days by the year 12000. These reasons make the Gregorian 
   calendar dates useless for astronomical purposes. 
   Julian days (Salve Julius Caesar!) are still used. Conversions to either 
   Julian days or Mayan days are necessary for any astronomical calculations. 
   
   The Gregorian dates must be represented here as day,month,year tuples */

greg2jdays (d::int, m::int, y::int) = 
   d+(153*m+2) div 5+365*y+y div 4-y div 100+y div 400-32045
   when a = (14-m) div 12; y = y+4800-a; m = m+12*a-3 end;

greg2pdays date = jday2pday (greg2jdays date);
greg2mdays date = jday2mday (greg2jdays date);

// gregorian date time -> psecs
greg2psecs g hms = (days2secs (greg2pdays g)) + (hms2secs hms); 

jdays2greg jd::int  = 
   (e-(153*m+2) div 5+1, m+3-12*(m div 10), 100*b+d-4800+m div 10)
   when a = jd+32044; b = (4*a+3) div 146097;
      c = a-146097*b div 4; d = (4*c+3) div 1461;
      e = c-1461*d div 4; m = (5*e+2) div 153 
   end;

pdays2greg pd::int = jdays2greg (pday2jday pd);
mdays2greg md::int = jdays2greg (mday2jday md);

//parse gmtime string and extract the components
gmparse wday psecs::int| gmparse wday psecs::bigint= asctime (gmtime psecs)!!(0..2);
gmparse mon psecs::int | gmparse mon psecs::bigint = asctime (gmtime psecs)!!(4..6);
gmparse day psecs::int | gmparse day psecs::bigint = asctime (gmtime psecs)!!(8..9);
gmparse utc psecs::int | gmparse utc psecs::bigint = asctime (gmtime psecs)!!(11..18);
gmparse year psecs::int| gmparse year psecs::bigint= asctime (gmtime psecs)!!(20..24);

/* phase of any cycle of 'length' from 'init' at time 'now' (must be same units)
   this is surprisingly accurate without computing the full orbital elements */
phase init::int length::int now::int | 
phase init::int length::int now::bigint |
phase init::int length::int now::double = ((now-init) dmod length)/length;

// same as above but returns dhms till the completion   
completion init::int length::int now::int |
completion init::int length::int now::bigint |
completion init::int length::int now::double = length-((now-init) dmod length);
/******************************************************************************/
cal = do puts [
    "******************************************************************",
    "*     Calendar and Clock, Copyright (c) 2008 by Libor Spacek     *",
    "******************************************************************",
    ((strftime "%X" (localtime time)) + "\t Local time"),
    ((gmparse utc time) + "\t UTC time"),
    ((strftime "%x" (localtime time)) + "\t Gregorian date today") ] $$
//   printf "%s   \t UTC Time in h:t:d:s\n"  
// (str (secs2htds (int (time mod secsinday)))) $$ // seconds since midnight
   printf "%12.4f   \t Julian day\n" jdaytoday $$
   printf "%12.4f   \t Mayan day\n" daytoday $$
   printf "%11.8f%%\t Percentage completion of the Mayan Cycle (MC)\n" 
   	complete $$
   printf "%8.4f   \t Days till the end of MC on %s\n" 
		(nextcycle,(str (jdays2greg (int (jdaytoday+nextcycle))))) $$
   printf "%6.4f   \t Days till the next full Moon on %s\n" 
		(nextfmoon,(str (jdays2greg (int (jdaytoday+nextfmoon))))) $$
   printf "%7.4f   \t Days to inf. conjunction of Venus on %s\n" 
		(nextvenus,(str (jdays2greg (int (jdaytoday+nextvenus))))) $$ 
	do puts [
   "******************************************************************",
   "* Mayan long count digits are [Baktuns:Katuns:Tuns:Vinals:Kins]  *",
   "* Baktun=20Katuns,Katun=20Tuns,Tun=18Vinals,Vinal=20Kins,Kin=day *",
   "* Mayan 13 Baktuns Cycle (MC) = 13*20*20*18*20 = 1,872,000 days  *", 
   "******************************************************************" ] $$
   printf "%s\t Mayan long count date today\n" longtoday $$
   printf "% s\t Mayan long countdown till the end of MC\n\n" mayanleft 
   when
   // let's do some calendric calculations
   nextfmoon = secs2days (completion fullmoon lunarmonth time); // in seconds
   nextvenus = secs2days (completion venusinf venussyn time);
   jdaytoday = mjday jdayposix time; // double julian day-could use gettimeofday
   daytoday  = mjday mdayposix time; // double mayan day-could use gettimeofday
   longtoday = str (days2mayan (int daytoday));
   nextcycle = completion 0 cycledays daytoday; // now in days
   mayanleft = str (days2mayan (succ (int nextcycle)));
   complete = 100.0*(phase 0 cycledays daytoday)
   end;
   
   cal;
   
