1 #include <u.h>
2 #include <libc.h>
3 #include "map.h"
4 
5 #define BIG 1.e15
6 #define HFUZZ .0001
7 
8 static double hcut[3] ;
9 static double kr[3] = { .5, -1., .5 };
10 static double ki[3] = { -1., 0., 1. }; 	/*to multiply by sqrt(3)/2*/
11 static double cr[3];
12 static double ci[3];
13 static struct place hem;
14 static struct coord twist;
15 static double  rootroot3, hkc;
16 static double w2;
17 static double rootk;
18 
19 static void
20 reflect(int i, double wr, double wi, double *x, double *y)
21 {
22 	double pr,pi,l;
23 	pr = cr[i]-wr;
24 	pi = ci[i]-wi;
25 	l = 2*(kr[i]*pr + ki[i]*pi);
26 	*x = wr + l*kr[i];
27 	*y = wi + l*ki[i];
28 }
29 
30 static int
31 Xhex(struct place *place, double *x, double *y)
32 {
33 	int ns;
34 	int i;
35 	double zr,zi;
36 	double sr,si,tr,ti,ur,ui,vr,vi,yr,yi;
37 	struct place p;
38 	copyplace(place,&p);
39 	ns = place->nlat.l >= 0;
40 	if(!ns) {
41 		p.nlat.l = -p.nlat.l;
42 		p.nlat.s = -p.nlat.s;
43 	}
44 	if(p.nlat.l<HFUZZ) {
45 		for(i=0;i<3;i++)
46 			if(fabs(reduce(p.wlon.l-hcut[i]))<HFUZZ) {
47 				if(i==2) {
48 					*x = 2*cr[0] - cr[1];
49 					*y = 0;
50 				} else {
51 					*x = cr[1];
52 					*y = 2*ci[2*i];
53 				}
54 				return(1);
55 			}
56 		p.nlat.l = HFUZZ;
57 		sincos(&p.nlat);
58 	}
59 	norm(&p,&hem,&twist);
60 	Xstereographic(&p,&zr,&zi);
61 	zr /= 2;
62 	zi /= 2;
63 	cdiv(1-zr,-zi,1+zr,zi,&sr,&si);
64 	csq(sr,si,&tr,&ti);
65 	ccubrt(1+3*tr,3*ti,&ur,&ui);
66 	csqrt(ur-1,ui,&vr,&vi);
67 	cdiv(rootroot3+vr,vi,rootroot3-vr,-vi,&yr,&yi);
68 	yr /= rootk;
69 	yi /= rootk;
70 	elco2(fabs(yr),yi,hkc,1.,1.,x,y);
71 	if(yr < 0)
72 		*x = w2 - *x;
73 	if(!ns) reflect(hcut[0]>place->wlon.l?0:
74 			hcut[1]>=place->wlon.l?1:
75 			2,*x,*y,x,y);
76 	return(1);
77 }
78 
79 proj
80 hex(void)
81 {
82 	int i;
83 	double t;
84 	double root3;
85 	double c,d;
86 	struct place p;
87 	hcut[2] = PI;
88 	hcut[1] = hcut[2]/3;
89 	hcut[0] = -hcut[1];
90 	root3 = sqrt(3.);
91 	rootroot3 = sqrt(root3);
92 	t = 15 -8*root3;
93 	hkc = t*(1-sqrt(1-1/(t*t)));
94 	elco2(BIG,0.,hkc,1.,1.,&w2,&t);
95 	w2 *= 2;
96 	rootk = sqrt(hkc);
97 	latlon(90.,90.,&hem);
98 	latlon(90.,0.,&p);
99 	Xhex(&p,&c,&t);
100 	latlon(0.,0.,&p);
101 	Xhex(&p,&d,&t);
102 	for(i=0;i<3;i++) {
103 		ki[i] *= root3/2;
104 		cr[i] = c + (c-d)*kr[i];
105 		ci[i] = (c-d)*ki[i];
106 	}
107 	deg2rad(0.,&twist);
108 	return(Xhex);
109 }
110 
111 int
112 hexcut(struct place *g, struct place *og, double *cutlon)
113 {
114 	int t,i;
115 	if(g->nlat.l>=-HFUZZ&&og->nlat.l>=-HFUZZ)
116 		return(1);
117 	for(i=0;i<3;i++) {
118 		t = ckcut(g,og,*cutlon=hcut[i]);
119 		if(t!=1) return(t);
120 	}
121 	return(1);
122 }