0707070035051032341006640000030000040000011211070533573712000001100000001722Makefile CFLAGS = -O
CC=cc
OBJ = map.o index.o symbol.o libmap.a
lplot: maplplot index.o symbol.o libmap.a
$(CC) $(OBJ) -lplot -lm -o map
v10: mapv10 index.o symbol.o libmap.a
$(CC) $(OBJ) -lm -o map
ps: maplplot plotPS.o index.o symbol.o libmap.a
$(CC) $(OBJ) plotPS.o -lm -o map
mapv10: map.c iplot.h
$(CC) -c $(CFLAGS) -DPLOT='"iplot.h"' map.c
maplplot: map.c plot.h
$(CC) -c $(CFLAGS) -DPLOT='"plot.h"' map.c
route: route.o libmap.a
$(CC) $(CFLAGS) route.o libmap.a -lm -o route
libmap.a:
cd libmap; make CC=$(CC) CFLAGS="-I.. $(CFLAGS)"
install: map
strip map
test -d /usr/lib/map || mkdir /usr/lib/map
cp map.sh /usr/bin/map
cp map mapdata/* /usr/lib/map
clean:
rm -f map route *.o libmap.a new.results
cd libmap; make clean
quicktest: v10
MAPPROG=./map MAPDIR=./mapdata map.sh mercator -l 0 10 0 10 -g -b >new.results
:
: If any "diff" output follows, it should show only roundoff differences.
diff test.results new.results || true
rm -f map map.o
0707070035051032331006640000030000040000010633140533562750600000700000007126README MATERIALS
Manuals are included in two forms.
formatted ASCII file:
map.man
input for troff -man:
map.5, map.7, proj.3, proj.5, route.1
The source directory has two subdirectories:
libmap source for all the projection subroutines
mapdata World Data Bank I, etc.
The source is written in ANSI C. Check CFLAGS and CC in the
Makefile and make them reflect your system's conventions.
QUICK TEST
For a quick test, to see whether you can compile and compute, try this
test, which makes "map" and checks its output for a bit of Africa:
make quicktest
PLOTTING FILTERS
Map produces output for a plotting filter (not included).
Unfortunately there is no Unix standard for plotting.
Here are ways to compile for various filters.
make For System V and SunOS filters plot(1)
or tplot(1).
make v10 For v10 research system plot(1), not
compatible with System V.
make ps PostScript. Maps are drawn in a 6.5-inch
square centered 1 inch from the top of an
8.5x11 page. To change, edit the
PostScript or plotPS.c.
As map uses only simple plotting features, it is usually
easy to interface to other plotting packages. See
iplot.h (used with "make v10") for details.
When you have a plotting filter, you can test map in the current
directory as follows. For sample arguments, see EXAMPLES on the
map(7) man page.
MAPDIR=./mapdata MAPPROG=./map map.sh arguments | filter
INSTALLATION
For real installation, examine the recipe for
make install
It puts the map shell script in /usr/bin and everything else
in /usr/lib/map. If you want to put things elsewhere, adjust
variables MAPDIR and MAPPROG in map.sh. map.sh will become
the command "map". After installation it will be run thus:
map arguments | filter
POSSIBLE PROBLEMS
Options
Option -f oes not work because World Data Bank II is
not in this distribution on account of its size (13Mb).
-y files are like v10, not Sys V, plot files.
The plot(5) man page from v10 is included for reference,
although most of it is irrelevant for this application.
Library
At least one version of tplot(1) has been seen to garble
the output. In this case the trouble went away by
compiling map with -l4014 instead of -lplot; see plot(3).
Collisions with a nonstandard library function sincos()
have been observed. sincos will have to be renamed
in map.h and all .c files.
Man pages
The man pages were produced with a -man macros slightly
different from Sys V. If you change font L to B throughout
you will probably get a passable result.
REFERENCES
Most standard texts on cartography discuss the major projections
and their uses. Some thorough works:
J. P. Snyder, An Album of Map Projections, US Geological Survey
Professional Paper 1453, USGPO Washington (1989)
J. P. Snyder, Map Projections - A Working Manual, US Geological
Survey Professional Paper 1395, USGPO Washington (1987)
D. H. Maling, Coordinate Systems and Map Projections, George
Philip and Son, London (1973)
J. A. Steers, An Introduction to the Study of Map Projections,
Univ. London Press (1970)
A short introduction, with 24 sample maps drawn by the "map"
program itself and the commands to draw them:
M. D. McIlroy, Projections: Mapmakers' Answers to the Riddle of
Presenting a Round Earth on Flat Paper, AT&T Bell Labs Computing
Science Tech. Report 140 (1987), order from Computing Science
Reports, Room 2C579, AT&T Bell Labs Murray Hill, NJ 07974
[In this report, option -s is used in an obsolete way. It
should always be replaced by -s2, with -s1 added to the
preceding command.]
Doug McIlroy
201-582-6050
research!doug
[email protected]
0707070035051032721006640000030000040000020633020533472714100001000000010405index.c #include "map.h"
/* p0;p1 evaluated to shut off "not used" diags */
static proj Yaitoff(double p0,double p1){p0;p1;return aitoff();}
static proj Yalbers(double p0,double p1){return albers(p0,p1);}
static proj Yazequalarea(double p0,double p1){p0;p1;return azequalarea();}
static proj Yazequidistant(double p0,double p1){p0;p1;return azequidistant();}
static proj Ybicentric(double p0,double p1){p1;return bicentric(p0);}
static proj Ybonne(double p0,double p1){p1;return bonne(p0);}
static proj Yconic(double p0,double p1){p1;return conic(p0);}
static proj Ycylequalarea(double p0,double p1){p1;return cylequalarea(p0);}
static proj Ycylindrical(double p0,double p1){p0;p1;return cylindrical();}
static proj Yelliptic(double p0,double p1){p1;return elliptic(p0);}
static proj Yfisheye(double p0,double p1){p1;return fisheye(p0);}
static proj Ygall(double p0,double p1){p1;return gall(p0);}
static proj Yglobular(double p0,double p1){p0;p1;return globular();}
static proj Ygnomonic(double p0,double p1){p0;p1;return gnomonic();}
static proj Yguyou(double p0,double p1){p0;p1;return guyou();}
static proj Yharrison(double p0,double p1){return harrison(p0,p1);}
static proj Yhex(double p0,double p1){p0;p1;return hex();}
static proj Yhoming(double p0,double p1){p1;return homing(p0);}
static proj Ylagrange(double p0,double p1){p0;p1;return lagrange();}
static proj Ylambert(double p0,double p1){return lambert(p0,p1);}
static proj Ylaue(double p0,double p1){p0;p1;return laue();}
static proj Ymecca(double p0,double p1){p1;return mecca(p0);}
static proj Ymercator(double p0,double p1){p0;p1;return mercator();}
static proj Ymollweide(double p0,double p1){p0;p1;return mollweide();}
static proj Ynewyorker(double p0,double p1){p1;return newyorker(p0);}
static proj Yorthographic(double p0,double p1){p0;p1;return orthographic();}
static proj Yperspective(double p0,double p1){p1;return perspective(p0);}
static proj Ypolyconic(double p0,double p1){p0;p1;return polyconic();}
static proj Yrectangular(double p0,double p1){p1;return rectangular(p0);}
static proj Ysimpleconic(double p0,double p1){return simpleconic(p0,p1);}
static proj Ysinusoidal(double p0,double p1){p0;p1;return sinusoidal();}
static proj Ysp_albers(double p0,double p1){return sp_albers(p0,p1);}
static proj Ysp_mercator(double p0,double p1){p0;p1;return sp_mercator();}
static proj Ysquare(double p0,double p1){p0;p1;return square();}
static proj Ystereographic(double p0,double p1){p0;p1;return stereographic();}
static proj Ytetra(double p0,double p1){p0;p1;return tetra();}
static proj Ytrapezoidal(double p0,double p1){return trapezoidal(p0,p1);}
static proj Yvandergrinten(double p0,double p1){p0;p1;return vandergrinten();}
struct index index[] = {
{"aitoff", Yaitoff, 0, picut, 0, 0},
{"albers", Yalbers, 2, picut, 3, 0},
{"azequalarea", Yazequalarea, 0, nocut, 1, 0},
{"azequidistant", Yazequidistant, 0, nocut, 1, 0},
{"bicentric", Ybicentric, 1, nocut, 0, 0},
{"bonne", Ybonne, 1, picut, 0, 0},
{"conic", Yconic, 1, picut, 0, 0},
{"cylequalarea", Ycylequalarea, 1, picut, 3, 0},
{"cylindrical", Ycylindrical, 0, picut, 0, 0},
{"elliptic", Yelliptic, 1, picut, 0, 0},
{"fisheye", Yfisheye, 1, nocut, 0, 0},
{"gall", Ygall, 1, picut, 3, 0},
{"globular", Yglobular, 0, picut, 0, 0},
{"gnomonic", Ygnomonic, 0, nocut, 0, 0},
{"guyou", Yguyou, 0, guycut, 0, 0},
{"harrison", Yharrison, 2, nocut, 0, 0},
{"hex", Yhex, 0, hexcut, 0, 0},
{"homing", Yhoming, 1, picut, 0, 0},
{"lagrange", Ylagrange,0,picut,0, 0},
{"lambert", Ylambert, 2, picut, 0, 0},
{"laue", Ylaue, 0, nocut, 0, 0},
{"mecca", Ymecca, 1, picut, 0, 0},
{"mercator", Ymercator, 0, picut, 0, 0},
{"mollweide", Ymollweide, 0, picut, 0, 0},
{"newyorker", Ynewyorker, 1, nocut, 0, 0},
{"orthographic", Yorthographic, 0, nocut, 0, 0},
{"perspective", Yperspective, 1, nocut, 0, 0},
{"polyconic", Ypolyconic, 0, picut, 0, 0},
{"rectangular", Yrectangular, 1, picut, 3, 0},
{"simpleconic", Ysimpleconic, 2, picut, 3, 0},
{"sinusoidal", Ysinusoidal, 0, picut, 0, 0},
{"sp_albers", Ysp_albers, 2, picut, 3, 1},
{"sp_mercator", Ysp_mercator, 0, picut, 0, 1},
{"square", Ysquare, 0, picut, 0, 0},
{"stereographic", Ystereographic, 0, nocut, 0, 0},
{"tetra", Ytetra, 0, tetracut, 0, 0},
{"trapezoidal", Ytrapezoidal, 2, picut, 3, 0},
{"vandergrinten", Yvandergrinten, 0, picut, 0, 0},
0
};
0707070035051032701006640000030000040000020035130417420602000001000000002016iplot.h /* Plotting functions for v8 and v9 systems */
/* This file is an alternative to plot.h */
/* open the plotting output */
#define openpl() printf("o\n")
/* close the plotting output */
#define closepl() printf("cl\n")
/* make sure the page or screen is clear */
#define erase() printf("e\n")
/* plot a point at _x,_y, which becomes current */
#define point(_x,_y) printf("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) printf("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) == ' ')printf("t \"%s\"\n",_s); else printf("t %s\n", _s); }
/* draw line from current point to _x,_y, which becomes current */
#define vec(_x,_y) printf("v %d %d\n", _x,_y)
/* _x,_y becomes current point */
#define move(_x, _y) printf("m %d %d\n", _x, _y)
/* specify style for drawing lines: "dotted", "solid", "dotdash" */
#define pen(_s) printf("pe %s\n", _s)
0707070035051031101006640000030000040000011401020525002214000002000000001206libmap/Makefile
# .o's ordered for a nonrandom library
OBJ= aitoff.o albers.o azequalarea.o elliptic.o azequidist.o \
bicentric.o bonne.o conic.o cylequalarea.o cylindrical.o fisheye.o \
gall.o guyou.o harrison.o hex.o homing.o lagrange.o lambert.o \
laue.o mercator.o mollweide.o newyorker.o polyconic.o simpleconic.o \
sinusoidal.o tetra.o perspective.o orthographic.o trapezoidal.o \
rectangular.o twocirc.o cuts.o ccubrt.o cubrt.o elco2.o complex.o zcoord.o
# ignore error on systems without ranlib
../libmap.a: ../map.h $(OBJ)
ar cr ../libmap.a $(OBJ)
ranlib ../libmap.a 2>/dev/null || true
clean:
rm -f *.o
%.o: %.c
$(CC) $(CFLAGS) -I.. -c $*.c
0707070035051031071006640000030000040000011231100533472766300002000000000554libmap/aitoff.c #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);
}
0707070035051031061006640000030000040000011400540533472766300002000000004555libmap/albers.c #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 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,r2;
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));
r2 = plat2.c/(n*sqrt(2-e2*plat2.s*plat2.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;
}
0707070035051031051006640000030000040000011400560533472766300002500000000361libmap/azequalarea.c #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);
}
0707070035051031041006640000030000040000011400600533472766400002400000000401libmap/azequidist.c #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);
}
0707070035051031031006640000030000040000011400620533472766400002300000000624libmap/bicentric.c #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);
}
0707070035051031021006640000030000040000011400640533472766400001700000001164libmap/bonne.c #include "map.h"
static struct coord stdpar;
static double r0;
static
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);
}
0707070035051031011006640000030000040000011400660470474632500002000000000301libmap/ccubrt.c #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);
}
0707070035051031001006640000030000040000011400700470474632500002100000002535libmap/complex.c #include "map.h"
/*complex divide, defensive against overflow from
* * and /, but not from + and -
* assumes underflow yields 0.0
* uses identities:
* (a + bi)/(c + di) = ((a + bd/c) + (b - ad/c)i)/(c + dd/c)
* (a + bi)/(c + di) = (b - ai)/(d - ci)
*/
void
cdiv(double a, double b, double c, double d, double *u, double *v)
{
double r,t;
if(fabs(c)<fabs(d)) {
t = -c; c = d; d = t;
t = -a; a = b; b = t;
}
r = d/c;
t = c + r*d;
*u = (a + r*b)/t;
*v = (b - r*a)/t;
}
void
cmul(double c1, double c2, double d1, double d2, double *e1, double *e2)
{
*e1 = c1*d1 - c2*d2;
*e2 = c1*d2 + c2*d1;
}
void
csq(double c1, double c2, double *e1, double *e2)
{
*e1 = c1*c1 - c2*c2;
*e2 = c1*c2*2;
}
/* complex square root
* assumes underflow yields 0.0
* uses these identities:
* sqrt(x+_iy) = sqrt(r(cos(t)+_isin(t))
* = sqrt(r)(cos(t/2)+_isin(t/2))
* cos(t/2) = sin(t)/2sin(t/2) = sqrt((1+cos(t)/2)
* sin(t/2) = sin(t)/2cos(t/2) = sqrt((1-cos(t)/2)
*/
void
csqrt(double c1, double c2, double *e1, double *e2)
{
double r,s;
double x,y;
x = fabs(c1);
y = fabs(c2);
if(x>=y) {
if(x==0) {
*e1 = *e2 = 0;
return;
}
r = x;
s = y/x;
} else {
r = y;
s = x/y;
}
r *= sqrt(1+ s*s);
if(c1>0) {
*e1 = sqrt((r+c1)/2);
*e2 = c2/(2* *e1);
} else {
*e2 = sqrt((r-c1)/2);
if(c2<0)
*e2 = -*e2;
*e1 = c2/(2* *e2);
}
}
0707070035051030771006640000030000040000011400720533472766500001700000000717libmap/conic.c #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);
}
0707070035051030761006640000030000040000011400740470474632600001700000000450libmap/cubrt.c #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);
}
0707070035051030451006640000030000040000011400760533472766500001600000001453libmap/cuts.c #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;
}
0707070035051030751006640000030000040000011401000533472766500002600000000502libmap/cylequalarea.c #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);
}
0707070035051030741006640000030000040000011401100533472766600002500000000376libmap/cylindrical.c #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);
}
0707070035051030731006640000030000040000011401120533472766600001700000004704libmap/elco2.c #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);
}
0707070035051030721006640000030000040000010612660533472766600002200000001122libmap/elliptic.c #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);
}
0707070035050603501006640000100000040000011210700533472766600002100000000566libmap/fisheye.c #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;
}
0707070035051030401006640000030000040000011401140533472766700001600000000737libmap/gall.c #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;
}
0707070035051030701006640000030000040000011402000533472766700001700000003271libmap/guyou.c #include "map.h"
static struct place gywhem, gyehem;
static struct coord gytwist;
static double gyconst, gykc, gyside;
static void
dosquare(double z1, double z2, double *x, double *y)
{
double w1,w2;
w1 = z1 -1;
if(fabs(w1*w1+z2*z2)>.000001) {
cdiv(z1+1,z2,w1,z2,&w1,&w2);
w1 *= gyconst;
w2 *= gyconst;
if(w1<0)
w1 = 0;
elco2(w1,w2,gykc,1.,1.,x,y);
} else {
*x = gyside;
*y = 0;
}
}
int
Xguyou(struct place *place, double *x, double *y)
{
int ew; /*which hemisphere*/
double z1,z2;
struct place pl;
ew = place->wlon.l<0;
copyplace(place,&pl);
norm(&pl,ew?&gyehem:&gywhem,&gytwist);
Xstereographic(&pl,&z1,&z2);
dosquare(z1/2,z2/2,x,y);
if(!ew)
*x -= gyside;
return(1);
}
proj
guyou(void)
{
double junk;
gykc = 1/(3+2*sqrt(2.));
gyconst = -(1+sqrt(2.));
elco2(-gyconst,0.,gykc,1.,1.,&gyside,&junk);
gyside *= 2;
latlon(0.,90.,&gywhem);
latlon(0.,-90.,&gyehem);
deg2rad(0.,&gytwist);
return(Xguyou);
}
int
guycut(struct place *g, struct place *og, double *cutlon)
{
int c;
c = picut(g,og,cutlon);
if(c!=1)
return(c);
*cutlon = 0.;
if(g->nlat.c<.7071||og->nlat.c<.7071)
return(ckcut(g,og,0.));
return(1);
}
static int
Xsquare(struct place *place, double *x, double *y)
{
double z1,z2;
double r, theta;
struct place p;
copyplace(place,&p);
if(place->nlat.l<0) {
p.nlat.l = -p.nlat.l;
p.nlat.s = -p.nlat.s;
}
if(p.nlat.l<FUZZ && fabs(p.wlon.l)>PI-FUZZ){
*y = -gyside/2;
*x = p.wlon.l>0?0:gyside;
return(1);
}
Xstereographic(&p,&z1,&z2);
r = sqrt(sqrt(hypot(z1,z2)/2));
theta = atan2(z1,-z2)/4;
dosquare(r*sin(theta),-r*cos(theta),x,y);
if(place->nlat.l<0)
*y = -gyside - *y;
return(1);
}
proj
square(void)
{
guyou();
return(Xsquare);
}
0707070035051030501006640000030000040000010631010533472766700002200000001241libmap/harrison.c #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 0;
t = a/d;
if(t < 0)
return 0;
if(v3*place->nlat.s < 1.)
return -1;
*y = t*p2*u2 + (v3-t*(v3-p3))*u3;
*x = t*p1;
if(*x * *x + *y * *y > 16)
return 0;
return 1;
}
proj
harrison(double r, double alpha)
{
u2 = cos(alpha*RAD);
u3 = sin(alpha*RAD);
v3 = r;
b = r*u2;
a = 1 + b;
if(r<1.001 || a<sqrt(r*r-1))
return 0;
return Xharrison;
}
0707070035051030671006640000030000040000011402020533472766700001500000004312libmap/hex.c #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;
register i;
double zr,zi;
double sr,si,tr,ti,ur,ui,vr,vi,yr,yi;
struct place p;
copyplace(place,&p);
ns = place->nlat.l >= 0;
if(!ns) {
p.nlat.l = -p.nlat.l;
p.nlat.s = -p.nlat.s;
}
if(p.nlat.l<HFUZZ) {
for(i=0;i<3;i++)
if(fabs(reduce(p.wlon.l-hcut[i]))<HFUZZ) {
if(i==2) {
*x = 2*cr[0] - cr[1];
*y = 0;
} else {
*x = cr[1];
*y = 2*ci[2*i];
}
return(1);
}
p.nlat.l = HFUZZ;
sincos(&p.nlat);
}
norm(&p,&hem,&twist);
Xstereographic(&p,&zr,&zi);
zr /= 2;
zi /= 2;
cdiv(1-zr,-zi,1+zr,zi,&sr,&si);
csq(sr,si,&tr,&ti);
ccubrt(1+3*tr,3*ti,&ur,&ui);
csqrt(ur-1,ui,&vr,&vi);
cdiv(rootroot3+vr,vi,rootroot3-vr,-vi,&yr,&yi);
yr /= rootk;
yi /= rootk;
elco2(fabs(yr),yi,hkc,1.,1.,x,y);
if(yr < 0)
*x = w2 - *x;
if(!ns) reflect(hcut[0]>place->wlon.l?0:
hcut[1]>=place->wlon.l?1:
2,*x,*y,x,y);
return(1);
}
proj
hex(void)
{
register 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);
}
0707070035051030511006640000030000040000011402040533472767300002000000002430libmap/homing.c #include "map.h"
static struct coord p0; /* standard parallel */
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 < .01)
return 0;
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? 0:
rad.c<0? -1:
1;
}
proj
mecca(double par)
{
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? -1: 1;
}
proj
homing(double par)
{
if(fabs(par)>80.)
return(0);
deg2rad(par,&p0);
return(Xhoming);
}
0707070035051030351006640000030000040000011402060533472767300002200000000663libmap/lagrange.c #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);
}
0707070035051030661006640000030000040000011402100533472767300002100000001576libmap/lambert.c #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);
}
0707070035051030651006640000030000040000011402120533472767500001600000000445libmap/laue.c #include "map.h"
static int
Xlaue(struct place *place, double *x, double *y)
{
double r;
if(place->nlat.l<PI/4+FUZZ)
return(-1);
r = tan(PI-2*place->nlat.l);
if(r>3)
return(-1);
*x = - r * place->wlon.s;
*y = - r * place->wlon.c;
return(1);
}
proj
laue(void)
{
return(Xlaue);
}
0707070035051030631006640000030000040000011402140533472767500002200000001026libmap/mercator.c #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);
}
0707070035051030621006640000030000040000011402160533472767500002300000000616libmap/mollweide.c #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);
}
0707070035051030521006640000030000040000010633000533472767600002300000000570libmap/newyorker.c #include "map.h"
static double a;
static int
Xnewyorker(struct place *place, double *x, double *y)
{
double r = PI/2 - place->nlat.l;
double s;
if(r<.001) /* cheat to plot center */
s = 0;
else if(r<a)
return -1;
else
s = log(r/a);
*x = -s * place->wlon.s;
*y = -s * place->wlon.c;
return(1);
}
proj
newyorker(double a0)
{
a = a0*RAD;
return(Xnewyorker);
}
0707070035051030611006640000030000040000011402220533472767600002600000000370libmap/orthographic.c #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);
}
0707070035051030601006640000030000040000011402240533472767600002500000001435libmap/perspective.c #include "map.h"
static double viewpt;
static int
Xperspective(struct place *place, double *x, double *y)
{
double r;
if(viewpt<=1. && 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(0);
if(fabs(viewpt)>1. && place->nlat.s<=1./viewpt)
return(0);
return(1);
}
int
Xstereographic(struct place *place, double *x, double *y)
{
viewpt = -1;
return(Xperspective(place, x, y));
}
proj
perspective(double radius)
{
viewpt = radius;
if(radius>=1000.)
return(Xorthographic);
if(fabs(radius-1.)<.0001)
return(0);
return(Xperspective);
}
proj
stereographic(void)
{
viewpt = -1.;
return(Xperspective);
}
proj
gnomonic(void)
{
viewpt = 0.;
return(Xperspective);
}
0707070035051030571006640000030000040000011402260533472767700002300000001047libmap/polyconic.c #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);
}
0707070035051030561006640000030000040000011402300533472770000002500000000426libmap/rectangular.c #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);
}
0707070035051030361006640000030000040000011402320533472770000002500000001143libmap/simpleconic.c #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;
}
0707070035051030551006640000030000040000011402360533472770000002400000000312libmap/sinusoidal.c #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);
}
0707070035051030541006640000030000040000011402400533472770100001700000011457libmap/tetra.c #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)
{
register 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);
}
0707070035051030711006640000030000040000011402440533472770100002500000001041libmap/trapezoidal.c #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;
}
0707070035051030431006640000030000040000011400500533472770100002100000003153libmap/twocirc.c #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;
}
0707070035051030531006640000030000040000011402460533501406400002000000004626libmap/zcoord.c #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;
}
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;
}
0707070035050172611006640000030000040000010336560447134002200000600000003321map.5 .TH MAP 5
.CT 1 inst_info graphics
.SH NAME
map \- digitized map formats
.SH DESCRIPTION
Files used by
.IR map (7)
are a sequence of structures of the form:
.PP
.EX
struct {
signed char patchlatitude;
signed char patchlongitude;
short n;
union {
struct {
short latitude;
short longitude;
} point[n];
struct {
short latitude;
short longitude;
struct {
signed char latdiff;
signed char londiff;
} point[\-n];
} highres;
} segment;
};
.EE
.PP
Fields
.L patchlatitude
and
.L patchlongitude
tell to what
10-degree by 10-degree
patch of the earth's surface a segment belongs.
Their values range from \-9 to 8 and from \-18 to 17,
respectively, and indicate the coordinates of the
southeast corner of the patch in units of 10 degrees.
.PP
Each segment of
.RB | n |
points is connected; consecutive segments
are not necessarily related.
Latitude and longitude
are measured in units of 0.0001 radian.
If
.B n
is negative, then
differences to the first and succeeding points
are measured in units of 0.00001 radian.
Latitude is counted positive to the north and
longitude positive to the west.
.PP
The patches are ordered lexicographically by
.L patchlatitude
then
.LR patchlongitude .
A printable
index to the first segment of each patch
in a file named
.I data
is kept in an associated file named
.IB data .x .
Each line of an index file contains
.L patchlatitude,
.L patchlongitude
and the byte position
of the patch
in the map file.
Both the map file and the index file are ordered by
patch latitude and longitude.
.PP
Shorts are stored in little-endian order, low byte first,
regardless of computer architecture.
To assure portability,
.I map
accesses them bytewise.
.SH "SEE ALSO"
.IR map (7),
.IR proj (3)
0707070035050172351006640000030000040000010270460533572105600000600000021604map.7 .TH MAP 7
.CT 1 inst_info
.SH NAME
map \- draw maps on various projections
.SH SYNOPSIS
.B map
.I projection
[
.I param ...
]
[
.I option ...
]
.PP
.SH DESCRIPTION
.I Map
prepares on the standard output a
map suitable for display by any
plotting filter described in
.IR plot (1).
A menu of projections is produced in response to an unknown
.IR projection .
For the meanings of
.I params
pertinent to particular projections
see
.IR proj (3).
.PP
The default data for
.I map
are world shorelines.
Option
.B -f
accesses the higher-resolution World Data Bank II.
.TP
.BR -f " [ \fIfeature\fR ... ]"
Features are ranked 1 (default) to 4 from major to minor.
Higher-numbered ranks include all lower-numbered ones.
Features are
.RS
.TF country[1-3]
.TP
.BR shore [ 1 - 4 ]
seacoasts, lakes, and islands; in the absence of
.BR -m ,
option
.B -f
automatically includes
.B shore1
.TP
.BR ilake [ 1 - 2 ]
intermittent lakes
.TP
.BR river [ 1 - 4 ]
rivers
.TP
.BR iriver [ 1 - 3 ]
intermittent rivers
.TP
.BR canal [ 1 - 3 ]
.BR 3 =irrigation
canals
.TP
.BR glacier
.TP
.BR iceshelf [ 12 ]
.TP
.BR reef
.TP
.BR saltpan [ 12 ]
.TP
.BR country [ 1 - 3 ]
.BR 2 =disputed
boundaries,
.BR 3 =indefinite
boundaries
.TP
.BR state
states and provinces (US and Canada only)
.PD
.RE
.PP
In other options
coordinates are in degrees, with north latitude
and west longitude counted as positive.
.TP 0
.BI -l " S N E W"
Set the southern and northern latitude
and the eastern and western longitude limits.
Missing arguments are filled out from the list
\-90, 90, \-180, 180.
.TP
.BI -k " S N E W
Set the scale as if for a map with limits
.B -l
.I "S N E W"
and no
.B -w
option.
.TP
.BI -o " lat lon rot"
Orient the map in a nonstandard position.
Imagine a transparent gridded sphere around the globe.
Turn the overlay about the North Pole
so that the Prime Meridian (longitude 0)
of the overlay coincides with meridian
.I lon
on the globe.
Then tilt the North Pole of the
overlay along its Prime Meridian to latitude
.I lat
on the globe.
Finally again turn the
overlay about its `North Pole' so
that its Prime Meridian coincides with the previous position
of meridian
.IR rot .
Project the map in
the standard form appropriate to the overlay, but presenting
information from the underlying globe.
Missing arguments are filled out from the list
90, 0, 0.
In the absence of
.BR \-o ,
the orientation is 90, 0,
.I m,
where
.I m
is the middle of the longitude range.
.TP
.BI -w " S N E W"
Window the map by the specified latitudes
and longitudes in the tilted, rotated coordinate system.
Missing arguments are filled out from the list \-90, 90, \-180, 180.
(It is wise to give an encompassing
.B -l
option with
.BR -w .
Otherwise for small windows computing time
varies inversely with area!)
.TP
.BI -d " n"
For speed, plot only every
.IR n th
point.
.TP
.B -r
Reverse left and right
(good for star charts and inside-out views).
.br
.ns
.TP
.B -s1
.br
.ns
.TP
.B -s2
Superpose. Outputs for a
.B -s1
map (no closing) and a
.B -s2
map (no opening) may be concatenated.
.TP
.BI -g " dlat dlon res"
Grid spacings are
.I dlat,
.I dlon.
Zero spacing means no grid.
Missing
.I dlat
is taken to be zero.
Missing
.I dlon
is taken the same as
.IR dlat .
Grid lines are drawn to a resolution of
.I res
(2\(de or less by default).
In the absence of
.BR \-g ,
grid spacing is 10\(de.
.TP
.BI -p " lat lon extent"
Position the point
.I lat, lon
at the center of a square plotting area.
Scale the map so that a side of the square is
.I extent
times the size of one degree of latitude
at the center.
By default maps are scaled and positioned
to fit within the plotting area.
An
.I extent
overrides option
.BR -k .
.TP
.BI -c " x y rot"
After all other positioning and scaling operations,
rotate the image
.I rot
degrees counterclockwise about the center
and move the center to position
.I x, y,
of the plotting area, whose nominal extent is
.RI \-1 \(<= x \(<= 1,
.RI \-1 \(<= y \(<= 1.
The map is clipped to this area.
Missing arguments are taken to be 0.
.TP
.BR -m " [ \fIfile\fP ... ]"
Use
map data from named files.
If no files are named, omit map data.
Files that cannot be found directly are looked up
a standard directory, which contains, in addition to the
data for
.BR -f ,
.RS
.LP
.TF counties
.TP
.B world
World Data Bank I from CIA (default)
.TP
.B states
US map from Census Bureau
.TP
.B counties
US map from Census Bureau
.PD
.RE
.IP
The environment variables
.B MAP
and
.B MAPDIR
change the default
map and default directory.
.TP
.BI -b " \fR[ \fPlat1 lon1 lat2 lon2 \fR... ]"
Suppress the drawing of the normal boundary
(defined by options
.BR -l
and
.BR -w ).
Coordinates, if present, define the vertices of a
polygon to which the map is clipped.
If only two vertices are given, they are taken to be the
diagonal of a rectangle.
To draw the polygon, give its vertices as a
.B -u
track.
.TP
.BI -t " file ..."
The arguments name ASCII files that
contain lists of points,
given as latitude-longitude pairs in degrees.
If the first file is named
.LR - ,
the standard input is taken instead.
The points of each list are plotted as connected `tracks'.
.IP
Points in a track file may be followed by label strings.
A label breaks the track.
A label may be prefixed by
\f5"\fR,
.LR : ,
or
.L !
and is terminated by a newline.
An unprefixed string or a string prefixed with
.L
"
is displayed at the designated point.
The first word of a
.L :
or
.L !
string names a special symbol (see option
.BR -y ).
An optional numerical second word is a scale factor
for the size of the symbol, 1 by default.
A
.L :
symbol is aligned with its top to the north; a
.L !
symbol is aligned vertically on the page.
.TP
.BI -u " file ..."
Same as
.BR -t ,
except the tracks are
unbroken lines.
.RB ( -t
tracks are dot-dash lines.)
.TP
.BI -y " file
The
.I file
contains
.IR plot (5)-style
data for
.L :
or
.L !
labels in
.B -t
or
.B -u
files.
Each symbol is defined by a comment
.BI : name
then a sequence of
.L m
and
.L v
commands.
Coordinates (0,0) fall on the plotting point.
Default scaling is as if the nominal plotting range were
.LR "ra -1 -1 1 1" ;
.L ra
commands in
.I file
change the scaling.
.SH EXAMPLES
.TP
.L
map perspective 1.025 -o 40.75 74
A view looking down on New York from 100 miles
(0.025 of the 4000-mile earth radius).
The job can be done faster by limiting the map so as not to `plot'
the invisible part of the world:
.LR "map perspective 1.025 -o 40.75 74 -l 20 60 30 100".
A circular border can be forced by adding option
.LR "-w 77.33" .
(Latitude 77.33\(de falls just inside a polar cap of
opening angle arccos(1/1.025) = 12.6804\(de.)
.TP
.L
map mercator -o 49.25 -106 180
A map whose `equator' is a great circle pasing east-west
through New York.
The pole of the map is placed 90\(de away (40.75+49.25=90)
on the
other side of the earth.
A 180\(de twist around the pole of the map arranges that the
Prime Meridian of the map runs from the pole of the
map over the North Pole to New York
instead of down the back side of the earth.
The same effect can be had from
.L
map mercator -o 130.75 74
.TP
.L
map albers 28 45 -l 20 50 60 130 -m states
A customary curved-latitude map of the United States.
.TP
.L
map albers 28 45 -l 20 50 60 130 -y yfile -t tfile
An example of tracks, labels, and symbols.
Arrows at New York and Miami are 8% and 12%
as long as the map is wide.
The contents of
.L yfile
and
.L tfile
are
.nf
.ft L
.ta 3i
ra -50 -50 50 50 25.77 80.20 :arrow 12
:arrow 25.77 80.20 Miami
m -1 0 25.77 80.20
v 0 0 35.00 74.02
v -.6 .3 40.67 74.02 !arrow 8
m -.6 -.3 40.67 74.02 " New York
v 0 0 34.05 118.25 Los Angeles
.ft
.TP
.L
map harrison 2 30 -l -90 90 120 240 -o 90 0 0
A fan view covering 60\(de on either
side of the Date Line, as seen from one earth radius
above the North Pole gazing at the
earth's limb, which is 30\(de off vertical.
Option
.B -o
overrides the default
.BR "-o 90 0 180" ,
which would rotate
the scene to behind the observer.
.SH FILES
All files in directory $MAPDIR
.TF counties
.TP
.F [1-4]??
World Data Bank II for option
.B -f
.TP
.BR world , states , counties
default and other maps for option
.B -m
.TP
.F *.x
map indexes
.TP
.F map
the program proper
.SH "SEE ALSO"
.IR map (5),
.IR proj (3),
.IR plot (1)
.SH DIAGNOSTICS
`Map seems to be empty'\(ema coarse survey found
zero extent within the
.B -l
and
.BR -w
bounds; for maps of limited extent
the grid resolution,
.I res,
or the limits may have to be refined.
.SH BUGS
The syntax of range specifications in
.B -y
files differs from that in options.
.br
Windows (option
.BR -w )
cannot cross the Date Line.
.br
No borders appear along edges arising from
visibility limits.
.br
Segments that cross a border are dropped, not clipped.
.br
Certain very long line segments are dropped on the assumption
that they were intended to go the other way around the world.
.br
Automatic scaling may miss the extreme points of
peculiarly shaped maps; use option
.B -p
to recover.
.br
Although
.I map
draws grid lines dotted and
.B -t
tracks dot-dashed, many plotting filters
cannot cope and make them solid.
0707070035051031151006640000030000040000020633030533574341400000600000055751map.c #include "map.h"
#define NTRACK 10 /* max number of -t and -u files */
#define NFILE 30 /* max number of map files */
#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 = "/usr/dict"; /* default map directory */
static char *file[NFILE+1] = { /* list of map files */
"world", /* default map */
0
};
extern struct index index[];
int halfwidth = HALFWIDTH;
static int (*cut)(struct place *, struct place *, double *);
static int poles;
static double orientation[3] = { 90., 0., 0. }; /* -o option */
static oriented; /* nonzero if -o option occurred */
static 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 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 mflag = 0; /* 1 if option -m occurred */
static int kflag = 0; /* 1 if option -k 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 track { /* options -t and -u */
int tracktyp; /* 't' or 'u' */
char *tracknam; /* name of input file */
} track[NTRACK];
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 track *);
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;
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] = 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<index[i].npar; k++)
fprintf(stderr," p%d", k);
fprintf(stderr,"\n");
}
exit(1);
found:
argv += 2;
argc -= 2;
cut = index[k].cut;
poles = index[k].poles;
for(i=0;i<index[k].npar;i++) {
if(i>=argc||option(argv[i])) {
fprintf(stderr,"%s needs %d params\n",index[k].name,index[k].npar);
exit(1);
}
params[i] = atof(argv[i]);
}
argv += i;
argc -= i;
while(argc>0&&option(argv[0])) {
argc--;
argv++;
switch(argv[-1][1]) {
case 'm':
i = 0;
if(!mflag) while(file[i]!=0)
i++;
for(i=0;i<NFILE&&argc>i&&!option(argv[i]);i++)
file[i] = argv[i];
file[i] = 0;
mflag++;
argc -= i;
argv += i;
break;
case 'b':
bflag = 0;
for(nvert=0;nvert<NVERT&&argc>=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':
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':
case 'u':
for(i=0;ntrack<NTRACK&&argc>i&&!option(argv[i]);i++) {
track[ntrack].tracktyp = argv[-1][1];
track[ntrack++].tracknam = argv[i];
}
argc -= i;
argv +=i;
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&&i<argc&&!option(argv[i]);i++)
orientation[i] = atof(argv[i]);
oriented++;
argv += i;
argc -= i;
break;
case 'l':
for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
limits[i] = atof(argv[i]);
argv += i;
argc -= i;
break;
case 'k':
kflag++;
for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
klimits[i] = atof(argv[i]);
argv += i;
argc -= i;
break;
case 'd':
if(argc>0&&!option(argv[0])) {
delta = atoi(argv[0]);
argv++;
argc--;
}
break;
case 'w':
windowed++;
for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
window[i] = atof(argv[i]);
argv += i;
argc -= i;
break;
case 'c':
for(i=0;i<3&&argc>i&&!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;
}
}
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];lat<klimits[1]+dd-FUZZ;lat+=dd) {
if(lat>klimits[1])
lat = klimits[1];
for(lon=klimits[2];lon<klimits[3]+dd-FUZZ;lon+=dd) {
if((kflag?posproj:normproj)
(lat,lon+(lon<klimits[3]?FUZZ:-FUZZ),
&x,&y)<=0)
continue;
if(x<xmin) xmin = x;
if(x>xmax) xmax = x;
if(y<ymin) ymin = y;
if(y>ymax) ymax = y;
}
}
} else {
for(i=0; i<nvert; i++) {
x = v[i].x;
y = v[i].y;
if(x<xmin) xmin = x;
if(x>xmax) xmax = x;
if(y<ymin) ymin = y;
if(y>ymax) 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);
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);
pen("solid");
if(bflag) {
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;i<ntrack;i++) {
longlines = LONGLINES;
satellite(&track[i]);
longlines = shortlines;
}
pen("solid");
for(i=0;file[i];i++)
getdata(file[i]);
move(right,bottom);
if(!s1flag)
closepl();
return 0;
}
int
normproj(double lat, double lon, double *x, double *y)
{
int i;
struct place geog;
latlon(lat,lon,&geog);
/*
printp(&geog);
*/
normalize(&geog);
if(!inwindow(&geog))
return(-1);
i = (*projection)(&geog,x,y);
if(rflag)
*x = -*x;
/*
printp(&geog);
fprintf(stderr,"%d %.3f %.3f\n",i,*x,*y);
*/
return(i);
}
int
posproj(double lat, double lon, double *x, double *y)
{
int i;
struct place geog;
latlon(lat,lon,&geog);
normalize(&geog);
i = (*projection)(&geog,x,y);
if(rflag)
*x = -*x;
return(i);
}
int
inwindow(struct place *geog)
{
if(geog->nlat.l<rwindow[0]||
geog->nlat.l>rwindow[1]||
geog->wlon.l<rwindow[2]||
geog->wlon.l>rwindow[3])
return(0);
else return(1);
}
int
inlimits(struct place *g)
{
if(rlimits[0]>g->nlat.l||
rlimits[1]<g->nlat.l)
return(0);
switch(limcase) {
case 0:
if(rlimits[2]+TWOPI>g->wlon.l&&
rlimits[3]<g->wlon.l)
return(0);
break;
case 1:
if(rlimits[2]>g->wlon.l||
rlimits[3]<g->wlon.l)
return(0);
break;
case 2:
if(rlimits[2]>g->wlon.l&&
rlimits[3]-TWOPI<g->wlon.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<hilat;lat+=10.)
for(lon=lolon;lon<hilon;lon+=10.) {
if(!seeable(lat,lon))
continue;
i = pnorm(lat);
j = pnorm(lon);
if((b=patch[i+9][j+18])&1)
continue;
fseek(ifile,b,0);
while((ip=getc(ifile))>=0&&(jp=getc(ifile))>=0){
if(ip!=(i&0377)||jp!=(j&0377))
break;
n = getshort(ifile);
conn = 0;
if(n > 0) { /* absolute coordinates */
for(k=0;k<n;k++){
kx = SCALERATIO*getshort(ifile);
ky = SCALERATIO*getshort(ifile);
if (((k%delta) != 0) && (k != (n-1)))
continue;
conv(kx,&g.nlat);
conv(ky,&g.wlon);
conn = plotpt(&g,conn);
}
} else { /* differential, scaled by SCALERATI0 */
n = -n;
kx = SCALERATIO*getshort(ifile);
ky = SCALERATIO*getshort(ifile);
for(k=0; k<n; k++) {
c = getc(ifile);
if(c&0200) c|= ~0177;
kx += c;
c = getc(ifile);
if(c&0200) c|= ~0177;
ky += c;
if(k%delta!=0&&k!=n-1)
continue;
conv(kx,&g.nlat);
conv(ky,&g.wlon);
conn = plotpt(&g,conn);
}
}
if(k==1) {
conv(kx,&g.nlat);
conv(ky,&g.wlon);
conn = plotpt(&g,conn);
}
}
}
fclose(ifile);
}
int
seeable(double lat0, double lon0)
{
double x, y;
double lat, lon;
for(lat=lat0;lat<=lat0+10;lat+=2*grid[2])
for(lon=lon0;lon<=lon0+10;lon+=2*grid[2])
if(normproj(lat,lon,&x,&y)>0)
return(1);
return(0);
}
void
satellite(struct track *t)
{
char sym[50];
char lbl[50];
double scale;
register conn;
double lat,lon;
struct place place;
static FILE *ifile = stdin;
if(t->tracknam[0]!='-'||t->tracknam[1]!=0) {
fclose(ifile);
if((ifile=fopen(t->tracknam,"r"))==NULL)
filerror("can't find track", t->tracknam);
}
pen(t->tracktyp=='t'?"dotdash":"solid");
for(;;) {
conn = 0;
while(!feof(ifile) && fscanf(ifile,"%f%f",&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 %f",sym,&scale) <= 1)
scale = 1;
if(plotpt(&place,conn?conn:-1)) {
int r = *lbl=='!'?0:rflag?-1:1;
if(putsym(&place,sym,scale,r) == 0)
text(lbl);
}
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);
exit(1);
}
void
filerror(char *s, char *f)
{
fprintf(stderr,"\r\n%s %s\n",s,f);
exit(1);
}
char *
mapindex(char *s)
{
char *t = malloc(strlen(s)+3);
strcpy(t,s);
strcat(t,".x");
return t;
}
#define NOPT 32767
static ox = NOPT, oy = NOPT;
int
cpoint(int xi, int yi, int conn)
{
int dx = abs(ox-xi);
int dy = abs(oy-yi);
if(xi<left||xi>=right || yi<bottom||yi>=top) {
ox = oy = NOPT;
return 0;
}
if(conn == -1) /* isolated plotting symbol */
;
else if(!conn)
move(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;
if(doproj(g,&kx,&ky)<=0) {
return(0);
}
ret = cpoint(kx,ky,conn);
return ret;
}
}
int
doproj(struct place *g, int *kx, int *ky)
{
double x,y,x1,y1;
/*fprintf(stderr,"dopr1 %f %f \n",g->nlat.l,g->wlon.l);*/
if((*projection)(g,&x,&y)<=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(1);
}
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)<=0 || doproj(&oldg,&okx,&oky)<=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)
{
#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))<FUZZ)
d1 = diddle(g1,lon,d2);
if((f2=fabs(d2))<FUZZ) {
d2 = diddle(g2,lon,d1);
if(doproj(g2,&kx,&ky)>0)
cpoint(kx,ky,0);
}
if(f1<FUZZ&&f2<FUZZ)
return(2);
if(f1>PI*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(lat0<tetrapt&&lat0+slat>tetrapt) {
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 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;
}
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<limits[1])
limits[1] = center + colat;
}
short
getshort(FILE *f)
{
int c, r;
c = getc(f);
r = (c | getc(f)<<8);
if (r&0x8000)
r |= ~0xFFFF; /* in case short > 16 bits */
return r;
}
double
fmin(double x, double y)
{
return(x<y?x:y);
}
double
fmax(double x, double y)
{
return(x>y?x:y);
}
void
clamp(double *px, double v)
{
*px = (v<0?fmax:fmin)(*px,v);
}
void
pathnames(void)
{
int i;
char *t, *indexfile;
FILE *f, *fx;
for(i=0; i<NFILE && file[i]; i++) {
if(*file[i]=='/')
continue;
indexfile = mapindex(file[i]);
/* ansi equiv of unix access() call */
f = fopen(file[i], "r");
fx = fopen(indexfile, "r");
if(f) fclose(f);
if(fx) fclose(fx);
free(indexfile);
if(f && fx)
continue;
t = malloc(strlen(file[i])+strlen(mapdir)+2);
strcpy(t,mapdir);
strcat(t,"/");
strcat(t,file[i]);
file[i] = t;
}
}
void
clipinit(void)
{
register i;
double s,t;
if(nvert<=0)
return;
for(i=0; i<nvert; i++) { /*convert latlon to xy*/
if(normproj(v[i].x,v[i].y,&v[i].x,&v[i].y)<=0)
error("invisible clipping vertex");
}
if(nvert==2) { /*rectangle with diag specified*/
nvert = 4;
v[2] = v[1];
v[1].x=v[0].x, v[1].y=v[2].y, v[3].x=v[2].x, v[3].y=v[0].y;
}
v[nvert] = v[0];
v[nvert+1] = v[1];
s = 0;
for(i=1; i<=nvert; i++) { /*test for convexity*/
t = (v[i-1].x-v[i].x)*(v[i+1].y-v[i].y) -
(v[i-1].y-v[i].y)*(v[i+1].x-v[i].x);
if(t<-FUZZ && s>=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; i<nvert; i++) { /*edge equation ax+by=c*/
e[i].a = s*(v[i+1].y - v[i].y);
e[i].b = s*(v[i].x - v[i+1].x);
e[i].c = s*(v[i].x*v[i+1].y - v[i].y*v[i+1].x);
}
}
int
inpoly(double x, double y)
{
register i;
for(i=0; i<nvert; i++) {
register struct edge *ei = &e[i];
double val = x*ei->a + 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 */
}
0707070035051031121006640000030000040000020636350533472776500000600000010246map.h #include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#ifdef PLOT
#include PLOT
#else
#include "iplot.h"
#endif
#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 */
};
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 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);
proj lagrange(void);
proj lambert(double, double);
proj laue(void);
proj loxodromic(double); /* not in library */
proj mecca(double);
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);
proj perspective(double);
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 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 *);
0707070035050601711006640000100000040000011211100533574360100001000000116632map.man
MAP(7) MAP(7)
NAME
map - draw maps on various projections
SYNOPSIS
map _p_r_o_j_e_c_t_i_o_n [ _p_a_r_a_m ... ] [ _o_p_t_i_o_n ... ]
DESCRIPTION
_M_a_p prepares on the standard output a map suitable for dis-
play by any plotting filter described in _p_l_o_t(1). A menu of
projections is produced in response to an unknown
_p_r_o_j_e_c_t_i_o_n. For the meanings of _p_a_r_a_m_s pertinent to partic-
ular projections see _p_r_o_j(3).
The default data for _m_a_p are world shorelines. Option -f
accesses the higher-resolution World Data Bank II.
-f [ _f_e_a_t_u_r_e ... ]
Features are ranked 1 (default) to 4 from major to
minor. Higher-numbered ranks include all lower-
numbered ones. Features are
shore[1-4] seacoasts, lakes, and islands; in the
absence of -m, option -f automatically
includes shore1
ilake[1-2] intermittent lakes
river[1-4] rivers
iriver[1-3] intermittent rivers
canal[1-3] 3=irrigation canals
glacier
iceshelf[12]
reef
saltpan[12]
country[1-3] 2=disputed boundaries, 3=indefinite
boundaries
state states and provinces (US and Canada only)
In other options coordinates are in degrees, with north lat-
itude and west longitude counted as positive.
-l _S _N _E _W
Set the southern and northern latitude and the eastern
and western longitude limits. Missing arguments are
filled out from the list -90, 90, -180, 180.
-k _S _N _E _W
Set the scale as if for a map with limits -l _S _N _E _W
and no -w option.
-o _l_a_t _l_o_n _r_o_t
Orient the map in a nonstandard position. Imagine a
transparent gridded sphere around the globe. Turn the
overlay about the North Pole so that the Prime Meridian
Page 1 Tenth Edition (printed 2/9/93)
MAP(7) MAP(7)
(longitude 0) of the overlay coincides with meridian
_l_o_n on the globe. Then tilt the North Pole of the
overlay along its Prime Meridian to latitude _l_a_t on the
globe. Finally again turn the overlay about its `North
Pole' so that its Prime Meridian coincides with the
previous position of meridian _r_o_t. Project the map in
the standard form appropriate to the overlay, but pre-
senting information from the underlying globe. Missing
arguments are filled out from the list 90, 0, 0. In
the absence of -o, the orientation is 90, 0, _m, where _m
is the middle of the longitude range.
-w _S _N _E _W
Window the map by the specified latitudes and longi-
tudes in the tilted, rotated coordinate system. Miss-
ing arguments are filled out from the list -90, 90,
-180, 180. (It is wise to give an encompassing -l
option with -w. Otherwise for small windows computing
time varies inversely with area!)
-d _n For speed, plot only every _nth point.
-r Reverse left and right (good for star charts and
inside-out views).
-s1
-s2 Superpose. Outputs for a -s1 map (no closing) and a
-s2 map (no opening) may be concatenated.
-g _d_l_a_t _d_l_o_n _r_e_s
Grid spacings are _d_l_a_t, _d_l_o_n. Zero spacing means no
grid. Missing _d_l_a_t is taken to be zero. Missing _d_l_o_n
is taken the same as _d_l_a_t. Grid lines are drawn to a
resolution of _r_e_s (28o9 or less by default). In the
absence of -g, grid spacing is 108o9.
-p _l_a_t _l_o_n _e_x_t_e_n_t
Position the point _l_a_t, _l_o_n at the center of a square
plotting area. Scale the map so that a side of the
square is _e_x_t_e_n_t times the size of one degree of lati-
tude at the center. By default maps are scaled and
positioned to fit within the plotting area. An _e_x_t_e_n_t
overrides option -k.
-c _x _y _r_o_t
After all other positioning and scaling operations,
rotate the image _r_o_t degrees counterclockwise about the
center and move the center to position _x, _y, of the
plotting area, whose nominal extent is -1_<x_<1, -1_<y_<1.
The map is clipped to this area. Missing arguments are
taken to be 0.
-m [ _f_i_l_e ... ]
Page 2 Tenth Edition (printed 2/9/93)
MAP(7) MAP(7)
Use map data from named files. If no files are named,
omit map data. Files that cannot be found directly are
looked up a standard directory, which contains, in
addition to the data for -f,
world World Data Bank I from CIA (default)
states US map from Census Bureau
counties US map from Census Bureau
The environment variables MAP and MAPDIR change the
default map and default directory.
-b [ _l_a_t_1 _l_o_n_1 _l_a_t_2 _l_o_n_2 ... ]
Suppress the drawing of the normal boundary (defined by
options -l and -w). Coordinates, if present, define
the vertices of a polygon to which the map is clipped.
If only two vertices are given, they are taken to be
the diagonal of a rectangle. To draw the polygon, give
its vertices as a -u track.
-t _f_i_l_e ...
The arguments name ASCII files that contain lists of
points, given as latitude-longitude pairs in degrees.
If the first file is named `-', the standard input is
taken instead. The points of each list are plotted as
connected `tracks'.
Points in a track file may be followed by label
strings. A label breaks the track. A label may be
prefixed by ", `:', or `!' and is terminated by a new-
line. An unprefixed string or a string prefixed with "
is displayed at the designated point. The first word
of a `:' or `!' string names a special symbol (see
option -y). An optional numerical second word is a
scale factor for the size of the symbol, 1 by default.
A `:' symbol is aligned with its top to the north; a
`!' symbol is aligned vertically on the page.
-u _f_i_l_e ...
Same as -t, except the tracks are unbroken lines. (-t
tracks are dot-dash lines.)
-y _f_i_l_e
The _f_i_l_e contains _p_l_o_t(5)-style data for `:' or `!'
labels in -t or -u files. Each symbol is defined by a
comment :_n_a_m_e then a sequence of `m' and `v' commands.
Coordinates (0,0) fall on the plotting point. Default
scaling is as if the nominal plotting range were `ra -1
-1 1 1'; `ra' commands in _f_i_l_e change the scaling.
EXAMPLES
map perspective 1.025 -o 40.75 74
Page 3 Tenth Edition (printed 2/9/93)
MAP(7) MAP(7)
A view looking down on New York from 100 miles (0.025
of the 4000-mile earth radius). The job can be done
faster by limiting the map so as not to `plot' the
invisible part of the world: `map perspective 1.025 -o
40.75 74 -l 20 60 30 100'. A circular border can be
forced by adding option `-w 77.33'. (Latitude 77.338o9
falls just inside a polar cap of opening angle arc-
cos(1/1.025) = 12.68048o9.)
map mercator -o 49.25 -106 180
A map whose `equator' is a great circle pasing east-
west through New York. The pole of the map is placed
908o9 away (40.75+49.25=90) on the other side of the
earth. A 1808o9 twist around the pole of the map
arranges that the Prime Meridian of the map runs from
the pole of the map over the North Pole to New York
instead of down the back side of the earth. The same
effect can be had from map mercator -o 130.75 74
map albers 28 45 -l 20 50 60 130 -m states
A customary curved-latitude map of the United States.
map albers 28 45 -l 20 50 60 130 -y yfile -t tfile
An example of tracks, labels, and symbols. Arrows at
New York and Miami are 8% and 12% as long as the map is
wide. The contents of `yfile' and `tfile' are
ra -50 -50 50 50 25.77 80.20 :arrow 12
:arrow 25.77 80.20 Miami
m -1 0 25.77 80.20
v 0 0 35.00 74.02
v -.6 .3 40.67 74.02 !arrow 8
m -.6 -.3 40.67 74.02 " New York
v 0 0 34.05 118.25 Los Angeles
map harrison 2 30 -l -90 90 120 240 -o 90 0 0
A fan view covering 608o9 on either
side of the Date Line, as seen from one earth radius
above the North Pole gazing at the
earth's limb, which is 308o9 off vertical.
Option
-o
overrides the default
-o 90 0 180,
which would rotate
the scene to behind the observer.
FILES
All files in directory $MAPDIR
[1-4]?? World Data Bank II for option -f
world,states,counties
default and other maps for option -m
Page 4 Tenth Edition (printed 2/9/93)
MAP(7) MAP(7)
*.x map indexes
map the program proper
SEE ALSO
_m_a_p(5), _p_r_o_j(3), _p_l_o_t(1)
DIAGNOSTICS
`Map seems to be empty'-a coarse survey found zero extent
within the -l and -w bounds; for maps of limited extent the
grid resolution, _r_e_s, or the limits may have to be refined.
BUGS
The syntax of range specifications in -y files differs from
that in options.
Windows (option -w) cannot cross the Date Line.
No borders appear along edges arising from visibility lim-
its.
Segments that cross a border are dropped, not clipped.
Certain very long line segments are dropped on the assump-
tion that they were intended to go the other way around the
world.
Automatic scaling may miss the extreme points of peculiarly
shaped maps; use option -p to recover.
Although _m_a_p draws grid lines dotted and -t tracks dot-
dashed, many plotting filters cannot cope and make them
solid.
Page 5 Tenth Edition (printed 2/9/93)
PROJ(3X) (bowell) PROJ(3X)
NAME
orient, normalize - map projections
SYNOPSIS
orient(lat, lon, rot)
double lat, lon, rot;
normalize(p)
struct place *p;
DESCRIPTION
Users of _m_a_p(7) may skip to the description of `Projection
generators' below.
The functions _o_r_i_e_n_t and _n_o_r_m_a_l_i_z_e plus a collection of map
projection generators are loaded by option -lmap of _l_d(1).
Most of them calculate maps for a spherical earth. Each map
projection is available in one standard form, into which
data must be normalized for transverse or nonpolar projec-
tions.
Each standard projection is displayed with the Prime Merid-
ian (longitude 0) being a straight vertical line, along
which North is up. The orientation of nonstandard projec-
tions is specified by _o_r_i_e_n_t. Imagine a transparent gridded
sphere around the globe. First turn the overlay about the
North Pole so that the Prime Meridian (longitude 0) of the
overlay coincides with meridian _l_o_n on the globe. Then tilt
the North Pole of the overlay along its Prime Meridian to
latitude _l_a_t on the globe. Finally again turn the overlay
about its `North Pole' so that its Prime Meridian coincides
with the previous position of (the overlay's) meridian _r_o_t.
Project the desired map in the standard form appropriate to
the overlay, but presenting information from the underlying
globe. It is not useful to use _o_r_i_e_n_t without using
_n_o_r_m_a_l_i_z_e.
_N_o_r_m_a_l_i_z_e converts latitude-longitude coordinates on the
globe to coordinates on the overlaid grid. The coordinates
and their sines and cosines are input to _n_o_r_m_a_l_i_z_e in a
place structure. Transformed coordinates and their sines
and cosines are returned in the same structure.
struct place {
double radianlat, sinlat, coslat;
double radianlon, sinlon, coslon;
};
The projection generators return a pointer to a function
that converts normalized coordinates to _x-_y coordinates for
the desired map, or 0 if the required projection is not
available. The returned function is exemplified by _p_r_o_j in
Page 6 Tenth Edition (printed 2/9/93)
PROJ(3X) (bowell) PROJ(3X)
this example:
struct place pt;
int (*proj)() = mercator();
double x, y;
orient(45.0, 30.0, 180.0); /* set coordinate rotation */
. . . /* fill in the pt structure */
normalize(&pt); /* rotate coordinates */
if((*proj)(&pt, &x, &y) > 0) /* project onto x,y plane */
plot(x, y);
The projection function (*proj)() returns 1 for a good
point, 0 for a point on a wrong sheet (e.g. the back of the
world in a perspective projection), and -1 for a point that
is deemed unplottable (e.g. points near the poles on a Mer-
cator projection).
Scaling may be determined from the _x-_y coordinates of
selected points. Latitudes and longitudes are measured in
degrees for ease of specification for _o_r_i_e_n_t and the projec-
tion generators but in radians for ease of calculation for
_n_o_r_m_a_l_i_z_e and _p_r_o_j. In either case latitude is measured pos-
itive north of the equator, and longitude positive west of
Greenwich. Radian longitude should be limited to the range
-_p_i <= _l_o_n < _p_i.
Projection generators
Equatorial projections centered on the Prime Meridian (lon-
gitude 0). Parallels are straight horizontal lines.
mercator() equally spaced straight meridians, confor-
mal, straight compass courses
sinusoidal() equally spaced parallels, equal-area, same
as _b_o_n_n_e(_0)
cylequalarea(lat0) equally spaced straight meridians,
equal-area, true scale on _l_a_t_0
cylindrical() central projection on tangent cylinder
rectangular(lat0) equally spaced parallels, equally
spaced straight meridians, true scale on _l_a_t_0
gall(lat0) parallels spaced stereographically on prime
meridian, equally spaced straight meridians, true scale
on _l_a_t_0
mollweide() (homalographic) equal-area, hemisphere is a
circle
Azimuthal projections centered on the North Pole. Parallels
are concentric circles. Meridians are equally spaced radial
lines.
azequidistant() equally spaced parallels, true dis-
tances from pole
azequalarea() equal-area
Page 7 Tenth Edition (printed 2/9/93)
PROJ(3X) (bowell) PROJ(3X)
gnomonic() central projection on tangent plane,
straight great circles
perspective(dist) viewed along earth's axis _d_i_s_t earth
radii from center of earth
orthographic() viewed from infinity
stereographic() conformal, projected from opposite pole
laue() _r_a_d_i_u_s = tan(2x_c_o_l_a_t_i_t_u_d_e), used in xray crys-
tallography
fisheye(n) stereographic seen from just inside medium
with refractive index n
newyorker(r) _r_a_d_i_u_s = log(_c_o_l_a_t_i_t_u_d_e/_r): extreme `fish-
eye' view from pedestal of radius _r degrees
Polar conic projections symmetric about the Prime Meridian.
Parallels are segments of concentric circles. Except in the
Bonne projection, meridians are equally spaced radial lines
orthogonal to the parallels.
conic(lat0) central projection on cone tangent at _l_a_t_0
simpleconic(lat0,lat1) equally spaced parallels, true
scale on _l_a_t_0 and _l_a_t_1
lambert(lat0,lat1) conformal, true scale on _l_a_t_0 and
_l_a_t_1
albers(lat0,lat1) equal-area, true scale on _l_a_t_0 and
_l_a_t_1
bonne(lat0) equally spaced parallels, equal-area, par-
allel _l_a_t_0 developed from tangent cone
Projections with bilateral symmetry about the Prime Meridian
and the equator.
polyconic() parallels developed from tangent cones,
equally spaced along Prime Meridian
aitoff() equal-area projection of globe onto 2-to-1
ellipse, based on _a_z_e_q_u_a_l_a_r_e_a
lagrange() conformal, maps whole sphere into a circle
bicentric(lon0) points plotted at true azimuth from two
centers on the equator at longitudes +__l_o_n_0, great cir-
cles are straight lines (a stretched gnomonic projec-
tion)
elliptic(lon0) points are plotted at true distance from
two centers on the equator at longitudes +__l_o_n_0
globular() hemisphere is circle, circular arc meridians
equally spaced on equator, circular arc parallels
equally spaced on 0- and 90-degree meridians
vandergrinten() sphere is circle, meridians as in
_g_l_o_b_u_l_a_r, circular arc parallels resemble _m_e_r_c_a_t_o_r
Doubly periodic conformal projections.
guyou() W and E hemispheres are square
square() world is square with Poles at diagonally oppo-
site corners
tetra() map on tetrahedron with edge tangent to Prime
Meridian at S Pole, unfolded into equilateral triangle
Page 8 Tenth Edition (printed 2/9/93)
PROJ(3X) (bowell) PROJ(3X)
hex() world is hexagon centered on N Pole, N and S
hemispheres are equilateral triangles
Miscellaneous projections.
harrison(dist,angle) oblique perspective from above the
North Pole, _d_i_s_t earth radii from center of earth,
looking along the Date Line _a_n_g_l_e degrees off vertical
trapezoidal(lat0,lat1) equally spaced parallels,
straight meridians equally spaced along parallels, true
scale at _l_a_t_0 and _l_a_t_1 on Prime Meridian
Retroazimuthal projections. At every point the angle
between vertical and a straight line to `Mecca', latitude
_l_a_t_0 on the prime meridian, is the true bearing of Mecca.
mecca(lat0) equally spaced vertical meridians
homing(lat0) distances to `Mecca' are true
Maps based on the spheroid. Of geodetic quality, these pro-
jections do not make sense for tilted orientations. For
descriptions, see corresponding maps above.
sp_mercator()
sp_albers(lat0,lat1)
SEE ALSO
_m_a_p(7), _m_a_p(5), _p_l_o_t(3)
BUGS
Only one projection and one orientation can be active at a
time.
The west-longitude-positive convention betrays Yankee chau-
vinism.
Page 9 Tenth Edition (printed 2/9/93)
MAP(5) MAP(5)
NAME
map - digitized map formats
DESCRIPTION
Files used by _m_a_p(7) are a sequence of structures of the
form:
struct {
signed char patchlatitude;
signed char patchlongitude;
short n;
union {
struct {
short latitude;
short longitude;
} point[n];
struct {
short latitude;
short longitude;
struct {
signed char latdiff;
signed char londiff;
} point[-n];
} highres;
} segment;
};
Fields `patchlatitude' and `patchlongitude' tell to what
10-degree by 10-degree patch of the earth's surface a seg-
ment belongs. Their values range from -9 to 8 and from -18
to 17, respectively, and indicate the coordinates of the
southeast corner of the patch in units of 10 degrees.
Each segment of |n| points is connected; consecutive seg-
ments are not necessarily related. Latitude and longitude
are measured in units of 0.0001 radian. If n is negative,
then differences to the first and succeeding points are mea-
sured in units of 0.00001 radian. Latitude is counted posi-
tive to the north and longitude positive to the west.
The patches are ordered lexicographically by `patchlatitude'
then `patchlongitude'. A printable index to the first seg-
ment of each patch in a file named _d_a_t_a is kept in an asso-
ciated file named _d_a_t_a.x. Each line of an index file con-
tains `patchlatitude,' `patchlongitude' and the byte posi-
tion of the patch in the map file. Both the map file and
the index file are ordered by patch latitude and longitude.
Shorts are stored in little-endian order, low byte first,
regardless of computer architecture. To assure portability,
_m_a_p accesses them bytewise.
Page 10 Tenth Edition (printed 2/9/93)
MAP(5) MAP(5)
SEE ALSO
_m_a_p(7), _p_r_o_j(3)
Page 11 Tenth Edition (printed 2/9/93)
PLOT(5) PLOT(5)
NAME
plot - graphics interface
DESCRIPTION
Files of this format are produced by routines described in
_p_l_o_t(3), and are interpreted for various devices by commands
described in _p_l_o_t(1). A graphics file is an ASCII stream of
instruction lines. Arguments are delimited by spaces, tabs,
or commas. Numbers may be floating point. Punctuation
marks (except `:') , spaces, and tabs at the beginning of
lines are ignored. Comments run from `:' to newline. Extra
letters appended to a valid instruction are ignored. Thus
`...line', `line', `li' all mean the same thing. Arguments
are interpreted as follows:
1. If an instruction requires no arguments, the rest of
the line is ignored.
2. If it requires a string argument, then all the line
after the first field separator is passed as argument.
Quote marks may be used to preserve leading blanks.
Strings may include newlines represented as `\n'.
3. Between numeric arguments alphabetic characters and
punctuation marks are ignored. Thus line from 5 6 to 7
8 draws a line from (5, 6) to (7, 8).
4. Instructions with numeric arguments remain in effect
until a new instruction is read. Such commands may
spill over many lines. Thus the following sequence will
draw a polygon with vertices (4.5, 6.77), (5.8, 5.6),
(7.8, 4.55), and (10.0, 3.6).
move 4.5 6.77
vec 5.8, 5.6 7.8
4.55 10.0, 3.6 4.5, 6.77
The instructions are executed in order. The last designated
point in a line, move, rmove, vec, rvec, arc, or point com-
mand becomes the `current point' (_X,_Y) for the next command.
Each of the following descriptions corresponds to a routine
in _p_l_o_t(3).
Open & Close
o _s_t_r_i_n_g Open plotting device. For _t_r_o_f_f, _s_t_r_i_n_g specifies
the size of the plot (default is `6i.')
cl Close plotting device.
Basic Plotting Commands
e Start another frame of output or erase the screen
on CRT terminals without scroll.
m _x _y (move) Current point becomes _x _y.
Page 12 Tenth Edition (printed 2/9/93)
PLOT(5) PLOT(5)
rm _d_x _d_y Current point becomes _X+_d_x _Y+_d_y.
poi _x _y Plot the point _x _y and make it the current point.
v _x _y Draw a vector from the current point to _x _y.
rv _d_x _d_y Draw vector from current point to X+dx Y+dy
li _x_1 _y_1 _x_2 _y_2
Draw a line from _x_1 _y_1 to _x_2 _y_2. Make the current
point _x_2 _y_2.
t _s_t_r_i_n_g Place the _s_t_r_i_n_g so that its first character is
centered on the current point (default). If
_s_t_r_i_n_g begins with `\C' (`\R'), it is centered
(right-adjusted) on the current point. A back-
slash at the beginning of the string may be
escaped with another backslash.
a _x_1 _y_1 _x_2 _y_2 _x_c _y_c _r
Draw a circular arc from _x_1 _y_1 to _x_2 _y_2 with cen-
ter _x_c _y_c and radius _r. If the radius is positive,
the arc is drawn counterclockwise; negative,
clockwise. The starting point is exact but the
ending point is approximate.
ci _x_c _y_c _r
Draw a circle centered at _x_c _y_c with radius _r. If
the range and frame parameters do not specify a
square, the `circle' will be elliptical.
di _x_c _y_c _r
Draw a disc centered at _x_c _y_c with radius _r using
the filling color (see cfill below). Only works
on the 5620; on other devices is the same as
circle.
bo _x_1 _y_1 _x_2 _y_2
Draw a box with lower left corner at _x_1 _y_1 and
upper right corner at _x_2 _y_2.
sb _x_1 _y_1 _x_2 _y_2
Draw a solid box with lower left corner at _x_1 _y_1
and upper right corner at _x_2 _y_2 using the filling
color (see cfill below).
par _x_1 _y_1 _x_2 _y_2 _x_g _y_g
Draw a parabola from _x_1 _y_1 to _x_2 _y_2 `guided' by _x_g
_y_g. The parabola passes through the midpoint of
the line joining _x_g _y_g with the midpoint of the
line joining _x_1 _y_1 and _x_2 _y_2 and is tangent to the
lines from _x_g _y_g to the endpoints.
pol { {_x_1 _y_1 ... _x_n _y_n} ... {_X_1 _Y_1 ... _X_m _Y_m} }
Draw polygons with vertices _x_1 _y_1 ... _x_n _y_n and _X_1
_Y_1 ... _X_m _Y_m. If only one polygon is specified,
the inner brackets are not needed.
fi { {_x_1 _y_1 ... _x_n _y_n} ... {_X_1 _Y_1 ... _X_m _Y_m} }
Fill a polygon. The arguments are the same as
those for pol except that the first vertex is
automatically repeated to close each polygon. The
polygons do not have to be connected. Enclosed
polygons appear as holes.
sp { {_x_1 _y_1 ... _x_n _y_n} ... {_X_1 _Y_1 ... _X_m _Y_m} }
Page 13 Tenth Edition (printed 2/9/93)
PLOT(5) PLOT(5)
Draw a parabolic spline guided by _x_1 _y_1 ... _x_n _y_n
with simple endpoints.
fsp { {_x_1 _y_1 ... _x_n _y_n} ... {_X_1 _Y_1 ... _X_m _Y_m} }
Draw a parabolic spline guided by _x_1 _y_1 ... _x_n _y_n
with double first endpoint.
lsp { {_x_1 _y_1 ... _x_n _y_n} ... {_X_1 _Y_1 ... _X_m _Y_m} }
Draw a parabolic spline guided by _x_1 _y_1 ... _x_n _y_n
with double last endpoint.
dsp { {_x_1 _y_1 ... _x_n _y_n} ... {_X_1 _Y_1 ... _X_m _Y_m} }
Draw a parabolic spline guided by _x_1 _y_1 ... _x_n _y_n
with double endpoints.
csp { {_x_1 _y_1 ... _x_n _y_n} ... {_X_1 _Y_1 ... _X_m _Y_m} }
in _f_i_l_e_n_a_m_e
(include) Take commands from _f_i_l_e_n_a_m_e.
de _s_t_r_i_n_g { _c_o_m_m_a_n_d_s }
Define _s_t_r_i_n_g as _c_o_m_m_a_n_d_s.
ca _s_t_r_i_n_g _s_c_a_l_e
Invoke commands defined as _s_t_r_i_n_g applying _s_c_a_l_e
to all coordinates.
Commands Controlling the Environment
co _s_t_r_i_n_g Draw lines with color _s_t_r_i_n_g. Available colors
depend on the device. _S_t_r_i_n_g may contain defini-
tions for several devices separated by `/'. Col-
ors possible for the various devices are:
pen black, red, green, blue, Tblack, Tred,
Tgreen, Tblue (assumes default carousel,
T=thick)
1-8 (pen number)
S_n_u_m_b_e_r character size as % of plotting area
troff
F_f_o_n_t
P_p_o_i_n_t _s_i_z_e
2621 H_c_h_a_r_a_c_t_e_r for plotting
pe _s_t_r_i_n_g Use _s_t_r_i_n_g as the style for drawing lines. Not
all pen styles are implemented for all devices.
_S_t_r_i_n_g may contain definitions for several devices
separated by `/'. The available pen styles are:
pen solid, dott[ed], short, long, dotd[ashed],
cdash, ddash
4014 solid , dott[ed], short, long, dotd[ashed],
ddash
troff
solid, dash only straight lines will be
dashed
5620 B_n_u_m_b_e_r line thickness
2621 H_c_h_a_r_a_c_t_e_r for plotting
cf _s_t_r_i_n_g Color for filling; may contain the definitions for
several devices. separated by `/'. The following
Page 14 Tenth Edition (printed 2/9/93)
PLOT(5) PLOT(5)
colors are available on the specified devices:
pen black, red, green, blue, Tblack, Tred,
Tgreen, Tblue
1-8 pen number
5620 B_t_e_x_t_u_r_e string with octal numbers for tex-
ture; see _t_y_p_e_s(9.5). The 16 words of tex-
ture should be followed by one word for the
mode used by _t_e_x_t_u_r_e(); see _b_i_t_b_l_t(9.3).
2621 H_c_h_a_r_a_c_t_e_r for filling
All /A_d_e_g_r_e_e_s slant of shading lines
/G_n_u_m_b_e_r gap between shading lines (in user
units)
ra _x_1 _y_1 _x_2 _y_2
The data will fall between _x_1 _y_1 and _x_2 _y_2. The
plot will be magnified or reduced to fit the
device as closely as possible.
Range settings that exactly fill the plotting area
with unity scaling appear below for devices sup-
ported by the filters of _p_l_o_t(1). The upper limit
is just outside the plotting area. In every case
the plotting area is taken to be square; points
outside may be displayable on devices with non-
square faces.
4014 `range 0. 0. 3120. 3120.'
troff `range 0. 0. 6144. 6144.'
2621 `range 0. 0. 22. 22.'
5620 range dependent on layer size
pen range dependent on paper size
fr _p_x_1 _p_y_1 _p_x_2 _p_y_2
Plot the data in the fraction of the display spec-
ified by _p_x_1 _p_y_1 for lower left corner and _p_x_2 _p_y_2
for upper right corner. Thus `frame .5 0 1. .5'
plots in the lower right quadrant of the display;
`frame 0. 1. 1. 0.' uses the whole display but
inverts the _y coordinates.
sa Save the current environment, and move to a new
one. The new environment inherits the old one.
There are 7 levels.
re Restore previous environment.
SEE ALSO
_p_l_o_t(1), _p_l_o_t(3), _g_r_a_p_h(1)
Page 15 Tenth Edition (printed 2/9/93)
ROUTE(1) ROUTE(1)
NAME
route - map orientations and great circle paths
SYNOPSIS
route [ -t ]
DESCRIPTION
_R_o_u_t_e without an option reads from the standard input a
sequence of points, expressed as latitude-longitude pairs.
For each pair of points, say A and B, _r_o_u_t_e prints on the
standard output
- A _m_a_p(7) orientation, expressed as a -o option, that
transforms the great circle from A to B into the equa-
tor. The transformed locations of A and B are located
on the equator at equal distances west and east of the
Prime Meridian.
- The transformed coordinates of the two points.
- The same information for a map orientation in which
point A appears to the east and B to the west of the
Prime Meridian.
Under option -t _r_o_u_t_e produces coordinates along the great
circle from A to B in the form of a track for use with _m_a_p
option -t.
Coordinates are expressed in degrees, with latitude positive
north of the equator and longitude positive west of Green-
wich.
EXAMPLES
echo 40.75 74 52 0 | route -t | map mercator -l 30 70 -10 80
-t - | plot
Produce a map of the North Atlantic with a great circle
plotted from New York to London.
route
40.75 74 52 0
Produces the output
-o 35.8283 -157.5308 -168.6148
A -0.0000 24.9943
B -0.0000 -24.9943
-o -35.8283 22.4692 168.6149
A -0.0000 -24.9943
B 0.0000 24.9943
From this we derive the following command to draw a
strip map 20 degrees wide along the great circle from
New York to London. The -w option windows the map in
Page 16 Tenth Edition (printed 2/9/93)
ROUTE(1) ROUTE(1)
the new coordinate system, where the equator runs along
the great circle and the two cities sit at 25W and 25E
on that `equator'. The -l option avoids surveying the
whole rest of the world for plottable points.
map mercator -o 35.83 -157.53 -168.61 -w -10 10 -30 30
-l 30 75 -10 80
SEE ALSO
_m_a_p(7)
Page 17 Tenth Edition (printed 2/9/93)
0707070035050675411007750000030000040000010322010533573267500000700000002667map.sh # MF -m -f or "", M map files, A other arguments
# MFLAG and FFLAG 0 or 1 for -m or -f ever used
FEATURE=no
MAPDIR=${MAPDIR-/usr/lib/map}
MAPPROG=${MAPPROG-$MAPDIR/map}
M= A= MF= MFLAG=0 FFLAG=0
for i in $*
do
case $i in
-m)
MF=-m MFLAG=1 ;;
-f)
MF=-f FFLAG=1 ;;
-*)
MF= A="$A $i" ;;
*)
case $MF$i in
-m*) M="$M $i" ;;
-friv*4) M="$M 201 202 203 204";;
-friv*3) M="$M 201 202 203";;
-friv*2) M="$M 201 202";;
-friv*) M="$M 201";;
-firiv*[34]) M="$M 206 207 208";;
-firiv*2) M="$M 206 207";;
-firiv*) M="$M 206";;
-fcoast*4|-fshore*4|-flake*4) M="$M 102 103 104";;
-fcoast*3|-fshore*3|-flake*3) M="$M 102 103";;
-fcoast*2|-fshore*2|-flake*2) M="$M 102";;
-fcoast*|-fshore*|-flake*) M="$M 101";;
-filake*[234]|-fishore*[234]) M="$M 106 107";;
-filake*|-fishore*) M="$M 106";;
-freef*) M="$M 108";;
-fcanal*[34]) M="$M 210 211 212";;
-fcanal*2) M="$M 210 211";;
-fcanal*) M="$M 210";;
-fglacier*) M="$M 115";;
-fstate*|-fprovince*) M="$M 401";;
-fcountr*[34]) M="$M 301 302 303";;
-fcountr*2) M="$M 301 302";;
-fcountr*) M="$M 301";;
-fsalt*[234]) M="$M 109 110";;
-fsalt*) M="$M 109";;
-fice*[234]|-fshel*[234]) M="$M 113 114";;
-fice*|-fshel*) M="$M 113";;
-f*) echo map: unknown feature $i 1>&2
exit 1 ;;
*)
A="$A $i" MF=
esac
esac
done
case "$MFLAG$FFLAG$M " in
00)
: ;;
1* | *" 101 "*)
M="-m $M" ;;
01*)
M="-m 101 $M"
esac
MAP=${MAP-world} MAPDIR=$MAPDIR $MAPPROG $A $M
0707070035050415411004640000030000040000011762600330024205300002100001013370mapdata/counties n<:oB:qH:sL: �7�7t7R74�6#�6
�6��6 ��6��6��6��6��6��6��6��6 ?�7E�7P�7U�7Y�7[�7]�7_�7a�7d�7e�7c�7b�7b�7 b�7b�7b8a8h8i8h8d8a8T8O8I8H8I8J
8J8I8M8 M8L8O8V"8^$8 �W8�X8�Z8�Z8�Y8�Y8�Y8�\8 �\8�_8�a8�a8�a8�`8�^8�Z8�W8�U8�T8�V8�Z8�a8 �a8�d8�f8�f8�h8 �|7��7��7��7��7��7��7��7��7��7��7��7��7��7��7��7��7��7��7��7��7 ��6��6��6��6��6��6��6 ��6��6��6��6��6��6 ��6��6��6��6��6��6��6��6{�6y�6w�6t�6n�6k�6g�6 g�6b�6\�6U�6P�6L�6I�6E�6 E�6@�6>�6;�65�62�64�63�6/�6/�64�65�64�6/�6-�6 -�6-�6+7'7'7*7+7+'7 ��8 �8�8�8�8�8�8�8�8�8�8�8�8�8�8!�8$�8*�8-�8 -�83�8;�8 6\7&U7L7 q�9o�9m�9k�9h�9f�9d�9d�9e�9h�9j�9i�9h�9f�9g�9g�9f�9`�9Y�9Q�9M�9J�9J�9K�9P�9R�9T�9 T�9T�9S�9O�9M�9M�9L�9G�9F�9F�9G :G:"