1 #include "map.h" 2 3 static struct place gywhem, gyehem; 4 static struct coord gytwist; 5 static double gyconst, gykc, gyside; 6 7 8 static void 9 dosquare(double z1, double z2, double *x, double *y) 10 { 11 double w1,w2; 12 w1 = z1 -1; 13 if(fabs(w1*w1+z2*z2)>.000001) { 14 cdiv(z1+1,z2,w1,z2,&w1,&w2); 15 w1 *= gyconst; 16 w2 *= gyconst; 17 if(w1<0) 18 w1 = 0; 19 elco2(w1,w2,gykc,1.,1.,x,y); 20 } else { 21 *x = gyside; 22 *y = 0; 23 } 24 } 25 26 int 27 Xguyou(struct place *place, double *x, double *y) 28 { 29 int ew; /*which hemisphere*/ 30 double z1,z2; 31 struct place pl; 32 ew = place->wlon.l<0; 33 copyplace(place,&pl); 34 norm(&pl,ew?&gyehem:&gywhem,&gytwist); 35 Xstereographic(&pl,&z1,&z2); 36 dosquare(z1/2,z2/2,x,y); 37 if(!ew) 38 *x -= gyside; 39 return(1); 40 } 41 42 proj 43 guyou(void) 44 { 45 double junk; 46 gykc = 1/(3+2*sqrt(2.)); 47 gyconst = -(1+sqrt(2.)); 48 elco2(-gyconst,0.,gykc,1.,1.,&gyside,&junk); 49 gyside *= 2; 50 latlon(0.,90.,&gywhem); 51 latlon(0.,-90.,&gyehem); 52 deg2rad(0.,&gytwist); 53 return(Xguyou); 54 } 55 56 int 57 guycut(struct place *g, struct place *og, double *cutlon) 58 { 59 int c; 60 c = picut(g,og,cutlon); 61 if(c!=1) 62 return(c); 63 *cutlon = 0.; 64 if(g->nlat.c<.7071||og->nlat.c<.7071) 65 return(ckcut(g,og,0.)); 66 return(1); 67 } 68 69 static int 70 Xsquare(struct place *place, double *x, double *y) 71 { 72 double z1,z2; 73 double r, theta; 74 struct place p; 75 copyplace(place,&p); 76 if(place->nlat.l<0) { 77 p.nlat.l = -p.nlat.l; 78 p.nlat.s = -p.nlat.s; 79 } 80 if(p.nlat.l<FUZZ && fabs(p.wlon.l)>PI-FUZZ){ 81 *y = -gyside/2; 82 *x = p.wlon.l>0?0:gyside; 83 return(1); 84 } 85 Xstereographic(&p,&z1,&z2); 86 r = sqrt(sqrt(hypot(z1,z2)/2)); 87 theta = atan2(z1,-z2)/4; 88 dosquare(r*sin(theta),-r*cos(theta),x,y); 89 if(place->nlat.l<0) 90 *y = -gyside - *y; 91 return(1); 92 } 93 94 proj 95 square(void) 96 { 97 guyou(); 98 return(Xsquare); 99 } 100