diff options
Diffstat (limited to 'src/cmd/map/libmap')
40 files changed, 2116 insertions, 0 deletions
diff --git a/src/cmd/map/libmap/aitoff.c b/src/cmd/map/libmap/aitoff.c new file mode 100644 index 00000000..83b777fb --- /dev/null +++ b/src/cmd/map/libmap/aitoff.c @@ -0,0 +1,26 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +#define Xaitwist Xaitpole.nlat +static struct place Xaitpole; + +static int +Xaitoff(struct place *place, double *x, double *y) +{ + struct place p; + copyplace(place,&p); + p.wlon.l /= 2.; + sincos(&p.wlon); + norm(&p,&Xaitpole,&Xaitwist); + Xazequalarea(&p,x,y); + *x *= 2.; + return(1); +} + +proj +aitoff(void) +{ + latlon(0.,0.,&Xaitpole); + return(Xaitoff); +} diff --git a/src/cmd/map/libmap/albers.c b/src/cmd/map/libmap/albers.c new file mode 100644 index 00000000..56e3e9f5 --- /dev/null +++ b/src/cmd/map/libmap/albers.c @@ -0,0 +1,117 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +/* For Albers formulas see Deetz and Adams "Elements of Map Projection", */ +/* USGS Special Publication No. 68, GPO 1921 */ + +static double r0sq, r1sq, d2, n, den, sinb1, sinb2; +static struct coord plat1, plat2; +static int southpole; + +static double num(double s) +{ + if(d2==0) + return(1); + s = d2*s*s; + return(1+s*(2./3+s*(3./5+s*(4./7+s*5./9)))); +} + +/* Albers projection for a spheroid, good only when N pole is fixed */ + +static int +Xspalbers(struct place *place, double *x, double *y) +{ + double r = sqrt(r0sq-2*(1-d2)*place->nlat.s*num(place->nlat.s)/n); + double t = n*place->wlon.l; + *y = r*cos(t); + *x = -r*sin(t); + if(!southpole) + *y = -*y; + else + *x = -*x; + return(1); +} + +/* lat1, lat2: std parallels; e2: squared eccentricity */ + +static proj albinit(double lat1, double lat2, double e2) +{ + double r1; + double t; + for(;;) { + if(lat1 < -90) + lat1 = -180 - lat1; + if(lat2 > 90) + lat2 = 180 - lat2; + if(lat1 <= lat2) + break; + t = lat1; lat1 = lat2; lat2 = t; + } + if(lat2-lat1 < 1) { + if(lat1 > 89) + return(azequalarea()); + return(0); + } + if(fabs(lat2+lat1) < 1) + return(cylequalarea(lat1)); + d2 = e2; + den = num(1.); + deg2rad(lat1,&plat1); + deg2rad(lat2,&plat2); + sinb1 = plat1.s*num(plat1.s)/den; + sinb2 = plat2.s*num(plat2.s)/den; + n = (plat1.c*plat1.c/(1-e2*plat1.s*plat1.s) - + plat2.c*plat2.c/(1-e2*plat2.s*plat2.s)) / + (2*(1-e2)*den*(sinb2-sinb1)); + r1 = plat1.c/(n*sqrt(1-e2*plat1.s*plat1.s)); + r1sq = r1*r1; + r0sq = r1sq + 2*(1-e2)*den*sinb1/n; + southpole = lat1<0 && plat2.c>plat1.c; + return(Xspalbers); +} + +proj +sp_albers(double lat1, double lat2) +{ + return(albinit(lat1,lat2,EC2)); +} + +proj +albers(double lat1, double lat2) +{ + return(albinit(lat1,lat2,0.)); +} + +static double scale = 1; +static double twist = 0; + +void +albscale(double x, double y, double lat, double lon) +{ + struct place place; + double alat, alon, x1,y1; + scale = 1; + twist = 0; + invalb(x,y,&alat,&alon); + twist = lon - alon; + deg2rad(lat,&place.nlat); + deg2rad(lon,&place.wlon); + Xspalbers(&place,&x1,&y1); + scale = sqrt((x1*x1+y1*y1)/(x*x+y*y)); +} + +void +invalb(double x, double y, double *lat, double *lon) +{ + int i; + double sinb_den, sinp; + x *= scale; + y *= scale; + *lon = atan2(-x,fabs(y))/(RAD*n) + twist; + sinb_den = (r0sq - x*x - y*y)*n/(2*(1-d2)); + sinp = sinb_den; + for(i=0; i<5; i++) + sinp = sinb_den/num(sinp); + *lat = asin(sinp)/RAD; +} diff --git a/src/cmd/map/libmap/azequalarea.c b/src/cmd/map/libmap/azequalarea.c new file mode 100644 index 00000000..6bae893d --- /dev/null +++ b/src/cmd/map/libmap/azequalarea.c @@ -0,0 +1,19 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +int +Xazequalarea(struct place *place, double *x, double *y) +{ + double r; + r = sqrt(1. - place->nlat.s); + *x = - r * place->wlon.s; + *y = - r * place->wlon.c; + return(1); +} + +proj +azequalarea(void) +{ + return(Xazequalarea); +} diff --git a/src/cmd/map/libmap/azequidist.c b/src/cmd/map/libmap/azequidist.c new file mode 100644 index 00000000..d26d33d4 --- /dev/null +++ b/src/cmd/map/libmap/azequidist.c @@ -0,0 +1,19 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +int +Xazequidistant(struct place *place, double *x, double *y) +{ + double colat; + colat = PI/2 - place->nlat.l; + *x = -colat * place->wlon.s; + *y = -colat * place->wlon.c; + return(1); +} + +proj +azequidistant(void) +{ + return(Xazequidistant); +} diff --git a/src/cmd/map/libmap/bicentric.c b/src/cmd/map/libmap/bicentric.c new file mode 100644 index 00000000..33fd8d19 --- /dev/null +++ b/src/cmd/map/libmap/bicentric.c @@ -0,0 +1,25 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static struct coord center; + +static int +Xbicentric(struct place *place, double *x, double *y) +{ + if(place->wlon.c<=.01||place->nlat.c<=.01) + return(-1); + *x = -center.c*place->wlon.s/place->wlon.c; + *y = place->nlat.s/(place->nlat.c*place->wlon.c); + return(*x**x+*y**y<=9); +} + +proj +bicentric(double l) +{ + l = fabs(l); + if(l>89) + return(0); + deg2rad(l,¢er); + return(Xbicentric); +} diff --git a/src/cmd/map/libmap/bonne.c b/src/cmd/map/libmap/bonne.c new file mode 100644 index 00000000..858f0d69 --- /dev/null +++ b/src/cmd/map/libmap/bonne.c @@ -0,0 +1,36 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static struct coord stdpar; +static double r0; + +static int +Xbonne(struct place *place, double *x, double *y) +{ + double r, alpha; + r = r0 - place->nlat.l; + if(r<.001) + if(fabs(stdpar.c)<1e-10) + alpha = place->wlon.l; + else if(fabs(place->nlat.c)==0) + alpha = 0; + else + alpha = place->wlon.l/(1+ + stdpar.c*stdpar.c*stdpar.c/place->nlat.c/3); + else + alpha = place->wlon.l * place->nlat.c / r; + *x = - r*sin(alpha); + *y = - r*cos(alpha); + return(1); +} + +proj +bonne(double par) +{ + if(fabs(par*RAD) < .01) + return(Xsinusoidal); + deg2rad(par, &stdpar); + r0 = stdpar.c/stdpar.s + stdpar.l; + return(Xbonne); +} diff --git a/src/cmd/map/libmap/ccubrt.c b/src/cmd/map/libmap/ccubrt.c new file mode 100644 index 00000000..8a7af566 --- /dev/null +++ b/src/cmd/map/libmap/ccubrt.c @@ -0,0 +1,13 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +void +ccubrt(double zr, double zi, double *wr, double *wi) +{ + double r, theta; + theta = atan2(zi,zr); + r = cubrt(hypot(zr,zi)); + *wr = r*cos(theta/3); + *wi = r*sin(theta/3); +} diff --git a/src/cmd/map/libmap/complex.c b/src/cmd/map/libmap/complex.c new file mode 100644 index 00000000..b3099fd4 --- /dev/null +++ b/src/cmd/map/libmap/complex.c @@ -0,0 +1,85 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +/*complex divide, defensive against overflow from + * * and /, but not from + and - + * assumes underflow yields 0.0 + * uses identities: + * (a + bi)/(c + di) = ((a + bd/c) + (b - ad/c)i)/(c + dd/c) + * (a + bi)/(c + di) = (b - ai)/(d - ci) +*/ +void +cdiv(double a, double b, double c, double d, double *u, double *v) +{ + double r,t; + if(fabs(c)<fabs(d)) { + t = -c; c = d; d = t; + t = -a; a = b; b = t; + } + r = d/c; + t = c + r*d; + *u = (a + r*b)/t; + *v = (b - r*a)/t; +} + +void +cmul(double c1, double c2, double d1, double d2, double *e1, double *e2) +{ + *e1 = c1*d1 - c2*d2; + *e2 = c1*d2 + c2*d1; +} + +void +csq(double c1, double c2, double *e1, double *e2) +{ + *e1 = c1*c1 - c2*c2; + *e2 = c1*c2*2; +} + +/* complex square root + * assumes underflow yields 0.0 + * uses these identities: + * sqrt(x+_iy) = sqrt(r(cos(t)+_isin(t)) + * = sqrt(r)(cos(t/2)+_isin(t/2)) + * cos(t/2) = sin(t)/2sin(t/2) = sqrt((1+cos(t)/2) + * sin(t/2) = sin(t)/2cos(t/2) = sqrt((1-cos(t)/2) +*/ +void +csqrt(double c1, double c2, double *e1, double *e2) +{ + double r,s; + double x,y; + x = fabs(c1); + y = fabs(c2); + if(x>=y) { + if(x==0) { + *e1 = *e2 = 0; + return; + } + r = x; + s = y/x; + } else { + r = y; + s = x/y; + } + r *= sqrt(1+ s*s); + if(c1>0) { + *e1 = sqrt((r+c1)/2); + *e2 = c2/(2* *e1); + } else { + *e2 = sqrt((r-c1)/2); + if(c2<0) + *e2 = -*e2; + *e1 = c2/(2* *e2); + } +} + + +void cpow(double c1, double c2, double *d1, double *d2, double pwr) +{ + double theta = pwr*atan2(c2,c1); + double r = pow(hypot(c1,c2), pwr); + *d1 = r*cos(theta); + *d2 = r*sin(theta); +} diff --git a/src/cmd/map/libmap/conic.c b/src/cmd/map/libmap/conic.c new file mode 100644 index 00000000..ba4430c6 --- /dev/null +++ b/src/cmd/map/libmap/conic.c @@ -0,0 +1,27 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static struct coord stdpar; + +static int +Xconic(struct place *place, double *x, double *y) +{ + double r; + if(fabs(place->nlat.l-stdpar.l) > 80.*RAD) + return(-1); + r = stdpar.c/stdpar.s - tan(place->nlat.l - stdpar.l); + *x = - r*sin(place->wlon.l * stdpar.s); + *y = - r*cos(place->wlon.l * stdpar.s); + if(r>3) return(0); + return(1); +} + +proj +conic(double par) +{ + if(fabs(par) <.1) + return(Xcylindrical); + deg2rad(par, &stdpar); + return(Xconic); +} diff --git a/src/cmd/map/libmap/cubrt.c b/src/cmd/map/libmap/cubrt.c new file mode 100644 index 00000000..fd508d2b --- /dev/null +++ b/src/cmd/map/libmap/cubrt.c @@ -0,0 +1,30 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +double +cubrt(double a) +{ + double x,y,x1; + if(a==0) + return(0.); + y = 1; + if(a<0) { + y = -y; + a = -a; + } + while(a<1) { + a *= 8; + y /= 2; + } + while(a>1) { + a /= 8; + y *= 2; + } + x = 1; + do { + x1 = x; + x = (2*x1+a/(x1*x1))/3; + } while(fabs(x-x1)>10.e-15); + return(x*y); +} diff --git a/src/cmd/map/libmap/cuts.c b/src/cmd/map/libmap/cuts.c new file mode 100644 index 00000000..ad9d6dc3 --- /dev/null +++ b/src/cmd/map/libmap/cuts.c @@ -0,0 +1,39 @@ +#include <u.h> +#include <libc.h> +#include "map.h" +extern void abort(void); + +/* these routines duplicate names found in map.c. they are +called from routines in hex.c, guyou.c, and tetra.c, which +are in turn invoked directly from map.c. this bad organization +arises from data hiding; only these three files know stuff +that's necessary for the proper handling of the unusual cuts +involved in these projections. + +the calling routines are not advertised as part of the library, +and the library duplicates should never get loaded, however they +are included to make the libary self-standing.*/ + +int +picut(struct place *g, struct place *og, double *cutlon) +{ + g; og; cutlon; + abort(); + return 0; +} + +int +ckcut(struct place *g1, struct place *g2, double lon) +{ + g1; g2; lon; + abort(); + return 0; +} + +double +reduce(double x) +{ + x; + abort(); + return 0; +} diff --git a/src/cmd/map/libmap/cylequalarea.c b/src/cmd/map/libmap/cylequalarea.c new file mode 100644 index 00000000..3cf222ff --- /dev/null +++ b/src/cmd/map/libmap/cylequalarea.c @@ -0,0 +1,24 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static double a; + +static int +Xcylequalarea(struct place *place, double *x, double *y) +{ + *x = - place->wlon.l * a; + *y = place->nlat.s; + return(1); +} + +proj +cylequalarea(double par) +{ + struct coord stdp0; + if(par > 89.0) + return(0); + deg2rad(par, &stdp0); + a = stdp0.c*stdp0.c; + return(Xcylequalarea); +} diff --git a/src/cmd/map/libmap/cylindrical.c b/src/cmd/map/libmap/cylindrical.c new file mode 100644 index 00000000..4d01bc23 --- /dev/null +++ b/src/cmd/map/libmap/cylindrical.c @@ -0,0 +1,19 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +int +Xcylindrical(struct place *place, double *x, double *y) +{ + if(fabs(place->nlat.l) > 80.*RAD) + return(-1); + *x = - place->wlon.l; + *y = place->nlat.s / place->nlat.c; + return(1); +} + +proj +cylindrical(void) +{ + return(Xcylindrical); +} diff --git a/src/cmd/map/libmap/elco2.c b/src/cmd/map/libmap/elco2.c new file mode 100644 index 00000000..b4c9bbf6 --- /dev/null +++ b/src/cmd/map/libmap/elco2.c @@ -0,0 +1,132 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +/* elliptic integral routine, R.Bulirsch, + * Numerische Mathematik 7(1965) 78-90 + * calculate integral from 0 to x+iy of + * (a+b*t^2)/((1+t^2)*sqrt((1+t^2)*(1+kc^2*t^2))) + * yields about D valid figures, where CC=10e-D + * for a*b>=0, except at branchpoints x=0,y=+-i,+-i/kc; + * there the accuracy may be reduced. + * fails for kc=0 or x<0 + * return(1) for success, return(0) for fail + * + * special case a=b=1 is equivalent to + * standard elliptic integral of first kind + * from 0 to atan(x+iy) of + * 1/sqrt(1-k^2*(sin(t))^2) where k^2=1-kc^2 +*/ + +#define ROOTINF 10.e18 +#define CC 1.e-6 + +int +elco2(double x, double y, double kc, double a, double b, double *u, double *v) +{ + double c,d,dn1,dn2,e,e1,e2,f,f1,f2,h,k,m,m1,m2,sy; + double d1[13],d2[13]; + int i,l; + if(kc==0||x<0) + return(0); + sy = y>0? 1: y==0? 0: -1; + y = fabs(y); + csq(x,y,&c,&e2); + d = kc*kc; + k = 1-d; + e1 = 1+c; + cdiv2(1+d*c,d*e2,e1,e2,&f1,&f2); + f2 = -k*x*y*2/f2; + csqr(f1,f2,&dn1,&dn2); + if(f1<0) { + f1 = dn1; + dn1 = -dn2; + dn2 = -f1; + } + if(k<0) { + dn1 = fabs(dn1); + dn2 = fabs(dn2); + } + c = 1+dn1; + cmul(e1,e2,c,dn2,&f1,&f2); + cdiv(x,y,f1,f2,&d1[0],&d2[0]); + h = a-b; + d = f = m = 1; + kc = fabs(kc); + e = a; + a += b; + l = 4; + for(i=1;;i++) { + m1 = (kc+m)/2; + m2 = m1*m1; + k *= f/(m2*4); + b += e*kc; + e = a; + cdiv2(kc+m*dn1,m*dn2,c,dn2,&f1,&f2); + csqr(f1/m1,k*dn2*2/f2,&dn1,&dn2); + cmul(dn1,dn2,x,y,&f1,&f2); + x = fabs(f1); + y = fabs(f2); + a += b/m1; + l *= 2; + c = 1 +dn1; + d *= k/2; + cmul(x,y,x,y,&e1,&e2); + k *= k; + + cmul(c,dn2,1+e1*m2,e2*m2,&f1,&f2); + cdiv(d*x,d*y,f1,f2,&d1[i],&d2[i]); + if(k<=CC) + break; + kc = sqrt(m*kc); + f = m2; + m = m1; + } + f1 = f2 = 0; + for(;i>=0;i--) { + f1 += d1[i]; + f2 += d2[i]; + } + x *= m1; + y *= m1; + cdiv2(1-y,x,1+y,-x,&e1,&e2); + e2 = x*2/e2; + d = a/(m1*l); + *u = atan2(e2,e1); + if(*u<0) + *u += PI; + a = d*sy/2; + *u = d*(*u) + f1*h; + *v = (-1-log(e1*e1+e2*e2))*a + f2*h*sy + a; + return(1); +} + +void +cdiv2(double c1, double c2, double d1, double d2, double *e1, double *e2) +{ + double t; + if(fabs(d2)>fabs(d1)) { + t = d1, d1 = d2, d2 = t; + t = c1, c1 = c2, c2 = t; + } + if(fabs(d1)>ROOTINF) + *e2 = ROOTINF*ROOTINF; + else + *e2 = d1*d1 + d2*d2; + t = d2/d1; + *e1 = (c1+t*c2)/(d1+t*d2); /* (c1*d1+c2*d2)/(d1*d1+d2*d2) */ +} + +/* complex square root of |x|+iy */ +void +csqr(double c1, double c2, double *e1, double *e2) +{ + double r2; + r2 = c1*c1 + c2*c2; + if(r2<=0) { + *e1 = *e2 = 0; + return; + } + *e1 = sqrt((sqrt(r2) + fabs(c1))/2); + *e2 = c2/(*e1*2); +} diff --git a/src/cmd/map/libmap/elliptic.c b/src/cmd/map/libmap/elliptic.c new file mode 100644 index 00000000..3f3b1d57 --- /dev/null +++ b/src/cmd/map/libmap/elliptic.c @@ -0,0 +1,35 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +struct coord center; + +static int +Xelliptic(struct place *place, double *x, double *y) +{ + double r1,r2; + r1 = acos(place->nlat.c*(place->wlon.c*center.c + - place->wlon.s*center.s)); + r2 = acos(place->nlat.c*(place->wlon.c*center.c + + place->wlon.s*center.s)); + *x = -(r1*r1 - r2*r2)/(4*center.l); + *y = (r1*r1+r2*r2)/2 - (center.l*center.l+*x**x); + if(*y < 0) + *y = 0; + *y = sqrt(*y); + if(place->nlat.l<0) + *y = -*y; + return(1); +} + +proj +elliptic(double l) +{ + l = fabs(l); + if(l>89) + return(0); + if(l<1) + return(Xazequidistant); + deg2rad(l,¢er); + return(Xelliptic); +} diff --git a/src/cmd/map/libmap/fisheye.c b/src/cmd/map/libmap/fisheye.c new file mode 100644 index 00000000..412d65e7 --- /dev/null +++ b/src/cmd/map/libmap/fisheye.c @@ -0,0 +1,26 @@ +#include <u.h> +#include <libc.h> +#include "map.h" +/* refractive fisheye, not logarithmic */ + +static double n; + +static int +Xfisheye(struct place *place, double *x, double *y) +{ + double r; + double u = sin(PI/4-place->nlat.l/2)/n; + if(fabs(u) > .97) + return -1; + r = tan(asin(u)); + *x = -r*place->wlon.s; + *y = -r*place->wlon.c; + return 1; +} + +proj +fisheye(double par) +{ + n = par; + return n<.1? 0: Xfisheye; +} diff --git a/src/cmd/map/libmap/gall.c b/src/cmd/map/libmap/gall.c new file mode 100644 index 00000000..84da651b --- /dev/null +++ b/src/cmd/map/libmap/gall.c @@ -0,0 +1,29 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static double scale; + +static int +Xgall(struct place *place, double *x, double *y) +{ + /* two ways to compute tan(place->nlat.l/2) */ + if(fabs(place->nlat.s)<.1) + *y = sin(place->nlat.l/2)/cos(place->nlat.l/2); + else + *y = (1-place->nlat.c)/place->nlat.s; + *x = -scale*place->wlon.l; + return 1; +} + +proj +gall(double par) +{ + double coshalf; + if(fabs(par)>80) + return 0; + par *= RAD; + coshalf = cos(par/2); + scale = cos(par)/(2*coshalf*coshalf); + return Xgall; +} diff --git a/src/cmd/map/libmap/gilbert.c b/src/cmd/map/libmap/gilbert.c new file mode 100644 index 00000000..173ffcd3 --- /dev/null +++ b/src/cmd/map/libmap/gilbert.c @@ -0,0 +1,51 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +int +Xgilbert(struct place *p, double *x, double *y) +{ +/* the interesting part - map the sphere onto a hemisphere */ + struct place q; + q.nlat.s = tan(0.5*(p->nlat.l)); + if(q.nlat.s > 1) q.nlat.s = 1; + if(q.nlat.s < -1) q.nlat.s = -1; + q.nlat.c = sqrt(1 - q.nlat.s*q.nlat.s); + q.wlon.l = p->wlon.l/2; + sincos(&q.wlon); +/* the dull part: present the hemisphere orthogrpahically */ + *y = q.nlat.s; + *x = -q.wlon.s*q.nlat.c; + return(1); +} + +proj +gilbert(void) +{ + return(Xgilbert); +} + +/* derivation of the interesting part: + map the sphere onto the plane by stereographic projection; + map the plane onto a half plane by sqrt; + map the half plane back to the sphere by stereographic + projection + + n,w are original lat and lon + r is stereographic radius + primes are transformed versions + + r = cos(n)/(1+sin(n)) + r' = sqrt(r) = cos(n')/(1+sin(n')) + + r'^2 = (1-sin(n')^2)/(1+sin(n')^2) = cos(n)/(1+sin(n)) + + this is a linear equation for sin n', with solution + + sin n' = (1+sin(n)-cos(n))/(1+sin(n)+cos(n)) + + use standard formula: tan x/2 = (1-cos x)/sin x = sin x/(1+cos x) + to show that the right side of the last equation is tan(n/2) +*/ + + 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); +} diff --git a/src/cmd/map/libmap/harrison.c b/src/cmd/map/libmap/harrison.c new file mode 100644 index 00000000..6b6003c2 --- /dev/null +++ b/src/cmd/map/libmap/harrison.c @@ -0,0 +1,40 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static double v3,u2,u3,a,b; /*v=view,p=obj,u=unit.y*/ + +static int +Xharrison(struct place *place, double *x, double *y) +{ + double p1 = -place->nlat.c*place->wlon.s; + double p2 = -place->nlat.c*place->wlon.c; + double p3 = place->nlat.s; + double d = b + u3*p2 - u2*p3; + double t; + if(d < .01) + return -1; + t = a/d; + if(v3*place->nlat.s < 1.) + return -1; + *y = t*p2*u2 + (v3-t*(v3-p3))*u3; + *x = t*p1; + if(t < 0) + return 0; + if(*x * *x + *y * *y > 16) + return -1; + return 1; +} + +proj +harrison(double r, double alpha) +{ + u2 = cos(alpha*RAD); + u3 = sin(alpha*RAD); + v3 = r; + b = r*u2; + a = 1 + b; + if(r<1.001 || a<sqrt(r*r-1)) + return 0; + return Xharrison; +} diff --git a/src/cmd/map/libmap/hex.c b/src/cmd/map/libmap/hex.c new file mode 100644 index 00000000..851f138f --- /dev/null +++ b/src/cmd/map/libmap/hex.c @@ -0,0 +1,122 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +#define BIG 1.e15 +#define HFUZZ .0001 + +static double hcut[3] ; +static double kr[3] = { .5, -1., .5 }; +static double ki[3] = { -1., 0., 1. }; /*to multiply by sqrt(3)/2*/ +static double cr[3]; +static double ci[3]; +static struct place hem; +static struct coord twist; +static double rootroot3, hkc; +static double w2; +static double rootk; + +static void +reflect(int i, double wr, double wi, double *x, double *y) +{ + double pr,pi,l; + pr = cr[i]-wr; + pi = ci[i]-wi; + l = 2*(kr[i]*pr + ki[i]*pi); + *x = wr + l*kr[i]; + *y = wi + l*ki[i]; +} + +static int +Xhex(struct place *place, double *x, double *y) +{ + int ns; + int i; + double zr,zi; + double sr,si,tr,ti,ur,ui,vr,vi,yr,yi; + struct place p; + copyplace(place,&p); + ns = place->nlat.l >= 0; + if(!ns) { + p.nlat.l = -p.nlat.l; + p.nlat.s = -p.nlat.s; + } + if(p.nlat.l<HFUZZ) { + for(i=0;i<3;i++) + if(fabs(reduce(p.wlon.l-hcut[i]))<HFUZZ) { + if(i==2) { + *x = 2*cr[0] - cr[1]; + *y = 0; + } else { + *x = cr[1]; + *y = 2*ci[2*i]; + } + return(1); + } + p.nlat.l = HFUZZ; + sincos(&p.nlat); + } + norm(&p,&hem,&twist); + Xstereographic(&p,&zr,&zi); + zr /= 2; + zi /= 2; + cdiv(1-zr,-zi,1+zr,zi,&sr,&si); + csq(sr,si,&tr,&ti); + ccubrt(1+3*tr,3*ti,&ur,&ui); + csqrt(ur-1,ui,&vr,&vi); + cdiv(rootroot3+vr,vi,rootroot3-vr,-vi,&yr,&yi); + yr /= rootk; + yi /= rootk; + elco2(fabs(yr),yi,hkc,1.,1.,x,y); + if(yr < 0) + *x = w2 - *x; + if(!ns) reflect(hcut[0]>place->wlon.l?0: + hcut[1]>=place->wlon.l?1: + 2,*x,*y,x,y); + return(1); +} + +proj +hex(void) +{ + int i; + double t; + double root3; + double c,d; + struct place p; + hcut[2] = PI; + hcut[1] = hcut[2]/3; + hcut[0] = -hcut[1]; + root3 = sqrt(3.); + rootroot3 = sqrt(root3); + t = 15 -8*root3; + hkc = t*(1-sqrt(1-1/(t*t))); + elco2(BIG,0.,hkc,1.,1.,&w2,&t); + w2 *= 2; + rootk = sqrt(hkc); + latlon(90.,90.,&hem); + latlon(90.,0.,&p); + Xhex(&p,&c,&t); + latlon(0.,0.,&p); + Xhex(&p,&d,&t); + for(i=0;i<3;i++) { + ki[i] *= root3/2; + cr[i] = c + (c-d)*kr[i]; + ci[i] = (c-d)*ki[i]; + } + deg2rad(0.,&twist); + return(Xhex); +} + +int +hexcut(struct place *g, struct place *og, double *cutlon) +{ + int t,i; + if(g->nlat.l>=-HFUZZ&&og->nlat.l>=-HFUZZ) + return(1); + for(i=0;i<3;i++) { + t = ckcut(g,og,*cutlon=hcut[i]); + if(t!=1) return(t); + } + return(1); +} diff --git a/src/cmd/map/libmap/homing.c b/src/cmd/map/libmap/homing.c new file mode 100644 index 00000000..366f69fe --- /dev/null +++ b/src/cmd/map/libmap/homing.c @@ -0,0 +1,121 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static struct coord p0; /* standard parallel */ + +int first; + +static double +trigclamp(double x) +{ + return x>1? 1: x<-1? -1: x; +} + +static struct coord az; /* azimuth of p0 seen from place */ +static struct coord rad; /* angular dist from place to p0 */ + +static int +azimuth(struct place *place) +{ + if(place->nlat.c < FUZZ) { + az.l = PI/2 + place->nlat.l - place->wlon.l; + sincos(&az); + rad.l = fabs(place->nlat.l - p0.l); + if(rad.l > PI) + rad.l = 2*PI - rad.l; + sincos(&rad); + return 1; + } + rad.c = trigclamp(p0.s*place->nlat.s + /* law of cosines */ + p0.c*place->nlat.c*place->wlon.c); + rad.s = sqrt(1 - rad.c*rad.c); + if(fabs(rad.s) < .001) { + az.s = 0; + az.c = 1; + } else { + az.s = trigclamp(p0.c*place->wlon.s/rad.s); /* sines */ + az.c = trigclamp((p0.s - rad.c*place->nlat.s) + /(rad.s*place->nlat.c)); + } + rad.l = atan2(rad.s, rad.c); + return 1; +} + +static int +Xmecca(struct place *place, double *x, double *y) +{ + if(!azimuth(place)) + return 0; + *x = -place->wlon.l; + *y = fabs(az.s)<.02? -az.c*rad.s/p0.c: *x*az.c/az.s; + return fabs(*y)>2? -1: + rad.c<0? 0: + 1; +} + +proj +mecca(double par) +{ + first = 1; + if(fabs(par)>80.) + return(0); + deg2rad(par,&p0); + return(Xmecca); +} + +static int +Xhoming(struct place *place, double *x, double *y) +{ + if(!azimuth(place)) + return 0; + *x = -rad.l*az.s; + *y = -rad.l*az.c; + return place->wlon.c<0? 0: 1; +} + +proj +homing(double par) +{ + first = 1; + if(fabs(par)>80.) + return(0); + deg2rad(par,&p0); + return(Xhoming); +} + +int +hlimb(double *lat, double *lon, double res) +{ + if(first) { + *lon = -90; + *lat = -90; + first = 0; + return 0; + } + *lat += res; + if(*lat <= 90) + return 1; + if(*lon == 90) + return -1; + *lon = 90; + *lat = -90; + return 0; +} + +int +mlimb(double *lat, double *lon, double res) +{ + int ret = !first; + if(fabs(p0.s) < .01) + return -1; + if(first) { + *lon = -180; + first = 0; + } else + *lon += res; + if(*lon > 180) + return -1; + *lat = atan(-cos(*lon*RAD)/p0.s*p0.c)/RAD; + return ret; +} diff --git a/src/cmd/map/libmap/lagrange.c b/src/cmd/map/libmap/lagrange.c new file mode 100644 index 00000000..02dd29eb --- /dev/null +++ b/src/cmd/map/libmap/lagrange.c @@ -0,0 +1,30 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static int +Xlagrange(struct place *place, double *x, double *y) +{ + double z1,z2; + double w1,w2,t1,t2; + struct place p; + copyplace(place,&p); + if(place->nlat.l<0) { + p.nlat.l = -p.nlat.l; + p.nlat.s = -p.nlat.s; + } + Xstereographic(&p,&z1,&z2); + csqrt(-z2/2,z1/2,&w1,&w2); + cdiv(w1-1,w2,w1+1,w2,&t1,&t2); + *y = -t1; + *x = t2; + if(place->nlat.l<0) + *y = -*y; + return(1); +} + +proj +lagrange(void) +{ + return(Xlagrange); +} diff --git a/src/cmd/map/libmap/lambert.c b/src/cmd/map/libmap/lambert.c new file mode 100644 index 00000000..6b688aa3 --- /dev/null +++ b/src/cmd/map/libmap/lambert.c @@ -0,0 +1,46 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static struct coord stdp0, stdp1; +static double k; + +static int +Xlambert(struct place *place, double *x, double *y) +{ + double r; + if(place->nlat.l < -80.*RAD) + return(-1); + if(place->nlat.l > 89.*RAD) + r = 0; /* slovenly */ + else + r = stdp0.c*exp(0.5*k*log( + (1+stdp0.s)*(1-place->nlat.s)/((1-stdp0.s)*(1+place->nlat.s)))); + if(stdp1.l<0.) + r = -r; + *x = - r*sin(k * place->wlon.l); + *y = - r*cos(k * place->wlon.l); + return(1); +} + +proj +lambert(double par0, double par1) +{ + double temp; + if(fabs(par0)>fabs(par1)){ + temp = par0; + par0 = par1; + par1 = temp; + } + deg2rad(par0, &stdp0); + deg2rad(par1, &stdp1); + if(fabs(par1+par0)<.1) + return(mercator()); + if(fabs(par1-par0)<.1) + return(perspective(-1.)); + if(fabs(par0)>89.5||fabs(par1)>89.5) + return(0); + k = 2*log(stdp1.c/stdp0.c)/log( + (1+stdp0.s)*(1-stdp1.s)/((1-stdp0.s)*(1+stdp1.s))); + return(Xlambert); +} diff --git a/src/cmd/map/libmap/laue.c b/src/cmd/map/libmap/laue.c new file mode 100644 index 00000000..06c5f3a7 --- /dev/null +++ b/src/cmd/map/libmap/laue.c @@ -0,0 +1,24 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + + +static int +Xlaue(struct place *place, double *x, double *y) +{ + double r; + if(place->nlat.l<PI/4+FUZZ) + return(-1); + r = tan(PI-2*place->nlat.l); + if(r>3) + return(-1); + *x = - r * place->wlon.s; + *y = - r * place->wlon.c; + return(1); +} + +proj +laue(void) +{ + return(Xlaue); +} diff --git a/src/cmd/map/libmap/lune.c b/src/cmd/map/libmap/lune.c new file mode 100644 index 00000000..dc58c6f0 --- /dev/null +++ b/src/cmd/map/libmap/lune.c @@ -0,0 +1,62 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +int Xstereographic(struct place *place, double *x, double *y); + +static struct place eastpole; +static struct place westpole; +static double eastx, easty; +static double westx, westy; +static double scale; +static double pwr; + +/* conformal map w = ((1+z)^A - (1-z)^A)/((1+z)^A + (1-z)^A), + where A<1, maps unit circle onto a convex lune with x= +-1 + mapping to vertices of angle A*PI at w = +-1 */ + +/* there are cuts from E and W poles to S pole, + in absence of a cut routine, error is returned for + points outside a polar cap through E and W poles */ + +static int Xlune(struct place *place, double *x, double *y) +{ + double stereox, stereoy; + double z1x, z1y, z2x, z2y; + double w1x, w1y, w2x, w2y; + double numx, numy, denx, deny; + if(place->nlat.l < eastpole.nlat.l-FUZZ) + return -1; + Xstereographic(place, &stereox, &stereoy); + stereox *= scale; + stereoy *= scale; + z1x = 1 + stereox; + z1y = stereoy; + z2x = 1 - stereox; + z2y = -stereoy; + cpow(z1x,z1y,&w1x,&w1y,pwr); + cpow(z2x,z2y,&w2x,&w2y,pwr); + numx = w1x - w2x; + numy = w1y - w2y; + denx = w1x + w2x; + deny = w1y + w2y; + cdiv(numx, numy, denx, deny, x, y); + return 1; +} + +proj +lune(double lat, double theta) +{ + deg2rad(lat, &eastpole.nlat); + deg2rad(-90.,&eastpole.wlon); + deg2rad(lat, &westpole.nlat); + deg2rad(90. ,&westpole.wlon); + Xstereographic(&eastpole, &eastx, &easty); + Xstereographic(&westpole, &westx, &westy); + if(fabs(easty)>FUZZ || fabs(westy)>FUZZ || + fabs(eastx+westx)>FUZZ) + abort(); + scale = 1/eastx; + pwr = theta/180; + return Xlune; +} diff --git a/src/cmd/map/libmap/mercator.c b/src/cmd/map/libmap/mercator.c new file mode 100644 index 00000000..0ab63653 --- /dev/null +++ b/src/cmd/map/libmap/mercator.c @@ -0,0 +1,36 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static int +Xmercator(struct place *place, double *x, double *y) +{ + if(fabs(place->nlat.l) > 80.*RAD) + return(-1); + *x = -place->wlon.l; + *y = 0.5*log((1+place->nlat.s)/(1-place->nlat.s)); + return(1); +} + +proj +mercator(void) +{ + return(Xmercator); +} + +static double ecc = ECC; + +static int +Xspmercator(struct place *place, double *x, double *y) +{ + if(Xmercator(place,x,y) < 0) + return(-1); + *y += 0.5*ecc*log((1-ecc*place->nlat.s)/(1+ecc*place->nlat.s)); + return(1); +} + +proj +sp_mercator(void) +{ + return(Xspmercator); +} diff --git a/src/cmd/map/libmap/mkfile b/src/cmd/map/libmap/mkfile new file mode 100644 index 00000000..e9031246 --- /dev/null +++ b/src/cmd/map/libmap/mkfile @@ -0,0 +1,50 @@ +<$PLAN9/src/mkhdr + +LIB=libmap.a +OFILES=aitoff.$O\ + albers.$O\ + azequalarea.$O\ + azequidist.$O\ + bicentric.$O\ + bonne.$O\ + ccubrt.$O\ + complex.$O\ + conic.$O\ + cubrt.$O\ + cylequalarea.$O\ + cylindrical.$O\ + elco2.$O\ + elliptic.$O\ + fisheye.$O\ + gall.$O\ + gilbert.$O\ + guyou.$O\ + harrison.$O\ + hex.$O\ + homing.$O\ + lagrange.$O\ + lambert.$O\ + laue.$O\ + lune.$O\ + mercator.$O\ + mollweide.$O\ + newyorker.$O\ + orthographic.$O\ + perspective.$O\ + polyconic.$O\ + rectangular.$O\ + simpleconic.$O\ + sinusoidal.$O\ + tetra.$O\ + trapezoidal.$O\ + twocirc.$O\ + zcoord.$O\ + +HFILES=../map.h\ + +<$PLAN9/src/mklib +CFLAGS=$CFLAGS -I.. + +nuke:V: + mk clean + rm -f libmap.a[$OS] diff --git a/src/cmd/map/libmap/mollweide.c b/src/cmd/map/libmap/mollweide.c new file mode 100644 index 00000000..3284c495 --- /dev/null +++ b/src/cmd/map/libmap/mollweide.c @@ -0,0 +1,25 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static int +Xmollweide(struct place *place, double *x, double *y) +{ + double z; + double w; + z = place->nlat.l; + if(fabs(z)<89.9*RAD) + do { /*newton for 2z+sin2z=pi*sin(lat)*/ + w = (2*z+sin(2*z)-PI*place->nlat.s)/(2+2*cos(2*z)); + z -= w; + } while(fabs(w)>=.00001); + *y = sin(z); + *x = - (2/PI)*cos(z)*place->wlon.l; + return(1); +} + +proj +mollweide(void) +{ + return(Xmollweide); +} diff --git a/src/cmd/map/libmap/newyorker.c b/src/cmd/map/libmap/newyorker.c new file mode 100644 index 00000000..370e3b37 --- /dev/null +++ b/src/cmd/map/libmap/newyorker.c @@ -0,0 +1,28 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static double a; + +static int +Xnewyorker(struct place *place, double *x, double *y) +{ + double r = PI/2 - place->nlat.l; + double s; + if(r<.001) /* cheat to plot center */ + s = 0; + else if(r<a) + return -1; + else + s = log(r/a); + *x = -s * place->wlon.s; + *y = -s * place->wlon.c; + return(1); +} + +proj +newyorker(double a0) +{ + a = a0*RAD; + return(Xnewyorker); +} diff --git a/src/cmd/map/libmap/orthographic.c b/src/cmd/map/libmap/orthographic.c new file mode 100644 index 00000000..7aac5b15 --- /dev/null +++ b/src/cmd/map/libmap/orthographic.c @@ -0,0 +1,35 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + + +int +Xorthographic(struct place *place, double *x, double *y) +{ + *x = - place->nlat.c * place->wlon.s; + *y = - place->nlat.c * place->wlon.c; + return(place->nlat.l<0.? 0 : 1); +} + +proj +orthographic(void) +{ + return(Xorthographic); +} + +int +olimb(double *lat, double *lon, double res) +{ + static int first = 1; + if(first) { + *lat = 0; + *lon = -180; + first = 0; + return 0; + } + *lon += res; + if(*lon <= 180) + return 1; + first = 1; + return -1; +} diff --git a/src/cmd/map/libmap/perspective.c b/src/cmd/map/libmap/perspective.c new file mode 100644 index 00000000..7ce6b0d6 --- /dev/null +++ b/src/cmd/map/libmap/perspective.c @@ -0,0 +1,84 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +#define ORTHRAD 1000 +static double viewpt; + +static int +Xperspective(struct place *place, double *x, double *y) +{ + double r; + if(viewpt<=1+FUZZ && fabs(place->nlat.s<=viewpt+.01)) + return(-1); + r = place->nlat.c*(viewpt - 1.)/(viewpt - place->nlat.s); + *x = - r*place->wlon.s; + *y = - r*place->wlon.c; + if(r>4.) + return(-1); + if(fabs(viewpt)>1 && place->nlat.s<1/viewpt || + fabs(viewpt)<=1 && place->nlat.s<viewpt) + return 0; + return(1); +} + +proj +perspective(double radius) +{ + viewpt = radius; + if(viewpt >= ORTHRAD) + return(Xorthographic); + if(fabs(viewpt-1.)<.0001) + return(0); + return(Xperspective); +} + + /* called from various conformal projections, + but not from stereographic itself */ +int +Xstereographic(struct place *place, double *x, double *y) +{ + double v = viewpt; + int retval; + viewpt = -1; + retval = Xperspective(place, x, y); + viewpt = v; + return retval; +} + +proj +stereographic(void) +{ + viewpt = -1.; + return(Xperspective); +} + +proj +gnomonic(void) +{ + viewpt = 0.; + return(Xperspective); +} + +int +plimb(double *lat, double *lon, double res) +{ + static int first = 1; + if(viewpt >= ORTHRAD) + return olimb(lat, lon, res); + if(first) { + first = 0; + *lon = -180; + if(fabs(viewpt) < .01) + *lat = 0; + else if(fabs(viewpt)<=1) + *lat = asin(viewpt)/RAD; + else + *lat = asin(1/viewpt)/RAD; + } else + *lon += res; + if(*lon <= 180) + return 1; + first = 1; + return -1; +} diff --git a/src/cmd/map/libmap/polyconic.c b/src/cmd/map/libmap/polyconic.c new file mode 100644 index 00000000..a07c97e3 --- /dev/null +++ b/src/cmd/map/libmap/polyconic.c @@ -0,0 +1,28 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +int +Xpolyconic(struct place *place, double *x, double *y) +{ + double r, alpha; + double lat2, lon2; + if(fabs(place->nlat.l) > .01) { + r = place->nlat.c / place->nlat.s; + alpha = place->wlon.l * place->nlat.s; + *y = place->nlat.l + r*(1 - cos(alpha)); + *x = - r*sin(alpha); + } else { + lon2 = place->wlon.l * place->wlon.l; + lat2 = place->nlat.l * place->nlat.l; + *y = place->nlat.l * (1+(lon2/2)*(1-(8+lon2)*lat2/12)); + *x = - place->wlon.l * (1-lat2*(3+lon2)/6); + } + return(1); +} + +proj +polyconic(void) +{ + return(Xpolyconic); +} diff --git a/src/cmd/map/libmap/rectangular.c b/src/cmd/map/libmap/rectangular.c new file mode 100644 index 00000000..d4a86c98 --- /dev/null +++ b/src/cmd/map/libmap/rectangular.c @@ -0,0 +1,22 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static double scale; + +static int +Xrectangular(struct place *place, double *x, double *y) +{ + *x = -scale*place->wlon.l; + *y = place->nlat.l; + return(1); +} + +proj +rectangular(double par) +{ + scale = cos(par*RAD); + if(scale<.1) + return 0; + return(Xrectangular); +} diff --git a/src/cmd/map/libmap/simpleconic.c b/src/cmd/map/libmap/simpleconic.c new file mode 100644 index 00000000..1ed1d1aa --- /dev/null +++ b/src/cmd/map/libmap/simpleconic.c @@ -0,0 +1,34 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static double r0, a; + +static int +Xsimpleconic(struct place *place, double *x, double *y) +{ + double r = r0 - place->nlat.l; + double t = a*place->wlon.l; + *x = -r*sin(t); + *y = -r*cos(t); + return 1; +} + +proj +simpleconic(double par0, double par1) +{ + struct coord lat0; + struct coord lat1; + deg2rad(par0,&lat0); + deg2rad(par1,&lat1); + if(fabs(lat0.l+lat1.l)<.01) + return rectangular(par0); + if(fabs(lat0.l-lat1.l)<.01) { + a = lat0.s/lat0.l; + r0 = lat0.c/lat0.s + lat0.l; + } else { + a = (lat1.c-lat0.c)/(lat0.l-lat1.l); + r0 = ((lat0.c+lat1.c)/a + lat1.l + lat0.l)/2; + } + return Xsimpleconic; +} diff --git a/src/cmd/map/libmap/sinusoidal.c b/src/cmd/map/libmap/sinusoidal.c new file mode 100644 index 00000000..6a79706e --- /dev/null +++ b/src/cmd/map/libmap/sinusoidal.c @@ -0,0 +1,17 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +int +Xsinusoidal(struct place *place, double *x, double *y) +{ + *x = - place->wlon.l * place->nlat.c; + *y = place->nlat.l; + return(1); +} + +proj +sinusoidal(void) +{ + return(Xsinusoidal); +} diff --git a/src/cmd/map/libmap/tetra.c b/src/cmd/map/libmap/tetra.c new file mode 100644 index 00000000..6bdef49b --- /dev/null +++ b/src/cmd/map/libmap/tetra.c @@ -0,0 +1,206 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +/* + * conformal map of earth onto tetrahedron + * the stages of mapping are + * (a) stereo projection of tetrahedral face onto + * isosceles curvilinear triangle with 3 120-degree + * angles and one straight side + * (b) map of this triangle onto half plane cut along + * 3 rays from the roots of unity to infinity + * formula (z^4+2*3^.5*z^2-1)/(z^4-2*3^.5*z^2-1) + * (c) do 3 times for each sector of plane: + * map of |arg z|<=pi/6, cut along z>1 into + * triangle |arg z|<=pi/6, Re z<=const, + * with upper side of cut going into upper half of + * of vertical side of triangle and lowere into lower + * formula int from 0 to z dz/sqrt(1-z^3) + * + * int from u to 1 3^.25*du/sqrt(1-u^3) = + F(acos((rt3-1+u)/(rt3+1-u)),sqrt(1/2+rt3/4)) + * int from 1 to u 3^.25*du/sqrt(u^3-1) = + * F(acos((rt3+1-u)/(rt3-1+u)),sqrt(1/2-rt3/4)) + * this latter formula extends analytically down to + * u=0 and is the basis of this routine, with the + * argument of complex elliptic integral elco2 + * being tan(acos...) + * the formula F(pi-x,k) = 2*F(pi/2,k)-F(x,k) is + * used to cross over into the region where Re(acos...)>pi/2 + * f0 and fpi are suitably scaled complete integrals +*/ + +#define TFUZZ 0.00001 + +static struct place tpole[4]; /* point of tangency of tetrahedron face*/ +static double tpoleinit[4][2] = { + 1., 0., + 1., 180., + -1., 90., + -1., -90. +}; +static struct tproj { + double tlat,tlon; /* center of stereo projection*/ + double ttwist; /* rotatn before stereo*/ + double trot; /*rotate after projection*/ + struct place projpl; /*same as tlat,tlon*/ + struct coord projtw; /*same as ttwist*/ + struct coord postrot; /*same as trot*/ +} tproj[4][4] = { +{/*00*/ {0.}, + /*01*/ {90., 0., 90., -90.}, + /*02*/ {0., 45., -45., 150.}, + /*03*/ {0., -45., -135., 30.} +}, +{/*10*/ {90., 0., -90., 90.}, + /*11*/ {0.}, + /*12*/ {0., 135., -135., -150.}, + /*13*/ {0., -135., -45., -30.} +}, +{/*20*/ {0., 45., 135., -30.}, + /*21*/ {0., 135., 45., -150.}, + /*22*/ {0.}, + /*23*/ {-90., 0., 180., 90.} +}, +{/*30*/ {0., -45., 45., -150.}, + /*31*/ {0., -135., 135., -30.}, + /*32*/ {-90., 0., 0., 90.}, + /*33*/ {0.} +}}; +static double tx[4] = { /*where to move facet after final rotation*/ + 0., 0., -1., 1. /*-1,1 to be sqrt(3)*/ +}; +static double ty[4] = { + 0., 2., -1., -1. +}; +static double root3; +static double rt3inv; +static double two_rt3; +static double tkc,tk,tcon; +static double f0r,f0i,fpir,fpii; + +static void +twhichp(struct place *g, int *p, int *q) +{ + int i,j,k; + double cosdist[4]; + struct place *tp; + for(i=0;i<4;i++) { + tp = &tpole[i]; + cosdist[i] = g->nlat.s*tp->nlat.s + + g->nlat.c*tp->nlat.c*( + g->wlon.s*tp->wlon.s + + g->wlon.c*tp->wlon.c); + } + j = 0; + for(i=1;i<4;i++) + if(cosdist[i] > cosdist[j]) + j = i; + *p = j; + k = j==0?1:0; + for(i=0;i<4;i++) + if(i!=j&&cosdist[i]>cosdist[k]) + k = i; + *q = k; +} + +int +Xtetra(struct place *place, double *x, double *y) +{ + int i,j; + struct place pl; + register struct tproj *tpp; + double vr, vi; + double br, bi; + double zr,zi,z2r,z2i,z4r,z4i,sr,si,tr,ti; + twhichp(place,&i,&j); + copyplace(place,&pl); + norm(&pl,&tproj[i][j].projpl,&tproj[i][j].projtw); + Xstereographic(&pl,&vr,&vi); + zr = vr/2; + zi = vi/2; + if(zr<=TFUZZ) + zr = TFUZZ; + csq(zr,zi,&z2r,&z2i); + csq(z2r,z2i,&z4r,&z4i); + z2r *= two_rt3; + z2i *= two_rt3; + cdiv(z4r+z2r-1,z4i+z2i,z4r-z2r-1,z4i-z2i,&sr,&si); + csqrt(sr-1,si,&tr,&ti); + cdiv(tcon*tr,tcon*ti,root3+1-sr,-si,&br,&bi); + if(br<0) { + br = -br; + bi = -bi; + if(!elco2(br,bi,tk,1.,1.,&vr,&vi)) + return 0; + vr = fpir - vr; + vi = fpii - vi; + } else + if(!elco2(br,bi,tk,1.,1.,&vr,&vi)) + return 0; + if(si>=0) { + tr = f0r - vi; + ti = f0i + vr; + } else { + tr = f0r + vi; + ti = f0i - vr; + } + tpp = &tproj[i][j]; + *x = tr*tpp->postrot.c + + ti*tpp->postrot.s + tx[i]; + *y = ti*tpp->postrot.c - + tr*tpp->postrot.s + ty[i]; + return(1); +} + +int +tetracut(struct place *g, struct place *og, double *cutlon) +{ + int i,j,k; + if((g->nlat.s<=-rt3inv&&og->nlat.s<=-rt3inv) && + (ckcut(g,og,*cutlon=0.)==2||ckcut(g,og,*cutlon=PI)==2)) + return(2); + twhichp(g,&i,&k); + twhichp(og,&j,&k); + if(i==j||i==0||j==0) + return(1); + return(0); +} + +proj +tetra(void) +{ + int i; + int j; + register struct place *tp; + register struct tproj *tpp; + double t; + root3 = sqrt(3.); + rt3inv = 1/root3; + two_rt3 = 2*root3; + tkc = sqrt(.5-.25*root3); + tk = sqrt(.5+.25*root3); + tcon = 2*sqrt(root3); + elco2(tcon/(root3-1),0.,tkc,1.,1.,&f0r,&f0i); + elco2(1.e15,0.,tk,1.,1.,&fpir,&fpii); + fpir *= 2; + fpii *= 2; + for(i=0;i<4;i++) { + tx[i] *= f0r*root3; + ty[i] *= f0r; + tp = &tpole[i]; + t = tp->nlat.s = tpoleinit[i][0]/root3; + tp->nlat.c = sqrt(1 - t*t); + tp->nlat.l = atan2(tp->nlat.s,tp->nlat.c); + deg2rad(tpoleinit[i][1],&tp->wlon); + for(j=0;j<4;j++) { + tpp = &tproj[i][j]; + latlon(tpp->tlat,tpp->tlon,&tpp->projpl); + deg2rad(tpp->ttwist,&tpp->projtw); + deg2rad(tpp->trot,&tpp->postrot); + } + } + return(Xtetra); +} + diff --git a/src/cmd/map/libmap/trapezoidal.c b/src/cmd/map/libmap/trapezoidal.c new file mode 100644 index 00000000..977cf60c --- /dev/null +++ b/src/cmd/map/libmap/trapezoidal.c @@ -0,0 +1,30 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static struct coord stdpar0, stdpar1; +static double k; +static double yeq; + +static int +Xtrapezoidal(struct place *place, double *x, double *y) +{ + *y = yeq + place->nlat.l; + *x = *y*k*place->wlon.l; + return 1; +} + +proj +trapezoidal(double par0, double par1) +{ + if(fabs(fabs(par0)-fabs(par1))<.1) + return rectangular(par0); + deg2rad(par0,&stdpar0); + deg2rad(par1,&stdpar1); + if(fabs(par1-par0) < .1) + k = stdpar1.s; + else + k = (stdpar1.c-stdpar0.c)/(stdpar0.l-stdpar1.l); + yeq = -stdpar1.l - stdpar1.c/k; + return Xtrapezoidal; +} diff --git a/src/cmd/map/libmap/twocirc.c b/src/cmd/map/libmap/twocirc.c new file mode 100644 index 00000000..0d0e48a4 --- /dev/null +++ b/src/cmd/map/libmap/twocirc.c @@ -0,0 +1,80 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +static double +quadratic(double a, double b, double c) +{ + double disc = b*b - 4*a*c; + return disc<0? 0: (-b-sqrt(disc))/(2*a); +} + +/* for projections with meridians being circles centered +on the x axis and parallels being circles centered on the y +axis. Find the intersection of the meridian thru (m,0), (90,0), +with the parallel thru (0,p), (p1,p2) */ + +static int +twocircles(double m, double p, double p1, double p2, double *x, double *y) +{ + double a; /* center of meridian circle, a>0 */ + double b; /* center of parallel circle, b>0 */ + double t,bb; + if(m > 0) { + twocircles(-m,p,p1,p2,x,y); + *x = -*x; + } else if(p < 0) { + twocircles(m,-p,p1,-p2,x,y); + *y = -*y; + } else if(p < .01) { + *x = m; + t = m/p1; + *y = p + (p2-p)*t*t; + } else if(m > -.01) { + *y = p; + *x = m - m*p*p; + } else { + b = p>=1? 1: p>.99? 0.5*(p+1 + p1*p1/(1-p)): + 0.5*(p*p-p1*p1-p2*p2)/(p-p2); + a = .5*(m - 1/m); + t = m*m-p*p+2*(b*p-a*m); + bb = b*b; + *x = quadratic(1+a*a/bb, -2*a + a*t/bb, + t*t/(4*bb) - m*m + 2*a*m); + *y = (*x*a+t/2)/b; + } + return 1; +} + +static int +Xglobular(struct place *place, double *x, double *y) +{ + twocircles(-2*place->wlon.l/PI, + 2*place->nlat.l/PI, place->nlat.c, place->nlat.s, x, y); + return 1; +} + +proj +globular(void) +{ + return Xglobular; +} + +static int +Xvandergrinten(struct place *place, double *x, double *y) +{ + double t = 2*place->nlat.l/PI; + double abst = fabs(t); + double pval = abst>=1? 1: abst/(1+sqrt(1-t*t)); + double p2 = 2*pval/(1+pval); + twocircles(-place->wlon.l/PI, pval, sqrt(1-p2*p2), p2, x, y); + if(t < 0) + *y = -*y; + return 1; +} + +proj +vandergrinten(void) +{ + return Xvandergrinten; +} diff --git a/src/cmd/map/libmap/zcoord.c b/src/cmd/map/libmap/zcoord.c new file mode 100644 index 00000000..7c3d3ad7 --- /dev/null +++ b/src/cmd/map/libmap/zcoord.c @@ -0,0 +1,143 @@ +#include <u.h> +#include <libc.h> +#include <stdio.h> +#include "map.h" + +static double cirmod(double); + +static struct place pole; /* map pole is tilted to here */ +static struct coord twist; /* then twisted this much */ +static struct place ipole; /* inverse transfrom */ +static struct coord itwist; + +void +orient(double lat, double lon, double theta) +{ + lat = cirmod(lat); + if(lat>90.) { + lat = 180. - lat; + lon -= 180.; + theta -= 180.; + } else if(lat < -90.) { + lat = -180. - lat; + lon -= 180.; + theta -= 180; + } + latlon(lat,lon,&pole); + deg2rad(theta, &twist); + latlon(lat,180.-theta,&ipole); + deg2rad(180.-lon, &itwist); +} + +void +latlon(double lat, double lon, struct place *p) +{ + lat = cirmod(lat); + if(lat>90.) { + lat = 180. - lat; + lon -= 180.; + } else if(lat < -90.) { + lat = -180. - lat; + lon -= 180.; + } + deg2rad(lat,&p->nlat); + deg2rad(lon,&p->wlon); +} + +void +deg2rad(double theta, struct coord *coord) +{ + theta = cirmod(theta); + coord->l = theta*RAD; + if(theta==90) { + coord->s = 1; + coord->c = 0; + } else if(theta== -90) { + coord->s = -1; + coord->c = 0; + } else + sincos(coord); +} + +static double +cirmod(double theta) +{ + while(theta >= 180.) + theta -= 360; + while(theta<-180.) + theta += 360.; + return(theta); +} + +void +sincos(struct coord *coord) +{ + coord->s = sin(coord->l); + coord->c = cos(coord->l); +} + +void +normalize(struct place *gg) +{ + norm(gg,&pole,&twist); +} + +void +invert(struct place *g) +{ + norm(g,&ipole,&itwist); +} + +void +norm(struct place *gg, struct place *pp, struct coord *tw) +{ + register struct place *g; /*geographic coords */ + register struct place *p; /* new pole in old coords*/ + struct place m; /* standard map coords*/ + g = gg; + p = pp; + if(p->nlat.s == 1.) { + if(p->wlon.l+tw->l == 0.) + return; + g->wlon.l -= p->wlon.l+tw->l; + } else { + if(p->wlon.l != 0) { + g->wlon.l -= p->wlon.l; + sincos(&g->wlon); + } + m.nlat.s = p->nlat.s * g->nlat.s + + p->nlat.c * g->nlat.c * g->wlon.c; + m.nlat.c = sqrt(1. - m.nlat.s * m.nlat.s); + m.nlat.l = atan2(m.nlat.s, m.nlat.c); + m.wlon.s = g->nlat.c * g->wlon.s; + m.wlon.c = p->nlat.c * g->nlat.s + - p->nlat.s * g->nlat.c * g->wlon.c; + m.wlon.l = atan2(m.wlon.s, - m.wlon.c) + - tw->l; + *g = m; + } + sincos(&g->wlon); + if(g->wlon.l>PI) + g->wlon.l -= 2*PI; + else if(g->wlon.l<-PI) + g->wlon.l += 2*PI; +} + +double +tan(double x) +{ + return(sin(x)/cos(x)); +} + +void +printp(struct place *g) +{ +printf("%.3f %.3f %.3f %.3f %.3f %.3f\n", +g->nlat.l,g->nlat.s,g->nlat.c,g->wlon.l,g->wlon.s,g->wlon.c); +} + +void +copyplace(struct place *g1, struct place *g2) +{ + *g2 = *g1; +} |