hex.c (2273B)
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 }