#include<u.h>
#include<libc.h>

double longitude = 2.0 + 20.0 / 60.0;
const double eccentricity = 0.0167;
const double obliquity = 23.44;
const double deg2rad = 2 * PI / 360;

/* Adjusts a time structure by delta min */
void
correct(Tm *t, double delta)
{
	int hours, min, sec;
	
	hours = delta / 60;
	min = delta - hours * 60;
	sec = (delta - (int)delta) * 60;
	t->sec += sec;
	if(t->sec >= 60) t->sec -= 60, t->min++;
	if(t->sec < 0) t->sec += 60, t->min--;
	t->min += min;
	if(t->min >= 60) t->min -= 60, t->hour++;
	if(t->min < 0) t->min += 60, t->hour--;
	t->hour += hours;
	if(t->hour >= 24) t->hour -= 24;
	if(t->hour < 0) t->hour += 24;
}

double
dsin(double x)
{
	return sin(x * deg2rad);
}

double
dcos(double x)
{
	return cos(x * deg2rad);
}

double
dtan(double x)
{
	return tan(x * deg2rad);
}

double
datan(double x)
{
	return atan(x) / deg2rad;
}

/* Hex formatter */
#pragma varargck type "H" Tm*
int
hexfmt(Fmt *f)
{
	int hours, min, sec;
	uint time;
	Tm *t;
	
	t = va_arg(f->args, Tm *);
	time = 0.5+(float)0x10000 * ((float)t->hour + (float)t->min / 60.0 + (float)t->sec / (60.0 * 60.0)) / 24.0;
	hours = time / 0x1000;
	min = time % 0x1000 / 0x10;
	sec = time % 0x10;
	return fmtprint(f, "%01X_%02X_%01X", hours, min, sec);
}

/* Duodecimal formatter */
#pragma varargck type "D" Tm*
int
duofmt(Fmt *f)
{
	Tm *t;
	t = va_arg(f->args, Tm*);
	return fmtprint(f, "%02d:%02d:%02d", t->hour, t->min, t->sec);
}

void
usage(void)
{
	fprint(2, "%s: [-l longitude] [-x] [-f]", argv0);
	exits("usage");
}

void
main(int argc, char *argv[])
{
	Tm *t;
	int hex, follow;
	char *fmt;
	
	fmt = "%D\n";
	hex = 0;
	follow = 0;
	
	double w, a, b, c, eot;
	
	ARGBEGIN{
	case 'l':
		longitude = atof(EARGF(usage()));
		break;
	case 'x':
		fmt = "%H\n";
		hex++;
		break;
	case 'f':
		follow++;
		break;
	case 'h':
	default:
		usage();
	}ARGEND;
	
	if(argc != 0)
		usage();
		
	fmtinstall('H', hexfmt);
	fmtinstall('D', duofmt);
	
beginning:
	/* Get local mean time */
	t = gmtime(time(0));
	correct(t, longitude * 4.0);
	
	/* Calculate equation of time */
	w = 360 / 365.24; // Mean angular orbit velocity
	a = w * (t->yday + 10); // angular orbital distance from solstice (10 days befor jan 1)
	b = a + (360 / PI) * eccentricity * dsin(w * (t->yday - 2)); // angular deflection of earth
	c = (a - datan(dtan(b) / dcos(obliquity))) / 180; // Difference between mean and corrected angles
	eot = 720 * (c - floor(c + 0.5));
	
	/* Alternative formulation */
	/* w = 2 * PI * (t->yday - 1) / 365.242;
	eot = 0.258 * cos(w) - 7.416 * sin(w) - 3.648 * cos(2 * w) - 9.228 * sin(2 * w);*/
	
	//print("%f\n", eot); 
	correct(t, -eot);
	print(fmt, t);
		
	if(follow){
		sleep(hex ? (1000 / 0.76) : 1000);
		goto beginning;
	}
	
	exits(0);
}