From 28994509cc11ac6a5443054dfae1fedfb69039bc Mon Sep 17 00:00:00 2001 From: rsc Date: Wed, 21 Apr 2004 22:19:33 +0000 Subject: Why not? --- src/cmd/map/index.c | 88 +++ src/cmd/map/iplot.h | 51 ++ src/cmd/map/libmap/aitoff.c | 26 + src/cmd/map/libmap/albers.c | 117 ++++ src/cmd/map/libmap/azequalarea.c | 19 + src/cmd/map/libmap/azequidist.c | 19 + src/cmd/map/libmap/bicentric.c | 25 + src/cmd/map/libmap/bonne.c | 36 ++ src/cmd/map/libmap/ccubrt.c | 13 + src/cmd/map/libmap/complex.c | 85 +++ src/cmd/map/libmap/conic.c | 27 + src/cmd/map/libmap/cubrt.c | 30 + src/cmd/map/libmap/cuts.c | 39 ++ src/cmd/map/libmap/cylequalarea.c | 24 + src/cmd/map/libmap/cylindrical.c | 19 + src/cmd/map/libmap/elco2.c | 132 ++++ src/cmd/map/libmap/elliptic.c | 35 ++ src/cmd/map/libmap/fisheye.c | 26 + src/cmd/map/libmap/gall.c | 29 + src/cmd/map/libmap/gilbert.c | 51 ++ src/cmd/map/libmap/guyou.c | 101 +++ src/cmd/map/libmap/harrison.c | 40 ++ src/cmd/map/libmap/hex.c | 122 ++++ src/cmd/map/libmap/homing.c | 121 ++++ src/cmd/map/libmap/lagrange.c | 30 + src/cmd/map/libmap/lambert.c | 46 ++ src/cmd/map/libmap/laue.c | 24 + src/cmd/map/libmap/lune.c | 62 ++ src/cmd/map/libmap/mercator.c | 36 ++ src/cmd/map/libmap/mkfile | 50 ++ src/cmd/map/libmap/mollweide.c | 25 + src/cmd/map/libmap/newyorker.c | 28 + src/cmd/map/libmap/orthographic.c | 35 ++ src/cmd/map/libmap/perspective.c | 84 +++ src/cmd/map/libmap/polyconic.c | 28 + src/cmd/map/libmap/rectangular.c | 22 + src/cmd/map/libmap/simpleconic.c | 34 + src/cmd/map/libmap/sinusoidal.c | 17 + src/cmd/map/libmap/tetra.c | 206 +++++++ src/cmd/map/libmap/trapezoidal.c | 30 + src/cmd/map/libmap/twocirc.c | 80 +++ src/cmd/map/libmap/zcoord.c | 143 +++++ src/cmd/map/map.c | 1226 +++++++++++++++++++++++++++++++++++++ src/cmd/map/map.h | 147 +++++ src/cmd/map/map.rc | 103 ++++ src/cmd/map/mapdemo.rc | 83 +++ src/cmd/map/mkfile | 32 + src/cmd/map/route.c | 131 ++++ src/cmd/map/sqrt.c | 52 ++ src/cmd/map/symbol.c | 192 ++++++ 50 files changed, 4221 insertions(+) create mode 100644 src/cmd/map/index.c create mode 100644 src/cmd/map/iplot.h create mode 100644 src/cmd/map/libmap/aitoff.c create mode 100644 src/cmd/map/libmap/albers.c create mode 100644 src/cmd/map/libmap/azequalarea.c create mode 100644 src/cmd/map/libmap/azequidist.c create mode 100644 src/cmd/map/libmap/bicentric.c create mode 100644 src/cmd/map/libmap/bonne.c create mode 100644 src/cmd/map/libmap/ccubrt.c create mode 100644 src/cmd/map/libmap/complex.c create mode 100644 src/cmd/map/libmap/conic.c create mode 100644 src/cmd/map/libmap/cubrt.c create mode 100644 src/cmd/map/libmap/cuts.c create mode 100644 src/cmd/map/libmap/cylequalarea.c create mode 100644 src/cmd/map/libmap/cylindrical.c create mode 100644 src/cmd/map/libmap/elco2.c create mode 100644 src/cmd/map/libmap/elliptic.c create mode 100644 src/cmd/map/libmap/fisheye.c create mode 100644 src/cmd/map/libmap/gall.c create mode 100644 src/cmd/map/libmap/gilbert.c create mode 100644 src/cmd/map/libmap/guyou.c create mode 100644 src/cmd/map/libmap/harrison.c create mode 100644 src/cmd/map/libmap/hex.c create mode 100644 src/cmd/map/libmap/homing.c create mode 100644 src/cmd/map/libmap/lagrange.c create mode 100644 src/cmd/map/libmap/lambert.c create mode 100644 src/cmd/map/libmap/laue.c create mode 100644 src/cmd/map/libmap/lune.c create mode 100644 src/cmd/map/libmap/mercator.c create mode 100644 src/cmd/map/libmap/mkfile create mode 100644 src/cmd/map/libmap/mollweide.c create mode 100644 src/cmd/map/libmap/newyorker.c create mode 100644 src/cmd/map/libmap/orthographic.c create mode 100644 src/cmd/map/libmap/perspective.c create mode 100644 src/cmd/map/libmap/polyconic.c create mode 100644 src/cmd/map/libmap/rectangular.c create mode 100644 src/cmd/map/libmap/simpleconic.c create mode 100644 src/cmd/map/libmap/sinusoidal.c create mode 100644 src/cmd/map/libmap/tetra.c create mode 100644 src/cmd/map/libmap/trapezoidal.c create mode 100644 src/cmd/map/libmap/twocirc.c create mode 100644 src/cmd/map/libmap/zcoord.c create mode 100644 src/cmd/map/map.c create mode 100644 src/cmd/map/map.h create mode 100755 src/cmd/map/map.rc create mode 100755 src/cmd/map/mapdemo.rc create mode 100644 src/cmd/map/mkfile create mode 100644 src/cmd/map/route.c create mode 100644 src/cmd/map/sqrt.c create mode 100644 src/cmd/map/symbol.c (limited to 'src/cmd/map') diff --git a/src/cmd/map/index.c b/src/cmd/map/index.c new file mode 100644 index 00000000..0ff75f5d --- /dev/null +++ b/src/cmd/map/index.c @@ -0,0 +1,88 @@ +#include +#include +#include "map.h" + +static proj Yaitoff(double p0, double p1){USED(p0); USED(p1); return aitoff();} +static proj Yalbers(double p0,double p1){USED(p0); USED(p1); return albers(p0,p1);} +static proj Yazequalarea(double p0, double p1){USED(p0); USED(p1); return azequalarea();} +static proj Yazequidistant(double p0, double p1){USED(p0); USED(p1); return azequidistant();} +static proj Ybicentric(double p0,double p1){USED(p0); USED(p1); return bicentric(p0);} +static proj Ybonne(double p0,double p1){USED(p0); USED(p1); return bonne(p0);} +static proj Yconic(double p0,double p1){USED(p0); USED(p1); return conic(p0);} +static proj Ycylequalarea(double p0,double p1){USED(p0); USED(p1); return cylequalarea(p0);} +static proj Ycylindrical(double p0, double p1){USED(p0); USED(p1); return cylindrical();} +static proj Yelliptic(double p0,double p1){USED(p0); USED(p1); return elliptic(p0);} +static proj Yfisheye(double p0,double p1){USED(p0); USED(p1); return fisheye(p0);} +static proj Ygall(double p0,double p1){USED(p0); USED(p1); return gall(p0);} +static proj Ygilbert(double p0, double p1){USED(p0); USED(p1); return gilbert();} +static proj Yglobular(double p0, double p1){USED(p0); USED(p1); return globular();} +static proj Ygnomonic(double p0, double p1){USED(p0); USED(p1); return gnomonic();} +static proj Yguyou(double p0, double p1){USED(p0); USED(p1); return guyou();} +static proj Yharrison(double p0,double p1){USED(p0); USED(p1); return harrison(p0,p1);} +static proj Yhex(double p0, double p1){USED(p0); USED(p1); return hex();} +static proj Yhoming(double p0,double p1){USED(p0); USED(p1); return homing(p0);} +static proj Ylagrange(double p0, double p1){USED(p0); USED(p1); return lagrange();} +static proj Ylambert(double p0,double p1){USED(p0); USED(p1); return lambert(p0,p1);} +static proj Ylaue(double p0, double p1){USED(p0); USED(p1); return laue();} +static proj Ylune(double p0,double p1){USED(p0); USED(p1); return lune(p0,p1);} +static proj Ymecca(double p0, double p1){USED(p0); USED(p1); return mecca(p0);} +static proj Ymercator(double p0, double p1){USED(p0); USED(p1); return mercator();} +static proj Ymollweide(double p0, double p1){USED(p0); USED(p1); return mollweide();} +static proj Ynewyorker(double p0,double p1){USED(p0); USED(p1); return newyorker(p0);} +static proj Yorthographic(double p0, double p1){USED(p0); USED(p1); return orthographic();} +static proj Yperspective(double p0,double p1){USED(p0); USED(p1); return perspective(p0);} +static proj Ypolyconic(double p0, double p1){USED(p0); USED(p1); return polyconic();} +static proj Yrectangular(double p0,double p1){USED(p0); USED(p1); return rectangular(p0);} +static proj Ysimpleconic(double p0,double p1){USED(p0); USED(p1); return simpleconic(p0,p1);} +static proj Ysinusoidal(double p0, double p1){USED(p0); USED(p1); return sinusoidal();} +static proj Ysp_albers(double p0,double p1){USED(p0); USED(p1); return sp_albers(p0,p1);} +static proj Ysp_mercator(double p0, double p1){USED(p0); USED(p1); return sp_mercator();} +static proj Ysquare(double p0, double p1){USED(p0); USED(p1); return square();} +static proj Ystereographic(double p0, double p1){USED(p0); USED(p1); return stereographic();} +static proj Ytetra(double p0, double p1){USED(p0); USED(p1); return tetra();} +static proj Ytrapezoidal(double p0,double p1){USED(p0); USED(p1); return trapezoidal(p0,p1);} +static proj Yvandergrinten(double p0, double p1){USED(p0); USED(p1); return vandergrinten();} + +struct index index[] = { + {"aitoff", Yaitoff, 0, picut, 0, 0, 0}, + {"albers", Yalbers, 2, picut, 3, 0, 0}, + {"azequalarea", Yazequalarea, 0, nocut, 1, 0, 0}, + {"azequidistant", Yazequidistant, 0, nocut, 1, 0, 0}, + {"bicentric", Ybicentric, 1, nocut, 0, 0, 0}, + {"bonne", Ybonne, 1, picut, 0, 0, 0}, + {"conic", Yconic, 1, picut, 0, 0, 0}, + {"cylequalarea", Ycylequalarea, 1, picut, 3, 0, 0}, + {"cylindrical", Ycylindrical, 0, picut, 0, 0, 0}, + {"elliptic", Yelliptic, 1, picut, 0, 0, 0}, + {"fisheye", Yfisheye, 1, nocut, 0, 0, 0}, + {"gall", Ygall, 1, picut, 3, 0, 0}, + {"gilbert", Ygilbert, 0, picut, 0, 0, 0}, + {"globular", Yglobular, 0, picut, 0, 0, 0}, + {"gnomonic", Ygnomonic, 0, nocut, 0, 0, plimb}, + {"guyou", Yguyou, 0, guycut, 0, 0, 0}, + {"harrison", Yharrison, 2, nocut, 0, 0, plimb}, + {"hex", Yhex, 0, hexcut, 0, 0, 0}, + {"homing", Yhoming, 1, nocut, 3, 0, hlimb}, + {"lagrange", Ylagrange,0,picut,0, 0, 0}, + {"lambert", Ylambert, 2, picut, 0, 0, 0}, + {"laue", Ylaue, 0, nocut, 0, 0, 0}, + {"lune", Ylune, 2, nocut, 0, 0, 0}, + {"mecca", Ymecca, 1, picut, 3, 0, mlimb}, + {"mercator", Ymercator, 0, picut, 3, 0, 0}, + {"mollweide", Ymollweide, 0, picut, 0, 0, 0}, + {"newyorker", Ynewyorker, 1, nocut, 0, 0, 0}, + {"orthographic", Yorthographic, 0, nocut, 0, 0, olimb}, + {"perspective", Yperspective, 1, nocut, 0, 0, plimb}, + {"polyconic", Ypolyconic, 0, picut, 0, 0, 0}, + {"rectangular", Yrectangular, 1, picut, 3, 0, 0}, + {"simpleconic", Ysimpleconic, 2, picut, 3, 0, 0}, + {"sinusoidal", Ysinusoidal, 0, picut, 0, 0, 0}, + {"sp_albers", Ysp_albers, 2, picut, 3, 1, 0}, + {"sp_mercator", Ysp_mercator, 0, picut, 0, 1, 0}, + {"square", Ysquare, 0, picut, 0, 0, 0}, + {"stereographic", Ystereographic, 0, nocut, 0, 0, 0}, + {"tetra", Ytetra, 0, tetracut, 0, 0, 0}, + {"trapezoidal", Ytrapezoidal, 2, picut, 3, 0, 0}, + {"vandergrinten", Yvandergrinten, 0, picut, 0, 0, 0}, + 0 +}; diff --git a/src/cmd/map/iplot.h b/src/cmd/map/iplot.h new file mode 100644 index 00000000..be472fff --- /dev/null +++ b/src/cmd/map/iplot.h @@ -0,0 +1,51 @@ +/* Plotting functions for v8 and v9 systems */ +/* This file is an alternative to plot.h */ + +/* open the plotting output */ +#define openpl() print("o\n") + +/* close the plotting output */ +#define closepl() print("cl\n") + +/* make sure the page or screen is clear */ +#define erase() print("e\n") + +/* plot a point at _x,_y, which becomes current */ +#define point(_x,_y) print("poi %d %d\n", _x,_y) + +/* coordinates to be assigned to lower left and upper right + corners of (square) plotting area */ +#define range(_x,_y,_X,_Y) print("ra %d %d %d %d\n", _x,_y,_X,_Y) + +/* place text, first letter at current point, which does not change */ +#define text(_s) {if(*(_s) == ' ')print("t \"%s\"\n",_s); else print("t %s\n", _s); } + +/* draw line from current point to _x,_y, which becomes current */ +#define vec(_x,_y) print("v %d %d\n", _x,_y) + +/* _x,_y becomes current point */ +#define move(_x, _y) print("m %d %d\n", _x, _y) + +/* specify style for drawing lines */ + +#define SOLID "solid" +#define DOTTED "dotted" +#define DASHED "dashed" +#define DOTDASH "dotdash" + +#define pen(_s) print("pe %s\n", _s) + +#define BLACK "z" +#define RED "r" +#define YELLOW "y" +#define GREEN "g" +#define BLUE "b" +#define CYAN "c" +#define MAGENTA "m" +#define WHITE "w" + +#define colorcode(_s) ((strcmp(_s,"black")==0)?BLACK:_s) + +#define colorx(_s) print("co %s\n", _s); /* funny name is all ken's fault */ + +#define comment(s,f) 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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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)=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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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.lPI-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 +#include +#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 +#include +#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.lplace->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 +#include +#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 +#include +#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 +#include +#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 +#include +#include "map.h" + + +static int +Xlaue(struct place *place, double *x, double *y) +{ + double r; + if(place->nlat.lnlat.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 +#include +#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 +#include +#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 +#include +#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 +#include +#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(rwlon.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 +#include +#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 +#include +#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= 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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#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 +#include +#include +#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; +} diff --git a/src/cmd/map/map.c b/src/cmd/map/map.c new file mode 100644 index 00000000..47109d43 --- /dev/null +++ b/src/cmd/map/map.c @@ -0,0 +1,1226 @@ +#include +#include +#include +#include "map.h" +#include "iplot.h" + +#define NVERT 20 /* max number of vertices in a -v polygon */ +#define HALFWIDTH 8192 /* output scaled to fit in -HALFWIDTH,HALFWIDTH */ +#define LONGLINES (HALFWIDTH*4) /* permissible segment lengths */ +#define SHORTLINES (HALFWIDTH/8) +#define SCALERATIO 10 /* of abs to rel data (see map(5)) */ +#define RESOL 2. /* coarsest resolution for tracing grid (degrees) */ +#define TWO_THRD 0.66666666666666667 + +int normproj(double, double, double *, double *); +int posproj(double, double, double *, double *); +int picut(struct place *, struct place *, double *); +double reduce(double); +short getshort(FILE *); +char *mapindex(char *); +proj projection; + + +static char *mapdir = "/lib/map"; /* default map directory */ +struct file { + char *name; + char *color; + char *style; +}; +static struct file dfltfile = { + "world", BLACK, SOLID /* default map */ +}; +static struct file *file = &dfltfile; /* list of map files */ +static int nfile = 1; /* length of list */ +static char *currcolor = BLACK; /* current color */ +static char *gridcolor = BLACK; +static char *bordcolor = BLACK; + +extern struct index index[]; +int halfwidth = HALFWIDTH; + +static int (*cut)(struct place *, struct place *, double *); +static int (*limb)(double*, double*, double); +static void dolimb(void); +static int onlimb; +static int poles; +static double orientation[3] = { 90., 0., 0. }; /* -o option */ +static int oriented; /* nonzero if -o option occurred */ +static int upright; /* 1 if orientation[0]==90, -1 if -90, else 0*/ +static int delta = 1; /* -d setting */ +static double limits[4] = { /* -l parameters */ + -90., 90., -180., 180. +}; +static double klimits[4] = { /* -k parameters */ + -90., 90., -180., 180. +}; +static int limcase; +static double rlimits[4]; /* limits expressed in radians */ +static double lolat, hilat, lolon, hilon; +static double window[4] = { /* option -w */ + -90., 90., -180., 180. +}; +static int windowed; /* nozero if option -w */ +static struct vert { double x, y; } v[NVERT+2]; /*clipping polygon*/ +static struct edge { double a, b, c; } e[NVERT]; /* coeffs for linear inequality */ +static int nvert; /* number of vertices in clipping polygon */ + +static double rwindow[4]; /* window, expressed in radians */ +static double params[2]; /* projection params */ +/* bounds on output values before scaling; found by coarse survey */ +static double xmin = 100.; +static double xmax = -100.; +static double ymin = 100.; +static double ymax = -100.; +static double xcent, ycent; +static double xoff, yoff; +double xrange, yrange; +static int left = -HALFWIDTH; +static int right = HALFWIDTH; +static int bottom = -HALFWIDTH; +static int top = HALFWIDTH; +static int longlines = SHORTLINES; /* drop longer segments */ +static int shortlines = SHORTLINES; +static int bflag = 1; /* 0 for option -b */ +static int s1flag = 0; /* 1 for option -s1 */ +static int s2flag = 0; /* 1 for option -s2 */ +static int rflag = 0; /* 1 for option -r */ +static int kflag = 0; /* 1 if option -k occurred */ +static int xflag = 0; /* 1 for option -x */ + int vflag = 1; /* -1 if option -v occurred */ +static double position[3]; /* option -p */ +static double center[3] = {0., 0., 0.}; /* option -c */ +static struct coord crot; /* option -c */ +static double grid[3] = { 10., 10., RESOL }; /* option -g */ +static double dlat, dlon; /* resolution for tracing grid in lat and lon */ +static double scaling; /* to compute final integer output */ +static struct file *track; /* options -t and -u */ +static int ntrack; /* number of tracks present */ +static char *symbolfile; /* option -y */ + +void clamp(double *px, double v); +void clipinit(void); +double diddle(struct place *, double, double); +double diddle(struct place *, double, double); +void dobounds(double, double, double, double, int); +void dogrid(double, double, double, double); +int duple(struct place *, double); +double fmax(double, double); +double fmin(double, double); +void getdata(char *); +int gridpt(double, double, int); +int inpoly(double, double); +int inwindow(struct place *); +void pathnames(void); +int pnorm(double); +void radbds(double *w, double *rw); +void revlon(struct place *, double); +void satellite(struct file *); +int seeable(double, double); +void windlim(void); +void realcut(void); + +int +option(char *s) +{ + + if(s[0]=='-' && (s[1]<'0'||s[1]>'9')) + return(s[1]!='.'&&s[1]!=0); + else + return(0); +} + +void +conv(int k, struct coord *g) +{ + g->l = (0.0001/SCALERATIO)*k; + sincos(g); +} + +int +main(int argc, char *argv[]) +{ + int i,k; + char *s, *t, *style; + double x, y; + double lat, lon; + double *wlim; + double dd; + if(sizeof(short)!=2) + abort(); /* getshort() won't work */ + s = getenv("MAP"); + if(s) + file[0].name = s; + s = getenv("MAPDIR"); + if(s) + mapdir = s; + if(argc<=1) + error("usage: map projection params options"); + for(k=0;index[k].name;k++) { + s = index[k].name; + t = argv[1]; + while(*s == *t){ + if(*s==0) goto found; + s++; + t++; + } + } + fprintf(stderr,"projections:\n"); + for(i=0;index[i].name;i++) { + fprintf(stderr,"%s",index[i].name); + for(k=0; k=argc||option(argv[i])) { + fprintf(stderr,"%s needs %d params\n",index[k].name,index[k].npar); + exits("error"); + } + params[i] = atof(argv[i]); + } + argv += i; + argc -= i; + while(argc>0&&option(argv[0])) { + argc--; + argv++; + switch(argv[-1][1]) { + case 'm': + if(file == &dfltfile) { + file = 0; + nfile = 0; + } + while(argc && !option(*argv)) { + file = realloc(file,(nfile+1)*sizeof(*file)); + file[nfile].name = *argv; + file[nfile].color = currcolor; + file[nfile].style = SOLID; + nfile++; + argv++; + argc--; + } + break; + case 'b': + bflag = 0; + for(nvert=0;nvert=2;nvert++) { + if(option(*argv)) + break; + v[nvert].x = atof(*argv++); + argc--; + if(option(*argv)) + break; + v[nvert].y = atof(*argv++); + argc--; + } + if(nvert>=NVERT) + error("too many clipping vertices"); + break; + case 'g': + gridcolor = currcolor; + for(i=0;i<3&&argc>i&&!option(argv[i]);i++) + grid[i] = atof(argv[i]); + switch(i) { + case 0: + grid[0] = grid[1] = 0.; + break; + case 1: + grid[1] = grid[0]; + } + argc -= i; + argv += i; + break; + case 't': + style = SOLID; + goto casetu; + case 'u': + style = DOTDASH; + casetu: + while(argc && !option(*argv)) { + track = realloc(track,(ntrack+1)*sizeof(*track)); + track[ntrack].name = *argv; + track[ntrack].color = currcolor; + track[ntrack].style = style; + ntrack++; + argv++; + argc--; + } + break; + case 'r': + rflag++; + break; + case 's': + switch(argv[-1][2]) { + case '1': + s1flag++; + break; + case 0: /* compatibility */ + case '2': + s2flag++; + } + break; + case 'o': + for(i=0;i<3&&i0&&!option(argv[0])) { + delta = atoi(argv[0]); + argv++; + argc--; + } + break; + case 'w': + bordcolor = currcolor; + windowed++; + for(i=0;ii&&!option(argv[i]);i++) + center[i] = atof(argv[i]); + argc -= i; + argv += i; + break; + case 'p': + for(i=0;i<3&&argc>i&&!option(argv[i]);i++) + position[i] = atof(argv[i]); + argc -= i; + argv += i; + if(i!=3||position[2]<=0) + error("incomplete positioning"); + break; + case 'y': + if(argc>0&&!option(argv[0])) { + symbolfile = argv[0]; + argc--; + argv++; + } + break; + case 'v': + if(index[k].limb == 0) + error("-v does not apply here"); + vflag = -1; + break; + case 'x': + xflag = 1; + break; + case 'C': + if(argc && !option(*argv)) { + currcolor = colorcode(*argv); + argc--; + argv++; + } + break; + } + } + if(argc>0) + error("error in arguments"); + pathnames(); + clamp(&limits[0],-90.); + clamp(&limits[1],90.); + clamp(&klimits[0],-90.); + clamp(&klimits[1],90.); + clamp(&window[0],-90.); + clamp(&window[1],90.); + radbds(limits,rlimits); + limcase = limits[2]<-180.?0: + limits[3]>180.?2: + 1; + if( + window[0]>=window[1]|| + window[2]>=window[3]|| + window[0]>90.|| + window[1]<-90.|| + window[2]>180.|| + window[3]<-180.) + error("unreasonable window"); + windlim(); + radbds(window,rwindow); + upright = orientation[0]==90? 1: orientation[0]==-90? -1: 0; + if(index[k].spheroid && !upright) + error("can't tilt the spheroid"); + if(limits[2]>limits[3]) + limits[3] += 360; + if(!oriented) + orientation[2] = (limits[2]+limits[3])/2; + orient(orientation[0],orientation[1],orientation[2]); + projection = (*index[k].prog)(params[0],params[1]); + if(projection == 0) + error("unreasonable projection parameters"); + clipinit(); + grid[0] = fabs(grid[0]); + grid[1] = fabs(grid[1]); + if(!kflag) + for(i=0;i<4;i++) + klimits[i] = limits[i]; + if(klimits[2]>klimits[3]) + klimits[3] += 360; + lolat = limits[0]; + hilat = limits[1]; + lolon = limits[2]; + hilon = limits[3]; + if(lolon>=hilon||lolat>=hilat||lolat<-90.||hilat>90.) + error("unreasonable limits"); + wlim = kflag? klimits: window; + dlat = fmin(hilat-lolat,wlim[1]-wlim[0])/16; + dlon = fmin(hilon-lolon,wlim[3]-wlim[2])/32; + dd = fmax(dlat,dlon); + while(grid[2]>fmin(dlat,dlon)/2) + grid[2] /= 2; + realcut(); + if(nvert<=0) { + for(lat=klimits[0];latklimits[1]) + lat = klimits[1]; + for(lon=klimits[2];lonxmax) xmax = x; + if(yymax) ymax = y; + } + } + } else { + for(i=0; ixmax) xmax = x; + if(yymax) ymax = y; + } + } + xrange = xmax - xmin; + yrange = ymax - ymin; + if(xrange<=0||yrange<=0) + error("map seems to be empty"); + scaling = 2; /*plotting area from -1 to 1*/ + if(position[2]!=0) { + if(posproj(position[0]-.5,position[1],&xcent,&ycent)==0|| + posproj(position[0]+.5,position[1],&x,&y)==0) + error("unreasonable position"); + scaling /= (position[2]*hypot(x-xcent,y-ycent)); + if(posproj(position[0],position[1],&xcent,&ycent)==0) + error("unreasonable position"); + } else { + scaling /= (xrange>yrange?xrange:yrange); + xcent = (xmin+xmax)/2; + ycent = (ymin+ymax)/2; + } + xoff = center[0]/scaling; + yoff = center[1]/scaling; + crot.l = center[2]*RAD; + sincos(&crot); + scaling *= HALFWIDTH*0.9; + if(symbolfile) + getsyms(symbolfile); + if(!s2flag) { + openpl(); + erase(); + } + range(left,bottom,right,top); + comment("grid",""); + colorx(gridcolor); + pen(DOTTED); + if(grid[0]>0.) + for(lat=ceil(lolat/grid[0])*grid[0]; + lat<=hilat;lat+=grid[0]) + dogrid(lat,lat,lolon,hilon); + if(grid[1]>0.) + for(lon=ceil(lolon/grid[1])*grid[1]; + lon<=hilon;lon+=grid[1]) + dogrid(lolat,hilat,lon,lon); + comment("border",""); + colorx(bordcolor); + pen(SOLID); + if(bflag) { + dolimb(); + dobounds(lolat,hilat,lolon,hilon,0); + dobounds(window[0],window[1],window[2],window[3],1); + } + lolat = floor(limits[0]/10)*10; + hilat = ceil(limits[1]/10)*10; + lolon = floor(limits[2]/10)*10; + hilon = ceil(limits[3]/10)*10; + if(lolon>hilon) + hilon += 360.; + /*do tracks first so as not to lose the standard input*/ + for(i=0;inlat.lnlat.l>rwindow[1]+FUZZ|| + geog->wlon.lwlon.l>rwindow[3]+FUZZ) + return(0); + else return(1); +} + +int +inlimits(struct place *g) +{ + if(rlimits[0]-FUZZ>g->nlat.l|| + rlimits[1]+FUZZnlat.l) + return(0); + switch(limcase) { + case 0: + if(rlimits[2]+TWOPI-FUZZ>g->wlon.l&& + rlimits[3]+FUZZwlon.l) + return(0); + break; + case 1: + if(rlimits[2]-FUZZ>g->wlon.l|| + rlimits[3]+FUZZwlon.l) + return(0); + break; + case 2: + if(rlimits[2]>g->wlon.l&& + rlimits[3]-TWOPI+FUZZwlon.l) + return(0); + break; + } + return(1); +} + + +long patch[18][36]; + +void +getdata(char *mapfile) +{ + char *indexfile; + int kx,ky,c; + int k; + long b; + long *p; + int ip, jp; + int n; + struct place g; + int i, j; + double lat, lon; + int conn; + FILE *ifile, *xfile; + + indexfile = mapindex(mapfile); + xfile = fopen(indexfile,"r"); + if(xfile==NULL) + filerror("can't find map index", indexfile); + free(indexfile); + for(i=0,p=patch[0];i<18*36;i++,p++) + *p = 1; + while(!feof(xfile) && fscanf(xfile,"%d%d%ld",&i,&j,&b)==3) + patch[i+9][j+18] = b; + fclose(xfile); + ifile = fopen(mapfile,"r"); + if(ifile==NULL) + filerror("can't find map data", mapfile); + for(lat=lolat;lat=0&&(jp=getc(ifile))>=0){ + if(ip!=(i&0377)||jp!=(j&0377)) + break; + n = getshort(ifile); + conn = 0; + if(n > 0) { /* absolute coordinates */ + kx = ky = 0; /* set */ + for(k=0;k0) + return(1); + return(0); +} + +void +satellite(struct file *t) +{ + char sym[50]; + char lbl[50]; + double scale; + int conn; + double lat,lon; + struct place place; + static FILE *ifile; + + if(ifile == nil) + ifile = stdin; + if(t->name[0]!='-'||t->name[1]!=0) { + fclose(ifile); + if((ifile=fopen(t->name,"r"))==NULL) + filerror("can't find track", t->name); + } + comment("track",t->name); + colorx(t->color); + pen(t->style); + for(;;) { + conn = 0; + while(!feof(ifile) && fscanf(ifile,"%lf%lf",&lat,&lon)==2){ + latlon(lat,lon,&place); + if(fscanf(ifile,"%1s",lbl) == 1) { + if(strchr("+-.0123456789",*lbl)==0) + break; + ungetc(*lbl,ifile); + } + conn = plotpt(&place,conn); + } + if(feof(ifile)) + return; + fscanf(ifile,"%[^\n]",lbl+1); + switch(*lbl) { + case '"': + if(plotpt(&place,conn)) + text(lbl+1); + break; + case ':': + case '!': + if(sscanf(lbl+1,"%s %lf",sym,&scale) <= 1) + scale = 1; + if(plotpt(&place,conn?conn:-1)) { + int r = *lbl=='!'?0:rflag?-1:1; + pen(SOLID); + if(putsym(&place,sym,scale,r) == 0) + text(lbl); + pen(t->style); + } + break; + default: + if(plotpt(&place,conn)) + text(lbl); + break; + } + } +} + +int +pnorm(double x) +{ + int i; + i = x/10.; + i %= 36; + if(i>=18) return(i-36); + if(i<-18) return(i+36); + return(i); +} + +void +error(char *s) +{ + fprintf(stderr,"map: \r\n%s\n",s); + exits("error"); +} + +void +filerror(char *s, char *f) +{ + fprintf(stderr,"\r\n%s %s\n",s,f); + exits("error"); +} + +char * +mapindex(char *s) +{ + char *t = malloc(strlen(s)+3); + strcpy(t,s); + strcat(t,".x"); + return t; +} + +#define NOPT 32767 +static int ox = NOPT; +static int oy = NOPT; + +int +cpoint(int xi, int yi, int conn) +{ + int dx = abs(ox-xi); + int dy = abs(oy-yi); + if(!xflag && (xi=right || yi=top)) { + ox = oy = NOPT; + return 0; + } + if(conn == -1) /* isolated plotting symbol */ + {} + else if(!conn) + point(xi,yi); + else { + if(dx+dy>longlines) { + ox = oy = NOPT; /* don't leap across cuts */ + return 0; + } + if(dx || dy) + vec(xi,yi); + } + ox = xi, oy = yi; + return dx+dy<=2? 2: 1; /* 2=very near; see dogrid */ +} + + +struct place oldg; + +int +plotpt(struct place *g, int conn) +{ + int kx,ky; + int ret; + double cutlon; + if(!inlimits(g)) { + return(0); +} + normalize(g); + if(!inwindow(g)) { + return(0); +} + switch((*cut)(g,&oldg,&cutlon)) { + case 2: + if(conn) { + ret = duple(g,cutlon)|duple(g,cutlon); + oldg = *g; + return(ret); + } + case 0: + conn = 0; + default: /* prevent diags about bad return value */ + case 1: + oldg = *g; + ret = doproj(g,&kx,&ky); + if(ret==0 || !onlimb && ret*vflag<=0) + return(0); + ret = cpoint(kx,ky,conn); + return ret; + } +} + +int +doproj(struct place *g, int *kx, int *ky) +{ + int i; + double x,y,x1,y1; +/*fprintf(stderr,"dopr1 %f %f \n",g->nlat.l,g->wlon.l);*/ + i = fixproj(g,&x,&y); + if(i == 0) + return(0); + if(rflag) + x = -x; +/*fprintf(stderr,"dopr2 %f %f\n",x,y);*/ + if(!inpoly(x,y)) { + return 0; +} + x1 = x - xcent; + y1 = y - ycent; + x = (x1*crot.c - y1*crot.s + xoff)*scaling; + y = (x1*crot.s + y1*crot.c + yoff)*scaling; + *kx = x + (x>0?.5:-.5); + *ky = y + (y>0?.5:-.5); + return(i); +} + +int +duple(struct place *g, double cutlon) +{ + int kx,ky; + int okx,oky; + struct place ig; + revlon(g,cutlon); + revlon(&oldg,cutlon); + ig = *g; + invert(&ig); + if(!inlimits(&ig)) + return(0); + if(doproj(g,&kx,&ky)*vflag<=0 || + doproj(&oldg,&okx,&oky)*vflag<=0) + return(0); + cpoint(okx,oky,0); + cpoint(kx,ky,1); + return(1); +} + +void +revlon(struct place *g, double cutlon) +{ + g->wlon.l = reduce(cutlon-reduce(g->wlon.l-cutlon)); + sincos(&g->wlon); +} + + +/* recognize problems of cuts + * move a point across cut to side of its predecessor + * if its very close to the cut + * return(0) if cut interrupts the line + * return(1) if line is to be drawn normally + * return(2) if line is so close to cut as to + * be properly drawn on both sheets +*/ + +int +picut(struct place *g, struct place *og, double *cutlon) +{ + *cutlon = PI; + return(ckcut(g,og,PI)); +} + +int +nocut(struct place *g, struct place *og, double *cutlon) +{ + USED(g); + USED(og); + USED(cutlon); +/* +#pragma ref g +#pragma ref og +#pragma ref cutlon +*/ + return(1); +} + +int +ckcut(struct place *g1, struct place *g2, double lon) +{ + double d1, d2; + double f1, f2; + int kx,ky; + d1 = reduce(g1->wlon.l -lon); + d2 = reduce(g2->wlon.l -lon); + if((f1=fabs(d1))0) + cpoint(kx,ky,0); + } + if(f1PI*TWO_THRD||f2>PI*TWO_THRD) + return(1); + return(d1*d2>=0); +} + +double +diddle(struct place *g, double lon, double d) +{ + double d1; + d1 = FUZZ/2; + if(d<0) + d1 = -d1; + g->wlon.l = reduce(lon+d1); + sincos(&g->wlon); + return(d1); +} + +double +reduce(double lon) +{ + if(lon>PI) + lon -= 2*PI; + else if(lon<-PI) + lon += 2*PI; + return(lon); +} + + +double tetrapt = 35.26438968; /* atan(1/sqrt(2)) */ + +void +dogrid(double lat0, double lat1, double lon0, double lon1) +{ + double slat,slon,tlat,tlon; + register int conn, oconn; + slat = tlat = slon = tlon = 0; + if(lat1>lat0) + slat = tlat = fmin(grid[2],dlat); + else + slon = tlon = fmin(grid[2],dlon);; + conn = oconn = 0; + while(lat0<=lat1&&lon0<=lon1) { + conn = gridpt(lat0,lon0,conn); + if(projection==Xguyou&&slat>0) { + if(lat0<-45&&lat0+slat>-45) + conn = gridpt(-45.,lon0,conn); + else if(lat0<45&&lat0+slat>45) + conn = gridpt(45.,lon0,conn); + } else if(projection==Xtetra&&slat>0) { + if(lat0<-tetrapt&&lat0+slat>-tetrapt) { + gridpt(-tetrapt-.001,lon0,conn); + conn = gridpt(-tetrapt+.001,lon0,0); + } + else if(lat0tetrapt) { + gridpt(tetrapt-.001,lon0,conn); + conn = gridpt(tetrapt+.001,lon0,0); + } + } + if(conn==0 && oconn!=0) { + if(slat+slon>.05) { + lat0 -= slat; /* steps too big */ + lon0 -= slon; /* or near bdry */ + slat /= 2; + slon /= 2; + conn = oconn = gridpt(lat0,lon0,conn); + } else + oconn = 0; + } else { + if(conn==2) { + slat = tlat; + slon = tlon; + conn = 1; + } + oconn = conn; + } + lat0 += slat; + lon0 += slon; + } + gridpt(lat1,lon1,conn); +} + +static int gridinv; /* nonzero when doing window bounds */ + +int +gridpt(double lat, double lon, int conn) +{ + struct place g; +/*fprintf(stderr,"%f %f\n",lat,lon);*/ + latlon(lat,lon,&g); + if(gridinv) + invert(&g); + return(plotpt(&g,conn)); +} + +/* win=0 ordinary grid lines, win=1 window lines */ + +void +dobounds(double lolat, double hilat, double lolon, double hilon, int win) +{ + gridinv = win; + if(lolat>-90 || win && (poles&1)!=0) + dogrid(lolat+FUZZ,lolat+FUZZ,lolon,hilon); + if(hilat<90 || win && (poles&2)!=0) + dogrid(hilat-FUZZ,hilat-FUZZ,lolon,hilon); + if(hilon-lolon<360 || win && cut==picut) { + dogrid(lolat,hilat,lolon+FUZZ,lolon+FUZZ); + dogrid(lolat,hilat,hilon-FUZZ,hilon-FUZZ); + } + gridinv = 0; +} + +static void +dolimb(void) +{ + double lat, lon; + double res = fmin(dlat, dlon)/4; + int conn = 0; + int newconn; + if(limb == 0) + return; + onlimb = gridinv = 1; + for(;;) { + newconn = (*limb)(&lat, &lon, res); + if(newconn == -1) + break; + conn = gridpt(lat, lon, conn*newconn); + } + onlimb = gridinv = 0; +} + + +void +radbds(double *w, double *rw) +{ + int i; + for(i=0;i<4;i++) + rw[i] = w[i]*RAD; + rw[0] -= FUZZ; + rw[1] += FUZZ; + rw[2] -= FUZZ; + rw[3] += FUZZ; +} + +void +windlim(void) +{ + double center = orientation[0]; + double colat; + if(center>90) + center = 180 - center; + if(center<-90) + center = -180 - center; + if(fabs(center)>90) + error("unreasonable orientation"); + colat = 90 - window[0]; + if(center-colat>limits[0]) + limits[0] = center - colat; + if(center+colat 16 bits */ + return r; +} + +double +fmin(double x, double y) +{ + return(xy?x:y); +} + +void +clamp(double *px, double v) +{ + *px = (v<0?fmax:fmin)(*px,v); +} + +void +pathnames(void) +{ + int i; + char *t, *indexfile, *name; + FILE *f, *fx; + for(i=0; i=0) s = 1; + if(t>FUZZ && s<=0) s = -1; + if(-FUZZ<=t&&t<=FUZZ || t*s>0) { + s = 0; + break; + } + } + if(s==0) + error("improper clipping polygon"); + for(i=0; ia + y*ei->b - ei->c; + if(val>10*FUZZ) + return(0); + } + return 1; +} + +void +realcut() +{ + struct place g; + double lat; + + if(cut != picut) /* punt on unusual cuts */ + return; + for(lat=window[0]; lat<=window[1]; lat+=grid[2]) { + g.wlon.l = PI; + sincos(&g.wlon); + g.nlat.l = lat*RAD; + sincos(&g.nlat); + if(!inwindow(&g)) { + break; +} + invert(&g); + if(inlimits(&g)) { + return; +} + } + longlines = shortlines = LONGLINES; + cut = nocut; /* not necessary; small eff. gain */ +} diff --git a/src/cmd/map/map.h b/src/cmd/map/map.h new file mode 100644 index 00000000..d96497c8 --- /dev/null +++ b/src/cmd/map/map.h @@ -0,0 +1,147 @@ +/* +#pragma lib "/sys/src/cmd/map/libmap/libmap.a$O" +#pragma src "/sys/src/cmd/map/libmap" +*/ + +#define index index0 +#ifndef PI +#define PI 3.1415926535897932384626433832795028841971693993751 +#endif + +#define TWOPI (2*PI) +#define RAD (PI/180) +double hypot(double, double); /* sqrt(a*a+b*b) */ +double tan(double); /* not in K&R library */ + +#define ECC .08227185422 /* eccentricity of earth */ +#define EC2 .006768657997 + +#define FUZZ .0001 +#define UNUSED 0.0 /* a dummy double parameter */ + +struct coord { + double l; /* lat or lon in radians*/ + double s; /* sin */ + double c; /* cos */ +}; +struct place { + struct coord nlat; + struct coord wlon; +}; + +typedef int (*proj)(struct place *, double *, double *); + +struct index { /* index of known projections */ + char *name; /* name of projection */ + proj (*prog)(double, double); + /* pointer to projection function */ + int npar; /* number of params */ + int (*cut)(struct place *, struct place *, double *); + /* function that handles cuts--eg longitude 180 */ + int poles; /*1 S pole is a line, 2 N pole is, 3 both*/ + int spheroid; /* poles must be at 90 deg if nonzero */ + int (*limb)(double *lat, double *lon, double resolution); + /* get next place on limb */ + /* return -1 if done, 0 at gap, else 1 */ +}; + + +proj aitoff(void); +proj albers(double, double); +int Xazequalarea(struct place *, double *, double *); +proj azequalarea(void); +int Xazequidistant(struct place *, double *, double *); +proj azequidistant(void); +proj bicentric(double); +proj bonne(double); +proj conic(double); +proj cylequalarea(double); +int Xcylindrical(struct place *, double *, double *); +proj cylindrical(void); +proj elliptic(double); +proj fisheye(double); +proj gall(double); +proj gilbert(void); +proj globular(void); +proj gnomonic(void); +int guycut(struct place *, struct place *, double *); +int Xguyou(struct place *, double *, double *); +proj guyou(void); +proj harrison(double, double); +int hexcut(struct place *, struct place *, double *); +proj hex(void); +proj homing(double); +int hlimb(double*, double*, double resolution); +proj lagrange(void); +proj lambert(double, double); +proj laue(void); +proj lune(double, double); +proj loxodromic(double); /* not in library */ +proj mecca(double); +int mlimb(double*, double*, double resolution); +proj mercator(void); +proj mollweide(void); +proj newyorker(double); +proj ortelius(double, double); /* not in library */ +int Xorthographic(struct place *place, double *x, double *y); +proj orthographic(void); +int olimb(double*, double*, double); +proj perspective(double); +int plimb(double*, double*, double resolution); +int Xpolyconic(struct place *, double *, double *); +proj polyconic(void); +proj rectangular(double); +proj simpleconic(double, double); +int Xsinusoidal(struct place *, double *, double *); +proj sinusoidal(void); +proj sp_albers(double, double); +proj sp_mercator(void); +proj square(void); +int Xstereographic(struct place *, double *, double *); +proj stereographic(void); +int Xtetra(struct place *, double *, double *); +int tetracut(struct place *, struct place *, double *); +proj tetra(void); +proj trapezoidal(double, double); +proj vandergrinten(void); +proj wreath(double, double); /* not in library */ + +void findxy(double, double *, double *); +void albscale(double, double, double, double); +void invalb(double, double, double *, double *); + +void cdiv(double, double, double, double, double *, double *); +void cmul(double, double, double, double, double *, double *); +void cpow(double, double, double *, double *, double); +void csq(double, double, double *, double *); +void csqrt(double, double, double *, double *); +void ccubrt(double, double, double *, double *); +double cubrt(double); +int elco2(double, double, double, double, double, double *, double *); +void cdiv2(double, double, double, double, double *, double *); +void csqr(double, double, double *, double *); + +void orient(double, double, double); +void latlon(double, double, struct place *); +void deg2rad(double, struct coord *); +void sincos(struct coord *); +void normalize(struct place *); +void invert(struct place *); +void norm(struct place *, struct place *, struct coord *); +void printp(struct place *); +void copyplace(struct place *, struct place *); + +int picut(struct place *, struct place *, double *); +int ckcut(struct place *, struct place *, double); +double reduce(double); + +void getsyms(char *); +int putsym(struct place *, char *, double, int); +void filerror(char *s, char *f); +void error(char *s); +int doproj(struct place *, int *, int *); +int cpoint(int, int, int); +int plotpt(struct place *, int); +int nocut(struct place *, struct place *, double *); + +extern int (*projection)(struct place *, double *, double *); diff --git a/src/cmd/map/map.rc b/src/cmd/map/map.rc new file mode 100755 index 00000000..eaa75c91 --- /dev/null +++ b/src/cmd/map/map.rc @@ -0,0 +1,103 @@ +#!/bin/rc + +rfork en + +# F FEATUREs, M map files, A other arguments +FEATURE=no + +if (~ $MAPPROG '') + MAPPROG=/bin/aux/mapd + +if (~ $MAPDIR '') + MAPDIR=/lib/map + +F=(); M=(); A=(); +for (i) { + switch ($FEATURE) { + case no + switch ($i) { + case -f + FEATURE=yes + F=($F) + case * + A=($A $i) + } + case yes + switch ($i) { + case -f + case -* + A=($A $i) + FEATURE=no + case riv*2 + F=($F 201 202) + case riv*3 + F=($F 201 202 203) + case riv*4 + F=($F 201 202 203 204) + case riv* + F=($F 201) + case iriv*2 + F=($F 206 207) + case iriv*[34] + F=($F 206 207 208) + case iriv* + F=($F 206) + case coast*2 shore*2 lake*2 + F=($F 102) + case coast*3 shore*3 lake*3 + F=($F 102 103) + case coast*4 shore*4 lake*4 + F=($F 102 103 104) + case coast* shore* lake* + case ilake*[234] ishore*[234] + F=($F 106 107) + case ilake* ishore* + F=($F 106) + case reef* + F=($F 108) + case canal*2 + F=($F 210 211) + case canal*[34] + F=($F 210 211 212) + case canal* + F=($F 210) + case glacier* + F=($F 115) + case state* province* + F=($F 401) + case countr*2 + F=($F 301 302) + case countr*[34] + F=($F 301 302 303) + case countr* + F=($F 301) + case salt*[234] + F=($F 109 110) + case salt* + F=($F 109) + case ice*[234] shel*[234] + F=($F 113 114) + case ice* shel* + F=($F 113) + case * + echo map: unknown feature $i >[1=2] + exits "unknown feature" + } + } +} + +for (j in $F) { + if (test -r $MAPDIR/$j) + M=($M $MAPDIR/$j) +} + +if (~ $F ?*) { + if (test -r $MAPDIR/101) + M=(101 $M) + M=(-m $M) +} + +if (~ $MAP '') + MAP=world + +MAP=$MAP MAPDIR=$MAPDIR $MAPPROG $A $M diff --git a/src/cmd/map/mapdemo.rc b/src/cmd/map/mapdemo.rc new file mode 100755 index 00000000..033969ab --- /dev/null +++ b/src/cmd/map/mapdemo.rc @@ -0,0 +1,83 @@ +#!/bin/rc + +fn demo {proj=$1; shift; + label=$1; shift; + { echo 'o' + echo 'ra -8192 -8492 8192 8492' + echo 'e' + echo 'm -8192 8192' + echo t $type + echo 'm -8192 -8192' + echo t $proj - $label + MAP=world MAPDIR=/lib/map map $proj $* -s -d 5 + } + sleep 5 +} + +rfork en +{ +type='Equatorial projections centered on long. 0. Parallels are straight lines.' + +demo mercator 'equally spaced straight meridians, conformal, straight compass courses' +demo sinusoidal 'equally spaced parallels, equal-area, same as bonne(0)' +demo cylequalarea 'equally spaced straight meridians, equal-area, true scale on Eq' 0 +demo cylindrical 'central projection on tangent cylinder' +demo rectangular 'equally spaced parallels, equally spaced straight meridians, true scale on Eq' 0 +demo gall 'parallels spaced stereographically on prime meridian, equally spaced straight meridians, true scale on Eq' 0 +demo mollweide '(homalographic) equal-area, hemisphere is a circle' +demo gilbert 'globe mapped conformally on hemisphere, viewed orthographically' + +type='Azimuthal: centered on the North Pole, Parallels are concentric circles, Meridians are equally spaced radial lines' + +demo azequidistant 'equally spaced parallels, true distances from pole' +demo azequalarea 'equal area' +demo gnomonic 'central projecton on tangent plane, straight great circles' +demo perspective 'viewed along earth''s axis 2 earth radii from center of earth' 2 +demo orthographic 'viewed from infinity' +demo stereographic 'conformal, projected from opposite pole' +demo laue 'radius = tan(2\(mu colatitude ), used in xray crystallography' +demo fisheye 'fisheye view of stereographic map, index of refraction 2' 2 -o 40.75 74 +demo newyorker 'New Yorker map from viewing pedestal of radius .5' .5 -o 40.75 74 + +type='Polar conic projections symmetric about the Prime Meridian. Parallels are segments of concentric circles.' + +demo conic 'central projection on cone tangent at 40' 40 +demo simpleconic 'equally spaced parallels, true scale on 20 and 50' 20 50 +demo lambert 'conformal, true scale on 20 and 50' 20 50 +demo albers 'equal-area, true scale on 20 and 50' 20 50 +demo bonne 'equally spaced parallels, equal-area, parallel 40 developed from tangent cone' 40 + +type='Projections with bilateral symmetry about the Prime Meridian and the equator.' + +demo polyconic 'parallels developed from tangent cones, equally spaced along Prime Meridian' +demo aitoff 'equal-area projection of globe onto 2-to-1 ellipse, based on azequalarea' +demo lagrange 'conformal, maps whole sphere into a circle' +demo bicentric 'points plotted at true azimuth from two centers on the equator at longitudes +-40, great circles are straight lines' 40 +demo elliptic 'points are plotted at true distance from two centers on the equator at longitudes +-40' 40 +demo globular 'hemisphere is circle, circular meridians and parallels' +demo vandergrinten 'sphere is circle, meridians as in globular, circular arc parallels resemble mercator' + +type='Doubly periodic conformal projections.' + +demo guyou 'W and E hemispheres are square' +demo square 'World is square with Poles at diagonally opposite corners' +demo tetra 'map on tetrahedron with edge tangent to Prime Meridian at S Pole, unfolded into equilateral triangle' +demo hex 'world is hexagon centered on N Pole, N and S hemispheres are equilateral +triangles' + +type='Retroazimuthal projections. Directions to center are true.' + +demo mecca 'equally spaced vertical meridians' 21.4 -o 90 -39.8 +demo homing 'distances to Mecca are true' 21.4 -o 90 -39.8 + +type='Miscellaneous projections.' + +demo harrison 'oblique perspective from above the North Pole, 2 earth radii from the earth, looking along the Date Line 40 degrees off vertical' 2 40 +demo trapezoidal 'equally spaced parallels, straight meridians equally spaced along parallels, true scale at 20 and 50 on Prime Meridian' 20 50 +demo lune 'conformal, polar cap above Eq is 60-degree lune' 0 60 + +type='Maps based on the spheroid' + +demo sp_mercator 'equally spaced straight meridians, conformal' +demo sp_albers 'equal-area, true scale on 20 and 50' 20 50 +} | plot diff --git a/src/cmd/map/mkfile b/src/cmd/map/mkfile new file mode 100644 index 00000000..4d272749 --- /dev/null +++ b/src/cmd/map/mkfile @@ -0,0 +1,32 @@ +<$PLAN9/src/mkhdr + +TARG=mapd +LIB=libmap/libmap.a +OFILES=map.$O\ + symbol.$O\ + index.$O\ + sqrt.$O\ + +HFILES=map.h\ + iplot.h\ + +<$PLAN9/src/mkone + + +$O.out:V: $OFILES $LIB + $LD $LDFLAGS -o $target $prereq + +$LIB:V: + cd libmap + mk install + +installall:V: + for(objtype in $CPUS) + mk install + cp map.rc /rc/bin/map + cp mapdemo.rc /rc/bin/mapdemo + +clean nuke:V: + rm -f *.[$OS] [$OS].out y.tab.? y.debug y.output $TARG + cd libmap; mk clean + diff --git a/src/cmd/map/route.c b/src/cmd/map/route.c new file mode 100644 index 00000000..c4c67134 --- /dev/null +++ b/src/cmd/map/route.c @@ -0,0 +1,131 @@ +#include +#include +#include "map.h" + +/* Given two lat-lon pairs, find an orientation for the + -o option of "map" that will place those two points + on the equator of a standard projection, equally spaced + about the prime meridian. + + -w and -l options are suggested also. + + Option -t prints out a series of + coordinates that follows the (great circle) track + in the original coordinate system, + followed by ". + This data is just right for map -t. + + Option -i inverts the map top-to-bottom. +*/ +struct place pole; +struct coord twist; +int track; +int inv = -1; + +extern void doroute(double, double, double, double, double); + +void +dorot(double a, double b, double *x, double *y, void (*f)(struct place *)) +{ + struct place g; + deg2rad(a,&g.nlat); + deg2rad(b,&g.wlon); + (*f)(&g); + *x = g.nlat.l/RAD; + *y = g.wlon.l/RAD; +} + +void +rotate(double a, double b, double *x, double *y) +{ + dorot(a,b,x,y,normalize); +} + +void +rinvert(double a, double b, double *x, double *y) +{ + dorot(a,b,x,y,invert); +} + +main(int argc, char **argv) +{ +#pragma ref argv + double an,aw,bn,bw; + ARGBEGIN { + case 't': + track = 1; + break; + + case 'i': + inv = 1; + break; + + default: + exits("route: bad option"); + } ARGEND; + if (argc<4) { + print("use route [-t] [-i] lat lon lat lon\n"); + exits("arg count"); + } + an = atof(argv[0]); + aw = atof(argv[1]); + bn = atof(argv[2]); + bw = atof(argv[3]); + doroute(inv*90.,an,aw,bn,bw); + return 0; +} + +void +doroute(double dir, double an, double aw, double bn, double bw) +{ + double an1,aw1,bn1,bw1,pn,pw; + double theta; + double cn,cw,cn1,cw1; + int i,n; + orient(an,aw,0.); + rotate(bn,bw,&bn1,&bw1); +/* printf("b %f %f\n",bn1,bw1);*/ + orient(an,aw,bw1); + rinvert(0.,dir,&pn,&pw); +/* printf("p %f %f\n",pn,pw);*/ + orient(pn,pw,0.); + rotate(an,aw,&an1,&aw1); + rotate(bn,bw,&bn1,&bw1); + theta = (aw1+bw1)/2; +/* printf("a %f %f \n",an1,aw1);*/ + orient(pn,pw,theta); + rotate(an,aw,&an1,&aw1); + rotate(bn,bw,&bn1,&bw1); + if(fabs(aw1-bw1)>180) + if(theta<0.) theta+=180; + else theta -= 180; + orient(pn,pw,theta); + rotate(an,aw,&an1,&aw1); + rotate(bn,bw,&bn1,&bw1); + if(!track) { + double dlat, dlon, t; + /* printf("A %.4f %.4f\n",an1,aw1); */ + /* printf("B %.4f %.4f\n",bn1,bw1); */ + cw1 = fabs(bw1-aw1); /* angular difference for map margins */ + /* while (aw<0.0) + aw += 360.; + while (bw<0.0) + bw += 360.; */ + dlon = fabs(aw-bw); + if (dlon>180) + dlon = 360-dlon; + dlat = fabs(an-bn); + printf("-o %.4f %.4f %.4f -w %.2f %.2f %.2f %.2f \n", + pn,pw,theta, -0.3*cw1, .3*cw1, -.6*cw1, .6*cw1); + + } else { + cn1 = 0; + n = 1 + fabs(bw1-aw1)/.2; + for(i=0;i<=n;i++) { + cw1 = aw1 + i*(bw1-aw1)/n; + rinvert(cn1,cw1,&cn,&cw); + printf("%f %f\n",cn,cw); + } + printf("\"\n"); + } +} diff --git a/src/cmd/map/sqrt.c b/src/cmd/map/sqrt.c new file mode 100644 index 00000000..1f933384 --- /dev/null +++ b/src/cmd/map/sqrt.c @@ -0,0 +1,52 @@ +/* + sqrt returns the square root of its floating + point argument. Newton's method. + + calls frexp +*/ + +#include +#include + +double +sqrt(double arg) +{ + double x, temp; + int exp, i; + + if(arg <= 0) { + if(arg < 0) + return 0.; + return 0; + } + x = frexp(arg, &exp); + while(x < 0.5) { + x *= 2; + exp--; + } + /* + * NOTE + * this wont work on 1's comp + */ + if(exp & 1) { + x *= 2; + exp--; + } + temp = 0.5 * (1.0+x); + + while(exp > 60) { + temp *= (1L<<30); + exp -= 60; + } + while(exp < -60) { + temp /= (1L<<30); + exp += 60; + } + if(exp >= 0) + temp *= 1L << (exp/2); + else + temp /= 1L << (-exp/2); + for(i=0; i<=4; i++) + temp = 0.5*(temp + arg/temp); + return temp; +} diff --git a/src/cmd/map/symbol.c b/src/cmd/map/symbol.c new file mode 100644 index 00000000..e8f3c680 --- /dev/null +++ b/src/cmd/map/symbol.c @@ -0,0 +1,192 @@ +#include +#include +#include +#include "map.h" +#include "iplot.h" + +#define NSYMBOL 20 + +enum flag { POINT,ENDSEG,ENDSYM }; +struct symb { + double x, y; + char name[10+1]; + enum flag flag; +} *symbol[NSYMBOL]; + +static int nsymbol; +static double halfrange = 1; +extern int halfwidth; +extern int vflag; + +static int getrange(FILE *); +static int getsymbol(FILE *, int); +static void setrot(struct place *, double, int); +static void dorot(struct symb *, double *, double *); + + +void +getsyms(char *file) +{ + FILE *sf = fopen(file,"r"); + if(sf==0) + filerror("cannot open", file); + while(nsymbolname)!=1) + return 0; + i = 0; + while(fscanf(sf,"%1s",s)!=EOF) { + switch(s[0]) { + case 'r': + if(!getrange(sf)) + break; + continue; + case 'm': + if(i>0) + symbol[n][i-1].flag = ENDSEG; + continue; + case ':': + ungetc(s[0],sf); + break; + default: + ungetc(s[0],sf); + case 'v': + if(fscanf(sf,"%lf %lf",&x,&y)!=2) + break; + sp[i].x = x*halfwidth/halfrange; + sp[i].y = y*halfwidth/halfrange; + sp[i].flag = POINT; + i++; + sp = symbol[n] = (struct symb*)realloc(symbol[n], + (i+1)*sizeof(struct symb)); + continue; + } + break; + } + if(i>0) + symbol[n][i-1].flag = ENDSYM; + else + symbol[n] = 0; + return 1; +} + +static int +getrange(FILE *sf) +{ + double x,y,xmin,ymin; + if(fscanf(sf,"%*s %lf %lf %lf %lf", + &xmin,&ymin,&x,&y)!=4) + return 0; + x -= xmin; + y -= ymin; + halfrange = (x>y? x: y)/2; + if(halfrange<=0) + error("bad ra command in -y file"); + return 1; +} + +/* r=0 upright;=1 normal;=-1 reverse*/ +int +putsym(struct place *p, char *name, double s, int r) +{ + int x,y,n; + struct symb *sp; + double dx,dy; + int conn = 0; + for(n=0; symbol[n]; n++) + if(strcmp(name,symbol[n]->name)==0) + break; + sp = symbol[n]; + if(sp==0) + return 0; + if(doproj(p,&x,&y)*vflag <= 0) + return 1; + setrot(p,s,r); + for(;;) { + dorot(sp,&dx,&dy); + conn = cpoint(x+(int)dx,y+(int)dy,conn); + switch(sp->flag) { + case ENDSEG: + conn = 0; + case POINT: + sp++; + continue; + case ENDSYM: + break; + } + break; + } + return 1; +} + +static double rot[2][2]; + +static void +setrot(struct place *p, double s, int r) +{ + double x0,y0,x1,y1; + struct place up; + up = *p; + up.nlat.l += .5*RAD; + sincos(&up.nlat); + if(r&&(*projection)(p,&x0,&y0)) { + if((*projection)(&up,&x1,&y1)<=0) { + up.nlat.l -= RAD; + sincos(&up.nlat); + if((*projection)(&up,&x1,&y1)<=0) + goto unit; + x1 = x0 - x1; + y1 = y0 - y1; + } else { + x1 -= x0; + y1 -= y0; + } + x1 = r*x1; + s /= hypot(x1,y1); + rot[0][0] = y1*s; + rot[0][1] = x1*s; + rot[1][0] = -x1*s; + rot[1][1] = y1*s; + } else { +unit: + rot[0][0] = rot[1][1] = s; + rot[0][1] = rot[1][0] = 0; + } +} + +static void +dorot(struct symb *sp, double *px, double *py) +{ + *px = rot[0][0]*sp->x + rot[0][1]*sp->y; + *py = rot[1][0]*sp->x + rot[1][1]*sp->y; +} -- cgit v1.2.3