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 #define BIG 1.e15
63e12c5d1SDavid du Colombier #define HFUZZ .0001
73e12c5d1SDavid du Colombier
8219b2ee8SDavid du Colombier static double hcut[3] ;
9219b2ee8SDavid du Colombier static double kr[3] = { .5, -1., .5 };
10219b2ee8SDavid du Colombier static double ki[3] = { -1., 0., 1. }; /*to multiply by sqrt(3)/2*/
11219b2ee8SDavid du Colombier static double cr[3];
12219b2ee8SDavid du Colombier static double ci[3];
133e12c5d1SDavid du Colombier static struct place hem;
143e12c5d1SDavid du Colombier static struct coord twist;
153e12c5d1SDavid du Colombier static double rootroot3, hkc;
16219b2ee8SDavid du Colombier static double w2;
173e12c5d1SDavid du Colombier static double rootk;
183e12c5d1SDavid du Colombier
193e12c5d1SDavid du Colombier static void
reflect(int i,double wr,double wi,double * x,double * y)20219b2ee8SDavid du Colombier reflect(int i, double wr, double wi, double *x, double *y)
213e12c5d1SDavid du Colombier {
22219b2ee8SDavid du Colombier double pr,pi,l;
233e12c5d1SDavid du Colombier pr = cr[i]-wr;
243e12c5d1SDavid du Colombier pi = ci[i]-wi;
253e12c5d1SDavid du Colombier l = 2*(kr[i]*pr + ki[i]*pi);
263e12c5d1SDavid du Colombier *x = wr + l*kr[i];
273e12c5d1SDavid du Colombier *y = wi + l*ki[i];
283e12c5d1SDavid du Colombier }
293e12c5d1SDavid du Colombier
303e12c5d1SDavid du Colombier static int
Xhex(struct place * place,double * x,double * y)31219b2ee8SDavid du Colombier Xhex(struct place *place, double *x, double *y)
323e12c5d1SDavid du Colombier {
333e12c5d1SDavid du Colombier int ns;
343e12c5d1SDavid du Colombier register i;
35219b2ee8SDavid du Colombier double zr,zi;
363e12c5d1SDavid du Colombier double sr,si,tr,ti,ur,ui,vr,vi,yr,yi;
373e12c5d1SDavid du Colombier struct place p;
383e12c5d1SDavid du Colombier copyplace(place,&p);
393e12c5d1SDavid du Colombier ns = place->nlat.l >= 0;
403e12c5d1SDavid du Colombier if(!ns) {
413e12c5d1SDavid du Colombier p.nlat.l = -p.nlat.l;
423e12c5d1SDavid du Colombier p.nlat.s = -p.nlat.s;
433e12c5d1SDavid du Colombier }
443e12c5d1SDavid du Colombier if(p.nlat.l<HFUZZ) {
453e12c5d1SDavid du Colombier for(i=0;i<3;i++)
463e12c5d1SDavid du Colombier if(fabs(reduce(p.wlon.l-hcut[i]))<HFUZZ) {
473e12c5d1SDavid du Colombier if(i==2) {
483e12c5d1SDavid du Colombier *x = 2*cr[0] - cr[1];
493e12c5d1SDavid du Colombier *y = 0;
503e12c5d1SDavid du Colombier } else {
513e12c5d1SDavid du Colombier *x = cr[1];
523e12c5d1SDavid du Colombier *y = 2*ci[2*i];
533e12c5d1SDavid du Colombier }
543e12c5d1SDavid du Colombier return(1);
553e12c5d1SDavid du Colombier }
563e12c5d1SDavid du Colombier p.nlat.l = HFUZZ;
573e12c5d1SDavid du Colombier sincos(&p.nlat);
583e12c5d1SDavid du Colombier }
593e12c5d1SDavid du Colombier norm(&p,&hem,&twist);
603e12c5d1SDavid du Colombier Xstereographic(&p,&zr,&zi);
613e12c5d1SDavid du Colombier zr /= 2;
623e12c5d1SDavid du Colombier zi /= 2;
633e12c5d1SDavid du Colombier cdiv(1-zr,-zi,1+zr,zi,&sr,&si);
643e12c5d1SDavid du Colombier csq(sr,si,&tr,&ti);
653e12c5d1SDavid du Colombier ccubrt(1+3*tr,3*ti,&ur,&ui);
663e12c5d1SDavid du Colombier csqrt(ur-1,ui,&vr,&vi);
673e12c5d1SDavid du Colombier cdiv(rootroot3+vr,vi,rootroot3-vr,-vi,&yr,&yi);
683e12c5d1SDavid du Colombier yr /= rootk;
693e12c5d1SDavid du Colombier yi /= rootk;
703e12c5d1SDavid du Colombier elco2(fabs(yr),yi,hkc,1.,1.,x,y);
713e12c5d1SDavid du Colombier if(yr < 0)
723e12c5d1SDavid du Colombier *x = w2 - *x;
733e12c5d1SDavid du Colombier if(!ns) reflect(hcut[0]>place->wlon.l?0:
743e12c5d1SDavid du Colombier hcut[1]>=place->wlon.l?1:
753e12c5d1SDavid du Colombier 2,*x,*y,x,y);
763e12c5d1SDavid du Colombier return(1);
773e12c5d1SDavid du Colombier }
783e12c5d1SDavid du Colombier
793e12c5d1SDavid du Colombier proj
hex(void)803e12c5d1SDavid du Colombier hex(void)
813e12c5d1SDavid du Colombier {
823e12c5d1SDavid du Colombier register i;
83219b2ee8SDavid du Colombier double t;
843e12c5d1SDavid du Colombier double root3;
85219b2ee8SDavid du Colombier double c,d;
863e12c5d1SDavid du Colombier struct place p;
873e12c5d1SDavid du Colombier hcut[2] = PI;
883e12c5d1SDavid du Colombier hcut[1] = hcut[2]/3;
893e12c5d1SDavid du Colombier hcut[0] = -hcut[1];
903e12c5d1SDavid du Colombier root3 = sqrt(3.);
913e12c5d1SDavid du Colombier rootroot3 = sqrt(root3);
923e12c5d1SDavid du Colombier t = 15 -8*root3;
933e12c5d1SDavid du Colombier hkc = t*(1-sqrt(1-1/(t*t)));
943e12c5d1SDavid du Colombier elco2(BIG,0.,hkc,1.,1.,&w2,&t);
953e12c5d1SDavid du Colombier w2 *= 2;
963e12c5d1SDavid du Colombier rootk = sqrt(hkc);
973e12c5d1SDavid du Colombier latlon(90.,90.,&hem);
983e12c5d1SDavid du Colombier latlon(90.,0.,&p);
993e12c5d1SDavid du Colombier Xhex(&p,&c,&t);
1003e12c5d1SDavid du Colombier latlon(0.,0.,&p);
1013e12c5d1SDavid du Colombier Xhex(&p,&d,&t);
1023e12c5d1SDavid du Colombier for(i=0;i<3;i++) {
1033e12c5d1SDavid du Colombier ki[i] *= root3/2;
1043e12c5d1SDavid du Colombier cr[i] = c + (c-d)*kr[i];
1053e12c5d1SDavid du Colombier ci[i] = (c-d)*ki[i];
1063e12c5d1SDavid du Colombier }
1073e12c5d1SDavid du Colombier deg2rad(0.,&twist);
1083e12c5d1SDavid du Colombier return(Xhex);
1093e12c5d1SDavid du Colombier }
1103e12c5d1SDavid du Colombier
1113e12c5d1SDavid du Colombier int
hexcut(struct place * g,struct place * og,double * cutlon)112219b2ee8SDavid du Colombier hexcut(struct place *g, struct place *og, double *cutlon)
1133e12c5d1SDavid du Colombier {
1143e12c5d1SDavid du Colombier int t,i;
1153e12c5d1SDavid du Colombier if(g->nlat.l>=-HFUZZ&&og->nlat.l>=-HFUZZ)
1163e12c5d1SDavid du Colombier return(1);
1173e12c5d1SDavid du Colombier for(i=0;i<3;i++) {
1183e12c5d1SDavid du Colombier t = ckcut(g,og,*cutlon=hcut[i]);
1193e12c5d1SDavid du Colombier if(t!=1) return(t);
1203e12c5d1SDavid du Colombier }
1213e12c5d1SDavid du Colombier return(1);
1223e12c5d1SDavid du Colombier }
123