diff options
Diffstat (limited to 'src/cmd/map/libmap/guyou.c')
-rw-r--r-- | src/cmd/map/libmap/guyou.c | 101 |
1 files changed, 101 insertions, 0 deletions
diff --git a/src/cmd/map/libmap/guyou.c b/src/cmd/map/libmap/guyou.c new file mode 100644 index 00000000..b37736f2 --- /dev/null +++ b/src/cmd/map/libmap/guyou.c @@ -0,0 +1,101 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static struct place gywhem, gyehem; +static struct coord gytwist; +static double gyconst, gykc, gyside; + + +static void +dosquare(double z1, double z2, double *x, double *y) +{ + double w1,w2; + w1 = z1 -1; + if(fabs(w1*w1+z2*z2)>.000001) { + cdiv(z1+1,z2,w1,z2,&w1,&w2); + w1 *= gyconst; + w2 *= gyconst; + if(w1<0) + w1 = 0; + elco2(w1,w2,gykc,1.,1.,x,y); + } else { + *x = gyside; + *y = 0; + } +} + +int +Xguyou(struct place *place, double *x, double *y) +{ + int ew; /*which hemisphere*/ + double z1,z2; + struct place pl; + ew = place->wlon.l<0; + copyplace(place,&pl); + norm(&pl,ew?&gyehem:&gywhem,&gytwist); + Xstereographic(&pl,&z1,&z2); + dosquare(z1/2,z2/2,x,y); + if(!ew) + *x -= gyside; + return(1); +} + +proj +guyou(void) +{ + double junk; + gykc = 1/(3+2*sqrt(2.)); + gyconst = -(1+sqrt(2.)); + elco2(-gyconst,0.,gykc,1.,1.,&gyside,&junk); + gyside *= 2; + latlon(0.,90.,&gywhem); + latlon(0.,-90.,&gyehem); + deg2rad(0.,&gytwist); + return(Xguyou); +} + +int +guycut(struct place *g, struct place *og, double *cutlon) +{ + int c; + c = picut(g,og,cutlon); + if(c!=1) + return(c); + *cutlon = 0.; + if(g->nlat.c<.7071||og->nlat.c<.7071) + return(ckcut(g,og,0.)); + return(1); +} + +static int +Xsquare(struct place *place, double *x, double *y) +{ + double z1,z2; + double r, theta; + struct place p; + copyplace(place,&p); + if(place->nlat.l<0) { + p.nlat.l = -p.nlat.l; + p.nlat.s = -p.nlat.s; + } + if(p.nlat.l<FUZZ && fabs(p.wlon.l)>PI-FUZZ){ + *y = -gyside/2; + *x = p.wlon.l>0?0:gyside; + return(1); + } + Xstereographic(&p,&z1,&z2); + r = sqrt(sqrt(hypot(z1,z2)/2)); + theta = atan2(z1,-z2)/4; + dosquare(r*sin(theta),-r*cos(theta),x,y); + if(place->nlat.l<0) + *y = -gyside - *y; + return(1); +} + +proj +square(void) +{ + guyou(); + return(Xsquare); +} |