1*59cc4ca5SDavid du Colombier #include <u.h>
2*59cc4ca5SDavid du Colombier #include <libc.h>
33e12c5d1SDavid du Colombier #include "map.h"
43e12c5d1SDavid du Colombier
53e12c5d1SDavid du Colombier static struct place gywhem, gyehem;
63e12c5d1SDavid du Colombier static struct coord gytwist;
7219b2ee8SDavid du Colombier static double gyconst, gykc, gyside;
83e12c5d1SDavid du Colombier
93e12c5d1SDavid du Colombier
103e12c5d1SDavid du Colombier static void
dosquare(double z1,double z2,double * x,double * y)11219b2ee8SDavid du Colombier dosquare(double z1, double z2, double *x, double *y)
123e12c5d1SDavid du Colombier {
133e12c5d1SDavid du Colombier double w1,w2;
143e12c5d1SDavid du Colombier w1 = z1 -1;
153e12c5d1SDavid du Colombier if(fabs(w1*w1+z2*z2)>.000001) {
163e12c5d1SDavid du Colombier cdiv(z1+1,z2,w1,z2,&w1,&w2);
173e12c5d1SDavid du Colombier w1 *= gyconst;
183e12c5d1SDavid du Colombier w2 *= gyconst;
193e12c5d1SDavid du Colombier if(w1<0)
203e12c5d1SDavid du Colombier w1 = 0;
213e12c5d1SDavid du Colombier elco2(w1,w2,gykc,1.,1.,x,y);
223e12c5d1SDavid du Colombier } else {
233e12c5d1SDavid du Colombier *x = gyside;
243e12c5d1SDavid du Colombier *y = 0;
253e12c5d1SDavid du Colombier }
263e12c5d1SDavid du Colombier }
273e12c5d1SDavid du Colombier
283e12c5d1SDavid du Colombier int
Xguyou(struct place * place,double * x,double * y)29219b2ee8SDavid du Colombier Xguyou(struct place *place, double *x, double *y)
303e12c5d1SDavid du Colombier {
313e12c5d1SDavid du Colombier int ew; /*which hemisphere*/
32219b2ee8SDavid du Colombier double z1,z2;
333e12c5d1SDavid du Colombier struct place pl;
343e12c5d1SDavid du Colombier ew = place->wlon.l<0;
353e12c5d1SDavid du Colombier copyplace(place,&pl);
363e12c5d1SDavid du Colombier norm(&pl,ew?&gyehem:&gywhem,&gytwist);
373e12c5d1SDavid du Colombier Xstereographic(&pl,&z1,&z2);
383e12c5d1SDavid du Colombier dosquare(z1/2,z2/2,x,y);
393e12c5d1SDavid du Colombier if(!ew)
403e12c5d1SDavid du Colombier *x -= gyside;
413e12c5d1SDavid du Colombier return(1);
423e12c5d1SDavid du Colombier }
433e12c5d1SDavid du Colombier
443e12c5d1SDavid du Colombier proj
guyou(void)453e12c5d1SDavid du Colombier guyou(void)
463e12c5d1SDavid du Colombier {
47219b2ee8SDavid du Colombier double junk;
483e12c5d1SDavid du Colombier gykc = 1/(3+2*sqrt(2.));
493e12c5d1SDavid du Colombier gyconst = -(1+sqrt(2.));
503e12c5d1SDavid du Colombier elco2(-gyconst,0.,gykc,1.,1.,&gyside,&junk);
513e12c5d1SDavid du Colombier gyside *= 2;
523e12c5d1SDavid du Colombier latlon(0.,90.,&gywhem);
533e12c5d1SDavid du Colombier latlon(0.,-90.,&gyehem);
543e12c5d1SDavid du Colombier deg2rad(0.,&gytwist);
553e12c5d1SDavid du Colombier return(Xguyou);
563e12c5d1SDavid du Colombier }
573e12c5d1SDavid du Colombier
583e12c5d1SDavid du Colombier int
guycut(struct place * g,struct place * og,double * cutlon)59219b2ee8SDavid du Colombier guycut(struct place *g, struct place *og, double *cutlon)
603e12c5d1SDavid du Colombier {
613e12c5d1SDavid du Colombier int c;
623e12c5d1SDavid du Colombier c = picut(g,og,cutlon);
633e12c5d1SDavid du Colombier if(c!=1)
643e12c5d1SDavid du Colombier return(c);
653e12c5d1SDavid du Colombier *cutlon = 0.;
663e12c5d1SDavid du Colombier if(g->nlat.c<.7071||og->nlat.c<.7071)
673e12c5d1SDavid du Colombier return(ckcut(g,og,0.));
683e12c5d1SDavid du Colombier return(1);
693e12c5d1SDavid du Colombier }
703e12c5d1SDavid du Colombier
713e12c5d1SDavid du Colombier static int
Xsquare(struct place * place,double * x,double * y)72219b2ee8SDavid du Colombier Xsquare(struct place *place, double *x, double *y)
733e12c5d1SDavid du Colombier {
74219b2ee8SDavid du Colombier double z1,z2;
75219b2ee8SDavid du Colombier double r, theta;
763e12c5d1SDavid du Colombier struct place p;
773e12c5d1SDavid du Colombier copyplace(place,&p);
783e12c5d1SDavid du Colombier if(place->nlat.l<0) {
793e12c5d1SDavid du Colombier p.nlat.l = -p.nlat.l;
803e12c5d1SDavid du Colombier p.nlat.s = -p.nlat.s;
813e12c5d1SDavid du Colombier }
823e12c5d1SDavid du Colombier if(p.nlat.l<FUZZ && fabs(p.wlon.l)>PI-FUZZ){
833e12c5d1SDavid du Colombier *y = -gyside/2;
843e12c5d1SDavid du Colombier *x = p.wlon.l>0?0:gyside;
853e12c5d1SDavid du Colombier return(1);
863e12c5d1SDavid du Colombier }
873e12c5d1SDavid du Colombier Xstereographic(&p,&z1,&z2);
883e12c5d1SDavid du Colombier r = sqrt(sqrt(hypot(z1,z2)/2));
893e12c5d1SDavid du Colombier theta = atan2(z1,-z2)/4;
903e12c5d1SDavid du Colombier dosquare(r*sin(theta),-r*cos(theta),x,y);
913e12c5d1SDavid du Colombier if(place->nlat.l<0)
923e12c5d1SDavid du Colombier *y = -gyside - *y;
933e12c5d1SDavid du Colombier return(1);
943e12c5d1SDavid du Colombier }
953e12c5d1SDavid du Colombier
963e12c5d1SDavid du Colombier proj
square(void)973e12c5d1SDavid du Colombier square(void)
983e12c5d1SDavid du Colombier {
993e12c5d1SDavid du Colombier guyou();
1003e12c5d1SDavid du Colombier return(Xsquare);
1013e12c5d1SDavid du Colombier }
102