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