#include "map.h"

#define BIG 1.e15
#define HFUZZ .0001

static float hcut[3] ;
static float kr[3] = { .5, -1., .5 };
static float ki[3] = { -1., 0., 1. }; 	/*to multiply by sqrt(3)/2*/
static float cr[3];
static float ci[3];
static struct place hem;
static struct coord twist;
static double  rootroot3, hkc;
static float w2;
static double rootk;

static void
reflect(int i, float wr, float wi, float *x, float *y)
{
	float pr,pi,l;
	pr = cr[i]-wr;
	pi = ci[i]-wi;
	l = 2*(kr[i]*pr + ki[i]*pi);
	*x = wr + l*kr[i];
	*y = wi + l*ki[i];
}

static int
Xhex(struct place *place, float *x, float *y)
{
	int ns;
	register i;
	float zr,zi;
	double sr,si,tr,ti,ur,ui,vr,vi,yr,yi;
	struct place p;
	copyplace(place,&p);
	ns = place->nlat.l >= 0;
	if(!ns) {
		p.nlat.l = -p.nlat.l;
		p.nlat.s = -p.nlat.s;
	}
	if(p.nlat.l<HFUZZ) {
		for(i=0;i<3;i++)
			if(fabs(reduce(p.wlon.l-hcut[i]))<HFUZZ) {
				if(i==2) {
					*x = 2*cr[0] - cr[1];
					*y = 0;
				} else {
					*x = cr[1];
					*y = 2*ci[2*i];
				}
				return(1);
			}
		p.nlat.l = HFUZZ;
		sincos(&p.nlat);
	}
	norm(&p,&hem,&twist);
	Xstereographic(&p,&zr,&zi);
	zr /= 2;
	zi /= 2;
	cdiv(1-zr,-zi,1+zr,zi,&sr,&si);
	csq(sr,si,&tr,&ti);
	ccubrt(1+3*tr,3*ti,&ur,&ui);
	csqrt(ur-1,ui,&vr,&vi);
	cdiv(rootroot3+vr,vi,rootroot3-vr,-vi,&yr,&yi);
	yr /= rootk;
	yi /= rootk;
	elco2(fabs(yr),yi,hkc,1.,1.,x,y);
	if(yr < 0)
		*x = w2 - *x;
	if(!ns) reflect(hcut[0]>place->wlon.l?0:
			hcut[1]>=place->wlon.l?1:
			2,*x,*y,x,y);
	return(1);
}

proj
hex(void)
{
	register i;
	float t;
	double root3;
	float c,d;
	struct place p;
	hcut[2] = PI;
	hcut[1] = hcut[2]/3;
	hcut[0] = -hcut[1];
	root3 = sqrt(3.);
	rootroot3 = sqrt(root3);
	t = 15 -8*root3;
	hkc = t*(1-sqrt(1-1/(t*t)));
	elco2(BIG,0.,hkc,1.,1.,&w2,&t);
	w2 *= 2;
	rootk = sqrt(hkc);
	latlon(90.,90.,&hem);
	latlon(90.,0.,&p);
	Xhex(&p,&c,&t);
	latlon(0.,0.,&p);
	Xhex(&p,&d,&t);
	for(i=0;i<3;i++) {
		ki[i] *= root3/2;
		cr[i] = c + (c-d)*kr[i];
		ci[i] = (c-d)*ki[i];
	}
	deg2rad(0.,&twist);
	return(Xhex);
}

int
hexcut(struct place *g, struct place *og, float *cutlon)
{
	int t,i;
	if(g->nlat.l>=-HFUZZ&&og->nlat.l>=-HFUZZ)
		return(1);
	for(i=0;i<3;i++) {
		t = ckcut(g,og,*cutlon=hcut[i]);
		if(t!=1) return(t);
	}
	return(1);
}