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

#include "sky.h"

// Munich, Germany
Obs defaultObs = {
	{degms2rad(48, 8, 14.74), degms2rad(11, 34, 31.76), 524.0,},
	NTP_kPa,
	NTP_K,
	0,
};

// Canon EOS 600D with a 24mm lens
Camera defaultCamera = {
	24.0,
	22.3,
	14.9,
	5344,
	3516,
	{0, 0, 0,},
};

double _sdG = 0;
double _cdG = 0;

double gammab_fw(double jd, double dut1, double dt);
double phib_fw(double jd, double dut1, double dt);
double psib_fw(double jd, double dut1, double dt);
double epsA_fw(double jd, double dut1, double dt);
void sXY(double jd, double dut1, double dt, double *s, double *X, double *Y);

Eqloc
gx2eq00(Gxloc gx)
{
	double sb, cb, slCPl, clCPl, sd, saaGcd, caaGcd;
	Eqloc eq;

	if(_sdG==0) _sdG = sin(deg2rad(27.12825));
	if(_cdG==0) _cdG = cos(deg2rad(27.12825));

	sb = sin(gx.b);
	cb = cos(gx.b);
	slCPl = sin(deg2rad(122.93192)-gx.l);
	clCPl = cos(deg2rad(122.93192)-gx.l);

	sd = sb*_sdG + cb*_cdG*clCPl;

	saaGcd = cb*slCPl;
	caaGcd = sb*_cdG - cb*_sdG*clCPl;

	eq.dec = asin(sd);
	eq.ra = atan2(saaGcd, caaGcd) + deg2rad(192.85948);
	while(eq.ra<-PIO2) eq.ra+=PI;
	while(eq.ra>PIO2) eq.ra -= PI;

	return eq;
}

Gxloc
eq2gx00(Eqloc eq)
{
	double sd, cd, saaG, caaG, sb, slCPlcb, clCPlcb;
	Gxloc gx;

	if(_sdG==0) _sdG = sin(deg2rad(27.12825));
	if(_cdG==0) _cdG = cos(deg2rad(27.12825));

	sd = sin(eq.dec);
	cd = cos(eq.dec);
	saaG = sin(eq.ra - deg2rad(192.85948));
	caaG = cos(eq.ra - deg2rad(192.85948));

	sb = sd*_sdG + cd*_cdG*caaG;

	slCPlcb = cd*saaG;
	clCPlcb = sd*_cdG - cd*_sdG*caaG;

	gx.b = asin(sb);
	gx.l = deg2rad(122.93192) - atan2(slCPlcb, clCPlcb);
	while(gx.l<-PIO2) gx.l+=PI;
	while(gx.l>PIO2) gx.l -= PI;

	return gx;
}

Hzloc
eq2hz(Eqloc eq, Gdloc geod, double jdutc, double dut1, double dt)
{
	double gst, /* tU, era, */ h, sph, cph, sd, cd, sh, ch, u, v, w;
	Hzloc hz;

	gst = gast(jdutc, dut1, dt);
	h = gst + geod.lng - eq.ra;
/*
	tU = jd+dut1-JD00;
	era = TWOPI * (tU -(long)tU + 0.7790572732640 + 0.00273781191135448 * tU);
	h = era + geod.lng - eq.ra;
*/
	sph = sin(geod.lat);
	cph = cos(geod.lat);
	sd = sin(eq.dec);
	cd = cos(eq.dec);
	sh = sin(h);
	ch = cos(h);

	u = sd*cph - cd*ch*sph;
	v = -cd*sh;
	w = sd*sph + cd*ch*cph;

	hz.az = atan2(v, u);
	hz.az = hz.az < 0.0 ? hz.az+TWOPI : hz.az;

	hz.z = acos(w);
	//hz.z = PIO2 - atan2(w, sqrt(u*u+v*v));

	return hz;
}

Eqloc
hz2eq(Hzloc hz, Gdloc geod, double jdutc, double dut1, double dt)
{
	double sph, cph, saz, caz, sz, cz, u, v, w, h, gst;
	Eqloc eq;

	sph = sin(geod.lat);
	cph = cos(geod.lat);
	saz = sin(hz.az);
	caz = cos(hz.az);
	sz = sin(hz.z);
	cz = cos(hz.z);

	u = cz*cph - caz*sz*sph;
	v = -saz*sz;
	w = cz*sph + caz*sz*cph;

	h = atan2(v, u);
	gst = gast(jdutc, dut1, dt);
	
	eq.ra = gst + geod.lng - h;
	eq.dec = asin(w);
	//eq.dec = atan2(w, sqrt(u*u+v*v));

	return eq;
}

Planar
hz2pxl(Camera *c, Hzloc hz)
{
	double dz, daz, tz, taz, sr, cr;
	Planar pt, ptr, px;

	dz = hz.z - (c->angle).z;
	daz = hz.az - (c->angle).az;

	tz = tan(dz);
	taz = tan(daz);

	sr = sin((c->angle).roll);
	cr = cos((c->angle).roll);

	pt.x = -c->f*tz;
	pt.y = c->f*taz;
	ptr.x = pt.x*cr-pt.y*sr;
	ptr.y = pt.x*sr+pt.y*cr;

	while(dz < -PI) dz += TWOPI;
	while(dz > PI) dz -= TWOPI;
	while(daz < -PI) daz += TWOPI;
	while(daz > PI) daz -= TWOPI;

	/* behind the camera */
	if(dz < -PIO2 || dz > PIO2 || daz < -PIO2 || daz > PIO2){
		px.x = NaN();
		px.y = NaN();
		return px;
	}

	px.x = ptr.x*c->wpx/c->w;
	px.y = ptr.y*c->hpx/c->h;
	return px;
}

Hzloc
pxl2hz(Camera *c, Planar px)
{
	double sr, cr, tz, taz;
	Planar pt, ptr;
	Hzloc hz;

	sr = sin(c->angle.roll);
	cr = cos(c->angle.roll);

	ptr.x = px.x*c->w/c->wpx;
	ptr.y = px.y*c->h/c->hpx;
	pt.x = ptr.x*cr+ptr.y*sr;
	pt.y = -ptr.x*sr+ptr.y*cr;

	tz = -pt.x / c->f;
	taz = -pt.y / c->f;

	hz.z = c->angle.z + atan(tz);
	hz.az = c->angle.az + atan(taz);

	return hz;
}

double
gast(double jd, double dut1, double dt)
{
	double tU, t, era, gmst;

	tU = jd+dut1-JD00;
	t = (tU+dt)/JDX100;

	/* earth rotation angle is used for the CIO-based hour angle */
	era = TWOPI * (tU -(long)tU + 0.7790572732640 + 0.00273781191135448 * tU);
	
	/* from Fukushima 2003 */
	gmst = era +
		(0.012911569 + t*4612.160517397) * DSEC +
		t*t*(1.391542507 + t*(-0.000124849 + t*(-0.000004991 + t*(-0.000000479)))) * DSEC;

//	/* GMST approx from Capitaine et al. 2005 */
//	gmst = era + DSEC * (((((
//		-0.0000000368)*t +
//		-0.000029956)*t +
//		-0.00000044)*t +
//		1.3915817)*t +
//		4612.156534)*t +
//		0.014506;

	return gmst+eqeqx(jd, dut1, dt);
}

double
eqeqx(double jd, double dut1, double dt)
{
	double tU, t, om, lm, ls, gs, gls, dpsi, epsbar;

	tU = jd+dut1-JD00;
	t = (tU+dt)/JDX100;

	/* longitude of the scending node of moon's orbit */
	om = (125.04452 + t*(-134.136261 + t*(0.0020708 + t/450000.0))) * DEG;

	/* mean longitude of the moon */
	lm = (218.31654591 + t*(481267.88134236 + t*(-0.00163 + t*(1./538841.0 - t/65194000.0)))) * DEG;

	/* mean longitude of the sun */
	ls = (280.46645 + t*(36000.76983 + t*0.0003032)) * DEG;

	/* mean anomaly of the sun */
	gs = (357.52910 + t*(35999.05030 + t*( -0.0001559 + t*(-0.00000048)))) * DEG;

	/* geometric longitude of the sun */
	gls = ls + ( (1.9146000 + t*(-0.004817 + t*(- 0.000014)))*sin(gs) + (0.019993 - 0.000101*t)*sin(2.0*gs) + 0.000290*sin(3.0*gs) ) * DEG;
	
	dpsi =  ( -17.20*sin(om) - 1.32*sin(gls+gls) - 0.23*sin(lm+lm) + 0.21*sin(om+om) ) * DSEC;
	epsbar = (84381.448 + t*(-46.8150 + t*(-0.00059 + t*0.001813))) * DSEC;

	return dpsi*cos(epsbar) + (0.00264*sin(om) + 0.000063*sin(om+om)) * DSEC;	
}

double
saemundsson(Obs *o, double z)
{
	return -h2rad(1.02/60.0) *
		o->P/101.0 *
		283.0/o->T /
		tan(PIO2-z+deg2rad(10.3)/
			(PIO2-z+deg2rad(5.11)));
}


/* functions below this line are not used at the moment but kept for future use */

double
gammab_fw(double jd, double dut1, double dt)
{
	double t;

	t = (jd+dut1+dt-JD00) / JDX100;

	/* from Wallace & Capitaine, 2006 (Fukushima-Williams equations) */
	return ((((((
		0.0000000260)*t - 
		0.000002788)*t -
		0.00031238)*t +
		0.4932044)*t +
		10.556378)*t -
		0.052928) * DSEC;
}

double
phib_fw(double jd, double dut1, double dt)
{
	double t;

	t = (jd+dut1+dt-JD00) / JDX100;

	/* from Wallace & Capitaine, 2006 (Fukushima-Williams equations) */
	return ((((((
		-0.0000000176)*t -
		0.000000440)*t +
		0.00053289)*t +
		0.0511268)*t -
		46.811016)*t +
		84381.412819) * DSEC;
}

double
psib_fw(double jd, double dut1, double dt)
{
	double t;

	t = (jd+dut1+dt-JD00) / JDX100;

	/* from Wallace & Capitaine, 2006 (Fukushima-Williams equations) */
	return ((((((
		-0.0000000148)*t +
		0.000026452)*t -
		0.00018522)*t -
		1.5584175)*t +
		5038.481484)*t -
		0.041775) * DSEC;
}

double
epsA_fw(double jd, double dut1, double dt)
{
	double t;

	t = (jd+dut1+dt-JD00) / JDX100;

	/* from Wallace & Capitaine, 2006 (Fukushima-Williams equations) */
	return ((((((
		-0.0000000434)*t -
		0.000000576)*t +
		0.00200340)*t - 
		0.0001831)*t -
		46.836769)*t +
		84381.406) * DSEC;
}

void
sXY(double jd, double dut1, double dt, double *s, double *X, double *Y)
{
	double tau, t, sXY2, Omega;

	tau = jd+dut1+dt-JD00;
	t = tau / JDX100;

	sXY2 = (((((15.62*t + 27.98)*t - 72574.11)*t - 122.68)*t + 3808.65)*t + 94.0);
	Omega = 2.182 - 9.242e-4*t;

	*X = 2.6603e-7*t - 33.2e-6*sin(Omega);
	*Y = -8.14e-14*t*t+44.6e-6*cos(Omega);
	
	*s = sXY2 - (*X)*(*Y)/2;

	return;	
}