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