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