1 #include <u.h>
2 #include <libc.h>
3 #include <stdio.h>
4 #include "map.h"
5 #include "iplot.h"
6 
7 #define NVERT 20	/* max number of vertices in a -v polygon */
8 #define HALFWIDTH 8192	/* output scaled to fit in -HALFWIDTH,HALFWIDTH */
9 #define LONGLINES (HALFWIDTH*4)	/* permissible segment lengths */
10 #define SHORTLINES (HALFWIDTH/8)
11 #define SCALERATIO 10	/* of abs to rel data (see map(5)) */
12 #define RESOL 2.	/* coarsest resolution for tracing grid (degrees) */
13 #define TWO_THRD 0.66666666666666667
14 
15 int normproj(double, double, double *, double *);
16 int posproj(double, double, double *, double *);
17 int picut(struct place *, struct place *, double *);
18 double reduce(double);
19 short getshort(FILE *);
20 char *mapindex(char *);
21 proj projection;
22 
23 
24 static char *mapdir = "#9/map";	/* default map directory */
25 struct file {
26 	char *name;
27 	char *color;
28 	char *style;
29 };
30 static struct file dfltfile = {
31 	"world", BLACK, SOLID	/* default map */
32 };
33 static struct file *file = &dfltfile;	/* list of map files */
34 static int nfile = 1;			/* length of list */
35 static char *currcolor = BLACK;		/* current color */
36 static char *gridcolor = BLACK;
37 static char *bordcolor = BLACK;
38 
39 extern struct index index[];
40 int halfwidth = HALFWIDTH;
41 
42 static int (*cut)(struct place *, struct place *, double *);
43 static int (*limb)(double*, double*, double);
44 static void dolimb(void);
45 static int onlimb;
46 static int poles;
47 static double orientation[3] = { 90., 0., 0. };	/* -o option */
48 static int oriented;	/* nonzero if -o option occurred */
49 static int upright;		/* 1 if orientation[0]==90, -1 if -90, else 0*/
50 static int delta = 1;	/* -d setting */
51 static double limits[4] = {	/* -l parameters */
52 	-90., 90., -180., 180.
53 };
54 static double klimits[4] = {	/* -k parameters */
55 	-90., 90., -180., 180.
56 };
57 static int limcase;
58 static double rlimits[4];	/* limits expressed in radians */
59 static double lolat, hilat, lolon, hilon;
60 static double window[4] = {	/* option -w */
61 	-90., 90., -180., 180.
62 };
63 static int windowed;	/* nozero if option -w */
64 static struct vert { double x, y; } v[NVERT+2];	/*clipping polygon*/
65 static struct edge { double a, b, c; } e[NVERT]; /* coeffs for linear inequality */
66 static int nvert;	/* number of vertices in clipping polygon */
67 
68 static double rwindow[4];	/* window, expressed in radians */
69 static double params[2];		/* projection params */
70 /* bounds on output values before scaling; found by coarse survey */
71 static double xmin = 100.;
72 static double xmax = -100.;
73 static double ymin = 100.;
74 static double ymax = -100.;
75 static double xcent, ycent;
76 static double xoff, yoff;
77 double xrange, yrange;
78 static int left = -HALFWIDTH;
79 static int right = HALFWIDTH;
80 static int bottom = -HALFWIDTH;
81 static int top = HALFWIDTH;
82 static int longlines = SHORTLINES; /* drop longer segments */
83 static int shortlines = SHORTLINES;
84 static int bflag = 1;	/* 0 for option -b */
85 static int s1flag = 0;	/* 1 for option -s1 */
86 static int s2flag = 0;	/* 1 for option -s2 */
87 static int rflag = 0;	/* 1 for option -r */
88 static int kflag = 0;	/* 1 if option -k occurred */
89 static int xflag = 0;	/* 1 for option -x */
90        int vflag = 1;	/* -1 if option -v occurred */
91 static double position[3];	/* option -p */
92 static double center[3] = {0., 0., 0.};	/* option -c */
93 static struct coord crot;		/* option -c */
94 static double grid[3] = { 10., 10., RESOL };	/* option -g */
95 static double dlat, dlon;	/* resolution for tracing grid in lat and lon */
96 static double scaling;	/* to compute final integer output */
97 static struct file *track;	/* options -t and -u */
98 static int ntrack;		/* number of tracks present */
99 static char *symbolfile;	/* option -y */
100 
101 void	clamp(double *px, double v);
102 void	clipinit(void);
103 double	diddle(struct place *, double, double);
104 double	diddle(struct place *, double, double);
105 void	dobounds(double, double, double, double, int);
106 void	dogrid(double, double, double, double);
107 int	duple(struct place *, double);
108 #define fmax map_fmax
109 #define fmin map_fmin
110 double	fmax(double, double);
111 double	fmin(double, double);
112 void	getdata(char *);
113 int	gridpt(double, double, int);
114 int	inpoly(double, double);
115 int	inwindow(struct place *);
116 void	pathnames(void);
117 int	pnorm(double);
118 void	radbds(double *w, double *rw);
119 void	revlon(struct place *, double);
120 void	satellite(struct file *);
121 int	seeable(double, double);
122 void	windlim(void);
123 void	realcut(void);
124 
125 int
126 option(char *s)
127 {
128 
129 	if(s[0]=='-' && (s[1]<'0'||s[1]>'9'))
130 		return(s[1]!='.'&&s[1]!=0);
131 	else
132 		return(0);
133 }
134 
135 void
136 conv(int k, struct coord *g)
137 {
138 	g->l = (0.0001/SCALERATIO)*k;
139 	sincos(g);
140 }
141 
142 int
143 main(int argc, char *argv[])
144 {
145 	int i,k;
146 	char *s, *t, *style;
147 	double x, y;
148 	double lat, lon;
149 	double *wlim;
150 	double dd;
151 	if(sizeof(short)!=2)
152 		abort();	/* getshort() won't work */
153 	mapdir = unsharp(mapdir);
154 	s = getenv("MAP");
155 	if(s)
156 		file[0].name = s;
157 	s = getenv("MAPDIR");
158 	if(s)
159 		mapdir = s;
160 	if(argc<=1)
161 		error("usage: map projection params options");
162 	for(k=0;index[k].name;k++) {
163 		s = index[k].name;
164 		t = argv[1];
165 		while(*s == *t){
166 			if(*s==0) goto found;
167 			s++;
168 			t++;
169 		}
170 	}
171 	fprintf(stderr,"projections:\n");
172 	for(i=0;index[i].name;i++)  {
173 		fprintf(stderr,"%s",index[i].name);
174 		for(k=0; k<index[i].npar; k++)
175 			fprintf(stderr," p%d", k);
176 		fprintf(stderr,"\n");
177 	}
178 	exits("error");
179 found:
180 	argv += 2;
181 	argc -= 2;
182 	cut = index[k].cut;
183 	limb = index[k].limb;
184 	poles = index[k].poles;
185 	for(i=0;i<index[k].npar;i++) {
186 		if(i>=argc||option(argv[i])) {
187 			fprintf(stderr,"%s needs %d params\n",index[k].name,index[k].npar);
188 			exits("error");
189 		}
190 		params[i] = atof(argv[i]);
191 	}
192 	argv += i;
193 	argc -= i;
194 	while(argc>0&&option(argv[0])) {
195 		argc--;
196 		argv++;
197 		switch(argv[-1][1]) {
198 		case 'm':
199 			if(file == &dfltfile) {
200 				file = 0;
201 				nfile = 0;
202 			}
203 			while(argc && !option(*argv)) {
204 				file = realloc(file,(nfile+1)*sizeof(*file));
205 				file[nfile].name = *argv;
206 				file[nfile].color = currcolor;
207 				file[nfile].style = SOLID;
208 				nfile++;
209 				argv++;
210 				argc--;
211 			}
212 			break;
213 		case 'b':
214 			bflag = 0;
215 			for(nvert=0;nvert<NVERT&&argc>=2;nvert++) {
216 				if(option(*argv))
217 					break;
218 				v[nvert].x = atof(*argv++);
219 				argc--;
220 				if(option(*argv))
221 					break;
222 				v[nvert].y = atof(*argv++);
223 				argc--;
224 			}
225 			if(nvert>=NVERT)
226 				error("too many clipping vertices");
227 			break;
228 		case 'g':
229 			gridcolor = currcolor;
230 			for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
231 				grid[i] = atof(argv[i]);
232 			switch(i) {
233 			case 0:
234 				grid[0] = grid[1] = 0.;
235 				break;
236 			case 1:
237 				grid[1] = grid[0];
238 			}
239 			argc -= i;
240 			argv += i;
241 			break;
242 		case 't':
243 			style = SOLID;
244 			goto casetu;
245 		case 'u':
246 			style = DOTDASH;
247 		casetu:
248 			while(argc && !option(*argv)) {
249 				track = realloc(track,(ntrack+1)*sizeof(*track));
250 				track[ntrack].name = *argv;
251 				track[ntrack].color = currcolor;
252 				track[ntrack].style = style;
253 				ntrack++;
254 				argv++;
255 				argc--;
256 			}
257 			break;
258 		case 'r':
259 			rflag++;
260 			break;
261 		case 's':
262 			switch(argv[-1][2]) {
263 			case '1':
264 				s1flag++;
265 				break;
266 			case 0:		/* compatibility */
267 			case '2':
268 				s2flag++;
269 			}
270 			break;
271 		case 'o':
272 			for(i=0;i<3&&i<argc&&!option(argv[i]);i++)
273 				orientation[i] = atof(argv[i]);
274 			oriented++;
275 			argv += i;
276 			argc -= i;
277 			break;
278 		case 'l':
279 			bordcolor = currcolor;
280 			for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
281 				limits[i] = atof(argv[i]);
282 			argv += i;
283 			argc -= i;
284 			break;
285 		case 'k':
286 			kflag++;
287 			for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
288 				klimits[i] = atof(argv[i]);
289 			argv += i;
290 			argc -= i;
291 			break;
292 		case 'd':
293 			if(argc>0&&!option(argv[0])) {
294 				delta = atoi(argv[0]);
295 				argv++;
296 				argc--;
297 			}
298 			break;
299 		case 'w':
300 			bordcolor = currcolor;
301 			windowed++;
302 			for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
303 				window[i] = atof(argv[i]);
304 			argv += i;
305 			argc -= i;
306 			break;
307 		case 'c':
308 			for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
309 				center[i] = atof(argv[i]);
310 			argc -= i;
311 			argv += i;
312 			break;
313 		case 'p':
314 			for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
315 				position[i] = atof(argv[i]);
316 			argc -= i;
317 			argv += i;
318 			if(i!=3||position[2]<=0)
319 				error("incomplete positioning");
320 			break;
321 		case 'y':
322 			if(argc>0&&!option(argv[0])) {
323 				symbolfile = argv[0];
324 				argc--;
325 				argv++;
326 			}
327 			break;
328 		case 'v':
329 			if(index[k].limb == 0)
330 				error("-v does not apply here");
331 			vflag = -1;
332 			break;
333 		case 'x':
334 			xflag = 1;
335 			break;
336 		case 'C':
337 			if(argc && !option(*argv)) {
338 				currcolor = colorcode(*argv);
339 				argc--;
340 				argv++;
341 			}
342 			break;
343 		}
344 	}
345 	if(argc>0)
346 		error("error in arguments");
347 	pathnames();
348 	clamp(&limits[0],-90.);
349 	clamp(&limits[1],90.);
350 	clamp(&klimits[0],-90.);
351 	clamp(&klimits[1],90.);
352 	clamp(&window[0],-90.);
353 	clamp(&window[1],90.);
354 	radbds(limits,rlimits);
355 	limcase = limits[2]<-180.?0:
356 		  limits[3]>180.?2:
357 		  1;
358 	if(
359 		window[0]>=window[1]||
360 		window[2]>=window[3]||
361 		window[0]>90.||
362 		window[1]<-90.||
363 		window[2]>180.||
364 		window[3]<-180.)
365 		error("unreasonable window");
366 	windlim();
367 	radbds(window,rwindow);
368 	upright = orientation[0]==90? 1: orientation[0]==-90? -1: 0;
369 	if(index[k].spheroid && !upright)
370 		error("can't tilt the spheroid");
371 	if(limits[2]>limits[3])
372 		limits[3] += 360;
373 	if(!oriented)
374 		orientation[2] = (limits[2]+limits[3])/2;
375 	orient(orientation[0],orientation[1],orientation[2]);
376 	projection = (*index[k].prog)(params[0],params[1]);
377 	if(projection == 0)
378 		error("unreasonable projection parameters");
379 	clipinit();
380 	grid[0] = fabs(grid[0]);
381 	grid[1] = fabs(grid[1]);
382 	if(!kflag)
383 		for(i=0;i<4;i++)
384 			klimits[i] = limits[i];
385 	if(klimits[2]>klimits[3])
386 		klimits[3] += 360;
387 	lolat = limits[0];
388 	hilat = limits[1];
389 	lolon = limits[2];
390 	hilon = limits[3];
391 	if(lolon>=hilon||lolat>=hilat||lolat<-90.||hilat>90.)
392 		error("unreasonable limits");
393 	wlim = kflag? klimits: window;
394 	dlat = fmin(hilat-lolat,wlim[1]-wlim[0])/16;
395 	dlon = fmin(hilon-lolon,wlim[3]-wlim[2])/32;
396 	dd = fmax(dlat,dlon);
397 	while(grid[2]>fmin(dlat,dlon)/2)
398 		grid[2] /= 2;
399 	realcut();
400 	if(nvert<=0) {
401 		for(lat=klimits[0];lat<klimits[1]+dd-FUZZ;lat+=dd) {
402 			if(lat>klimits[1])
403 				lat = klimits[1];
404 			for(lon=klimits[2];lon<klimits[3]+dd-FUZZ;lon+=dd) {
405 				i = (kflag?posproj:normproj)
406 					(lat,lon+(lon<klimits[3]?FUZZ:-FUZZ),
407 				   &x,&y);
408 				if(i*vflag <= 0)
409 					continue;
410 				if(x<xmin) xmin = x;
411 				if(x>xmax) xmax = x;
412 				if(y<ymin) ymin = y;
413 				if(y>ymax) ymax = y;
414 			}
415 		}
416 	} else {
417 		for(i=0; i<nvert; i++) {
418 			x = v[i].x;
419 			y = v[i].y;
420 			if(x<xmin) xmin = x;
421 			if(x>xmax) xmax = x;
422 			if(y<ymin) ymin = y;
423 			if(y>ymax) ymax = y;
424 		}
425 	}
426 	xrange = xmax - xmin;
427 	yrange = ymax - ymin;
428 	if(xrange<=0||yrange<=0)
429 		error("map seems to be empty");
430 	scaling = 2;	/*plotting area from -1 to 1*/
431 	if(position[2]!=0) {
432 		if(posproj(position[0]-.5,position[1],&xcent,&ycent)==0||
433 		   posproj(position[0]+.5,position[1],&x,&y)==0)
434 			error("unreasonable position");
435 		scaling /= (position[2]*hypot(x-xcent,y-ycent));
436 		if(posproj(position[0],position[1],&xcent,&ycent)==0)
437 			error("unreasonable position");
438 	} else {
439 		scaling /= (xrange>yrange?xrange:yrange);
440 		xcent = (xmin+xmax)/2;
441 		ycent = (ymin+ymax)/2;
442 	}
443 	xoff = center[0]/scaling;
444 	yoff = center[1]/scaling;
445 	crot.l = center[2]*RAD;
446 	sincos(&crot);
447 	scaling *= HALFWIDTH*0.9;
448 	if(symbolfile)
449 		getsyms(symbolfile);
450 	if(!s2flag) {
451 		openpl();
452 		erase();
453 	}
454 	range(left,bottom,right,top);
455 	comment("grid","");
456 	colorx(gridcolor);
457 	pen(DOTTED);
458 	if(grid[0]>0.)
459 		for(lat=ceil(lolat/grid[0])*grid[0];
460 		    lat<=hilat;lat+=grid[0])
461 			dogrid(lat,lat,lolon,hilon);
462 	if(grid[1]>0.)
463 		for(lon=ceil(lolon/grid[1])*grid[1];
464 		    lon<=hilon;lon+=grid[1])
465 			dogrid(lolat,hilat,lon,lon);
466 	comment("border","");
467 	colorx(bordcolor);
468 	pen(SOLID);
469 	if(bflag) {
470 		dolimb();
471 		dobounds(lolat,hilat,lolon,hilon,0);
472 		dobounds(window[0],window[1],window[2],window[3],1);
473 	}
474 	lolat = floor(limits[0]/10)*10;
475 	hilat = ceil(limits[1]/10)*10;
476 	lolon = floor(limits[2]/10)*10;
477 	hilon = ceil(limits[3]/10)*10;
478 	if(lolon>hilon)
479 		hilon += 360.;
480 	/*do tracks first so as not to lose the standard input*/
481 	for(i=0;i<ntrack;i++) {
482 		longlines = LONGLINES;
483 		satellite(&track[i]);
484 		longlines = shortlines;
485 	}
486 	for(i=0;i<nfile;i++) {
487 		comment("mapfile",file[i].name);
488 		colorx(file[i].color);
489 		pen(file[i].style);
490 		getdata(file[i].name);
491 	}
492 	move(right,bottom);
493 	if(!s1flag)
494 		closepl();
495 	return 0;
496 }
497 
498 /* Out of perverseness (really to recover from a dubious,
499    but documented, convention) the returns from projection
500    functions (-1 unplottable, 0 wrong sheet, 1 good) are
501    recoded into -1 wrong sheet, 0 unplottable, 1 good. */
502 
503 int
504 fixproj(struct place *g, double *x, double *y)
505 {
506 	int i = (*projection)(g,x,y);
507 	return i<0? 0: i==0? -1: 1;
508 }
509 
510 int
511 normproj(double lat, double lon, double *x, double *y)
512 {
513 	int i;
514 	struct place geog;
515 	latlon(lat,lon,&geog);
516 /*
517 	printp(&geog);
518 */
519 	normalize(&geog);
520 	if(!inwindow(&geog))
521 		return(-1);
522 	i = fixproj(&geog,x,y);
523 	if(rflag)
524 		*x = -*x;
525 /*
526 	printp(&geog);
527 	fprintf(stderr,"%d %.3f %.3f\n",i,*x,*y);
528 */
529 	return(i);
530 }
531 
532 int
533 posproj(double lat, double lon, double *x, double *y)
534 {
535 	int i;
536 	struct place geog;
537 	latlon(lat,lon,&geog);
538 	normalize(&geog);
539 	i = fixproj(&geog,x,y);
540 	if(rflag)
541 		*x = -*x;
542 	return(i);
543 }
544 
545 int
546 inwindow(struct place *geog)
547 {
548 	if(geog->nlat.l<rwindow[0]-FUZZ||
549 	   geog->nlat.l>rwindow[1]+FUZZ||
550 	   geog->wlon.l<rwindow[2]-FUZZ||
551 	   geog->wlon.l>rwindow[3]+FUZZ)
552 		return(0);
553 	else return(1);
554 }
555 
556 int
557 inlimits(struct place *g)
558 {
559 	if(rlimits[0]-FUZZ>g->nlat.l||
560 	   rlimits[1]+FUZZ<g->nlat.l)
561 		return(0);
562 	switch(limcase) {
563 	case 0:
564 		if(rlimits[2]+TWOPI-FUZZ>g->wlon.l&&
565 		   rlimits[3]+FUZZ<g->wlon.l)
566 			return(0);
567 		break;
568 	case 1:
569 		if(rlimits[2]-FUZZ>g->wlon.l||
570 		   rlimits[3]+FUZZ<g->wlon.l)
571 			return(0);
572 		break;
573 	case 2:
574 		if(rlimits[2]>g->wlon.l&&
575 		   rlimits[3]-TWOPI+FUZZ<g->wlon.l)
576 			return(0);
577 		break;
578 	}
579 	return(1);
580 }
581 
582 
583 long patch[18][36];
584 
585 void
586 getdata(char *mapfile)
587 {
588 	char *indexfile;
589 	int kx,ky,c;
590 	int k;
591 	long b;
592 	long *p;
593 	int ip, jp;
594 	int n;
595 	struct place g;
596 	int i, j;
597 	double lat, lon;
598 	int conn;
599 	FILE *ifile, *xfile;
600 
601 	indexfile = mapindex(mapfile);
602 	xfile = fopen(indexfile,"r");
603 	if(xfile==NULL)
604 		filerror("can't find map index", indexfile);
605 	free(indexfile);
606 	for(i=0,p=patch[0];i<18*36;i++,p++)
607 		*p = 1;
608 	while(!feof(xfile) && fscanf(xfile,"%d%d%ld",&i,&j,&b)==3)
609 		patch[i+9][j+18] = b;
610 	fclose(xfile);
611 	ifile = fopen(mapfile,"r");
612 	if(ifile==NULL)
613 		filerror("can't find map data", mapfile);
614 	for(lat=lolat;lat<hilat;lat+=10.)
615 		for(lon=lolon;lon<hilon;lon+=10.) {
616 			if(!seeable(lat,lon))
617 				continue;
618 			i = pnorm(lat);
619 			j = pnorm(lon);
620 			if((b=patch[i+9][j+18])&1)
621 				continue;
622 			fseek(ifile,b,0);
623 			while((ip=getc(ifile))>=0&&(jp=getc(ifile))>=0){
624 				if(ip!=(i&0377)||jp!=(j&0377))
625 					break;
626 				n = getshort(ifile);
627 				conn = 0;
628 				if(n > 0) {	/* absolute coordinates */
629 					kx = ky = 0;	/* set */
630 					for(k=0;k<n;k++){
631 						kx = SCALERATIO*getshort(ifile);
632 						ky = SCALERATIO*getshort(ifile);
633 						if (((k%delta) != 0) && (k != (n-1)))
634 							continue;
635 						conv(kx,&g.nlat);
636 						conv(ky,&g.wlon);
637 						conn = plotpt(&g,conn);
638 					}
639 				} else {	/* differential, scaled by SCALERATI0 */
640 					n = -n;
641 					kx = SCALERATIO*getshort(ifile);
642 					ky = SCALERATIO*getshort(ifile);
643 					for(k=0; k<n; k++) {
644 						c = getc(ifile);
645 						if(c&0200) c|= ~0177;
646 						kx += c;
647 						c = getc(ifile);
648 						if(c&0200) c|= ~0177;
649 						ky += c;
650 						if(k%delta!=0&&k!=n-1)
651 							continue;
652 						conv(kx,&g.nlat);
653 						conv(ky,&g.wlon);
654 						conn = plotpt(&g,conn);
655 					}
656 				}
657 				if(k==1) {
658 					conv(kx,&g.nlat);
659 					conv(ky,&g.wlon);
660 					plotpt(&g,conn);
661 				}
662 			}
663 		}
664 	fclose(ifile);
665 }
666 
667 int
668 seeable(double lat0, double lon0)
669 {
670 	double x, y;
671 	double lat, lon;
672 	for(lat=lat0;lat<=lat0+10;lat+=2*grid[2])
673 		for(lon=lon0;lon<=lon0+10;lon+=2*grid[2])
674 			if(normproj(lat,lon,&x,&y)*vflag>0)
675 				return(1);
676 	return(0);
677 }
678 
679 void
680 satellite(struct file *t)
681 {
682 	char sym[50];
683 	char lbl[50];
684 	double scale;
685 	int conn;
686 	double lat,lon;
687 	struct place place;
688 	static FILE *ifile;
689 
690 	if(ifile == nil)
691 		ifile = stdin;
692 	if(t->name[0]!='-'||t->name[1]!=0) {
693 		fclose(ifile);
694 		if((ifile=fopen(t->name,"r"))==NULL)
695 			filerror("can't find track", t->name);
696 	}
697 	comment("track",t->name);
698 	colorx(t->color);
699 	pen(t->style);
700 	for(;;) {
701 		conn = 0;
702 		while(!feof(ifile) && fscanf(ifile,"%lf%lf",&lat,&lon)==2){
703 			latlon(lat,lon,&place);
704 			if(fscanf(ifile,"%1s",lbl) == 1) {
705 				if(strchr("+-.0123456789",*lbl)==0)
706 					break;
707 				ungetc(*lbl,ifile);
708 			}
709 			conn = plotpt(&place,conn);
710 		}
711 		if(feof(ifile))
712 			return;
713 		fscanf(ifile,"%[^\n]",lbl+1);
714 		switch(*lbl) {
715 		case '"':
716 			if(plotpt(&place,conn))
717 				text(lbl+1);
718 			break;
719 		case ':':
720 		case '!':
721 			if(sscanf(lbl+1,"%s %lf",sym,&scale) <= 1)
722 				scale = 1;
723 			if(plotpt(&place,conn?conn:-1)) {
724 				int r = *lbl=='!'?0:rflag?-1:1;
725 				pen(SOLID);
726 				if(putsym(&place,sym,scale,r) == 0)
727 					text(lbl);
728 				pen(t->style);
729 			}
730 			break;
731 		default:
732 			if(plotpt(&place,conn))
733 				text(lbl);
734 			break;
735 		}
736 	}
737 }
738 
739 int
740 pnorm(double x)
741 {
742 	int i;
743 	i = x/10.;
744 	i %= 36;
745 	if(i>=18) return(i-36);
746 	if(i<-18) return(i+36);
747 	return(i);
748 }
749 
750 void
751 error(char *s)
752 {
753 	fprintf(stderr,"map: \r\n%s\n",s);
754 	exits("error");
755 }
756 
757 void
758 filerror(char *s, char *f)
759 {
760 	fprintf(stderr,"\r\n%s %s\n",s,f);
761 	exits("error");
762 }
763 
764 char *
765 mapindex(char *s)
766 {
767 	char *t = malloc(strlen(s)+3);
768 	strcpy(t,s);
769 	strcat(t,".x");
770 	return t;
771 }
772 
773 #define NOPT 32767
774 static int ox = NOPT;
775 static int oy = NOPT;
776 
777 int
778 cpoint(int xi, int yi, int conn)
779 {
780 	int dx = abs(ox-xi);
781 	int dy = abs(oy-yi);
782 	if(!xflag && (xi<left||xi>=right || yi<bottom||yi>=top)) {
783 		ox = oy = NOPT;
784 		return 0;
785 	}
786 	if(conn == -1)		/* isolated plotting symbol */
787 		{}
788 	else if(!conn)
789 		point(xi,yi);
790 	else {
791 		if(dx+dy>longlines) {
792 			ox = oy = NOPT;	/* don't leap across cuts */
793 			return 0;
794 		}
795 		if(dx || dy)
796 			vec(xi,yi);
797 	}
798 	ox = xi, oy = yi;
799 	return dx+dy<=2? 2: 1;	/* 2=very near; see dogrid */
800 }
801 
802 
803 struct place oldg;
804 
805 int
806 plotpt(struct place *g, int conn)
807 {
808 	int kx,ky;
809 	int ret;
810 	double cutlon;
811 	if(!inlimits(g)) {
812 		return(0);
813 }
814 	normalize(g);
815 	if(!inwindow(g)) {
816 		return(0);
817 }
818 	switch((*cut)(g,&oldg,&cutlon)) {
819 	case 2:
820 		if(conn) {
821 			ret = duple(g,cutlon)|duple(g,cutlon);
822 			oldg = *g;
823 			return(ret);
824 		}
825 	case 0:
826 		conn = 0;
827 	default:	/* prevent diags about bad return value */
828 	case 1:
829 		oldg = *g;
830 		ret = doproj(g,&kx,&ky);
831 		if(ret==0 || !onlimb && ret*vflag<=0)
832 			return(0);
833 		ret = cpoint(kx,ky,conn);
834 		return ret;
835 	}
836 }
837 
838 int
839 doproj(struct place *g, int *kx, int *ky)
840 {
841 	int i;
842 	double x,y,x1,y1;
843 /*fprintf(stderr,"dopr1 %f %f \n",g->nlat.l,g->wlon.l);*/
844 	i = fixproj(g,&x,&y);
845 	if(i == 0)
846 		return(0);
847 	if(rflag)
848 		x = -x;
849 /*fprintf(stderr,"dopr2 %f %f\n",x,y);*/
850 	if(!inpoly(x,y)) {
851 		return 0;
852 }
853 	x1 = x - xcent;
854 	y1 = y - ycent;
855 	x = (x1*crot.c - y1*crot.s + xoff)*scaling;
856 	y = (x1*crot.s + y1*crot.c + yoff)*scaling;
857 	*kx = x + (x>0?.5:-.5);
858 	*ky = y + (y>0?.5:-.5);
859 	return(i);
860 }
861 
862 int
863 duple(struct place *g, double cutlon)
864 {
865 	int kx,ky;
866 	int okx,oky;
867 	struct place ig;
868 	revlon(g,cutlon);
869 	revlon(&oldg,cutlon);
870 	ig = *g;
871 	invert(&ig);
872 	if(!inlimits(&ig))
873 		return(0);
874 	if(doproj(g,&kx,&ky)*vflag<=0 ||
875 	   doproj(&oldg,&okx,&oky)*vflag<=0)
876 		return(0);
877 	cpoint(okx,oky,0);
878 	cpoint(kx,ky,1);
879 	return(1);
880 }
881 
882 void
883 revlon(struct place *g, double cutlon)
884 {
885 	g->wlon.l = reduce(cutlon-reduce(g->wlon.l-cutlon));
886 	sincos(&g->wlon);
887 }
888 
889 
890 /*	recognize problems of cuts
891  *	move a point across cut to side of its predecessor
892  *	if its very close to the cut
893  *	return(0) if cut interrupts the line
894  *	return(1) if line is to be drawn normally
895  *	return(2) if line is so close to cut as to
896  *	be properly drawn on both sheets
897 */
898 
899 int
900 picut(struct place *g, struct place *og, double *cutlon)
901 {
902 	*cutlon = PI;
903 	return(ckcut(g,og,PI));
904 }
905 
906 int
907 nocut(struct place *g, struct place *og, double *cutlon)
908 {
909 	USED(g);
910 	USED(og);
911 	USED(cutlon);
912 /*
913 #pragma	ref g
914 #pragma	ref og
915 #pragma	ref cutlon
916 */
917 	return(1);
918 }
919 
920 int
921 ckcut(struct place *g1, struct place *g2, double lon)
922 {
923 	double d1, d2;
924 	double f1, f2;
925 	int kx,ky;
926 	d1 = reduce(g1->wlon.l -lon);
927 	d2 = reduce(g2->wlon.l -lon);
928 	if((f1=fabs(d1))<FUZZ)
929 		d1 = diddle(g1,lon,d2);
930 	if((f2=fabs(d2))<FUZZ) {
931 		d2 = diddle(g2,lon,d1);
932 		if(doproj(g2,&kx,&ky)*vflag>0)
933 			cpoint(kx,ky,0);
934 	}
935 	if(f1<FUZZ&&f2<FUZZ)
936 		return(2);
937 	if(f1>PI*TWO_THRD||f2>PI*TWO_THRD)
938 		return(1);
939 	return(d1*d2>=0);
940 }
941 
942 double
943 diddle(struct place *g, double lon, double d)
944 {
945 	double d1;
946 	d1 = FUZZ/2;
947 	if(d<0)
948 		d1 = -d1;
949 	g->wlon.l = reduce(lon+d1);
950 	sincos(&g->wlon);
951 	return(d1);
952 }
953 
954 double
955 reduce(double lon)
956 {
957 	if(lon>PI)
958 		lon -= 2*PI;
959 	else if(lon<-PI)
960 		lon += 2*PI;
961 	return(lon);
962 }
963 
964 
965 double tetrapt = 35.26438968;	/* atan(1/sqrt(2)) */
966 
967 void
968 dogrid(double lat0, double lat1, double lon0, double lon1)
969 {
970 	double slat,slon,tlat,tlon;
971 	register int conn, oconn;
972 	slat = tlat = slon = tlon = 0;
973 	if(lat1>lat0)
974 		slat = tlat = fmin(grid[2],dlat);
975 	else
976 		slon = tlon = fmin(grid[2],dlon);;
977 	conn = oconn = 0;
978 	while(lat0<=lat1&&lon0<=lon1) {
979 		conn = gridpt(lat0,lon0,conn);
980 		if(projection==Xguyou&&slat>0) {
981 			if(lat0<-45&&lat0+slat>-45)
982 				conn = gridpt(-45.,lon0,conn);
983 			else if(lat0<45&&lat0+slat>45)
984 				conn = gridpt(45.,lon0,conn);
985 		} else if(projection==Xtetra&&slat>0) {
986 			if(lat0<-tetrapt&&lat0+slat>-tetrapt) {
987 				gridpt(-tetrapt-.001,lon0,conn);
988 				conn = gridpt(-tetrapt+.001,lon0,0);
989 			}
990 			else if(lat0<tetrapt&&lat0+slat>tetrapt) {
991 				gridpt(tetrapt-.001,lon0,conn);
992 				conn = gridpt(tetrapt+.001,lon0,0);
993 			}
994 		}
995 		if(conn==0 && oconn!=0) {
996 			if(slat+slon>.05) {
997 				lat0 -= slat;	/* steps too big */
998 				lon0 -= slon;	/* or near bdry */
999 				slat /= 2;
1000 				slon /= 2;
1001 				conn = oconn = gridpt(lat0,lon0,conn);
1002 			} else
1003 				oconn = 0;
1004 		} else {
1005 			if(conn==2) {
1006 				slat = tlat;
1007 				slon = tlon;
1008 				conn = 1;
1009 			}
1010 			oconn = conn;
1011 	 	}
1012 		lat0 += slat;
1013 		lon0 += slon;
1014 	}
1015 	gridpt(lat1,lon1,conn);
1016 }
1017 
1018 static int gridinv;		/* nonzero when doing window bounds */
1019 
1020 int
1021 gridpt(double lat, double lon, int conn)
1022 {
1023 	struct place g;
1024 /*fprintf(stderr,"%f %f\n",lat,lon);*/
1025 	latlon(lat,lon,&g);
1026 	if(gridinv)
1027 		invert(&g);
1028 	return(plotpt(&g,conn));
1029 }
1030 
1031 /* win=0 ordinary grid lines, win=1 window lines */
1032 
1033 void
1034 dobounds(double lolat, double hilat, double lolon, double hilon, int win)
1035 {
1036 	gridinv = win;
1037 	if(lolat>-90 || win && (poles&1)!=0)
1038 		dogrid(lolat+FUZZ,lolat+FUZZ,lolon,hilon);
1039 	if(hilat<90 || win && (poles&2)!=0)
1040 		dogrid(hilat-FUZZ,hilat-FUZZ,lolon,hilon);
1041 	if(hilon-lolon<360 || win && cut==picut) {
1042 		dogrid(lolat,hilat,lolon+FUZZ,lolon+FUZZ);
1043 		dogrid(lolat,hilat,hilon-FUZZ,hilon-FUZZ);
1044 	}
1045 	gridinv = 0;
1046 }
1047 
1048 static void
1049 dolimb(void)
1050 {
1051 	double lat, lon;
1052 	double res = fmin(dlat, dlon)/4;
1053 	int conn = 0;
1054 	int newconn;
1055 	if(limb == 0)
1056 		return;
1057 	onlimb = gridinv = 1;
1058 	for(;;) {
1059 		newconn = (*limb)(&lat, &lon, res);
1060 		if(newconn == -1)
1061 			break;
1062 		conn = gridpt(lat, lon, conn*newconn);
1063 	}
1064 	onlimb = gridinv = 0;
1065 }
1066 
1067 
1068 void
1069 radbds(double *w, double *rw)
1070 {
1071 	int i;
1072 	for(i=0;i<4;i++)
1073 		rw[i] = w[i]*RAD;
1074 	rw[0] -= FUZZ;
1075 	rw[1] += FUZZ;
1076 	rw[2] -= FUZZ;
1077 	rw[3] += FUZZ;
1078 }
1079 
1080 void
1081 windlim(void)
1082 {
1083 	double center = orientation[0];
1084 	double colat;
1085 	if(center>90)
1086 		center = 180 - center;
1087 	if(center<-90)
1088 		center = -180 - center;
1089 	if(fabs(center)>90)
1090 		error("unreasonable orientation");
1091 	colat = 90 - window[0];
1092 	if(center-colat>limits[0])
1093 		limits[0] = center - colat;
1094 	if(center+colat<limits[1])
1095 		limits[1] = center + colat;
1096 }
1097 
1098 
1099 short
1100 getshort(FILE *f)
1101 {
1102 	int c, r;
1103 	c = getc(f);
1104 	r = (c | getc(f)<<8);
1105 	if (r&0x8000)
1106 		r |= ~0xFFFF;	/* in case short > 16 bits */
1107 	return r;
1108 }
1109 
1110 double
1111 fmin(double x, double y)
1112 {
1113 	return(x<y?x:y);
1114 }
1115 
1116 double
1117 fmax(double x, double y)
1118 {
1119 	return(x>y?x:y);
1120 }
1121 
1122 void
1123 clamp(double *px, double v)
1124 {
1125 	*px = (v<0?fmax:fmin)(*px,v);
1126 }
1127 
1128 void
1129 pathnames(void)
1130 {
1131 	int i;
1132 	char *t, *indexfile, *name;
1133 	FILE *f, *fx;
1134 	for(i=0; i<nfile; i++) {
1135 		name = file[i].name;
1136 		if(*name=='/')
1137 			continue;
1138 		indexfile = mapindex(name);
1139 			/* ansi equiv of unix access() call */
1140 		f = fopen(name, "r");
1141 		fx = fopen(indexfile, "r");
1142 		if(f) fclose(f);
1143 		if(fx) fclose(fx);
1144 		free(indexfile);
1145 		if(f && fx)
1146 			continue;
1147 		t = malloc(strlen(name)+strlen(mapdir)+2);
1148 		strcpy(t,mapdir);
1149 		strcat(t,"/");
1150 		strcat(t,name);
1151 		file[i].name = t;
1152 	}
1153 }
1154 
1155 void
1156 clipinit(void)
1157 {
1158 	int i;
1159 	double s,t;
1160 	if(nvert<=0)
1161 		return;
1162 	for(i=0; i<nvert; i++) {	/*convert latlon to xy*/
1163 		if(normproj(v[i].x,v[i].y,&v[i].x,&v[i].y)==0)
1164 			error("invisible clipping vertex");
1165 	}
1166 	if(nvert==2) {			/*rectangle with diag specified*/
1167 		nvert = 4;
1168 		v[2] = v[1];
1169 		v[1].x=v[0].x, v[1].y=v[2].y, v[3].x=v[2].x, v[3].y=v[0].y;
1170 	}
1171 	v[nvert] = v[0];
1172 	v[nvert+1] = v[1];
1173 	s = 0;
1174 	for(i=1; i<=nvert; i++) {	/*test for convexity*/
1175 		t = (v[i-1].x-v[i].x)*(v[i+1].y-v[i].y) -
1176 		    (v[i-1].y-v[i].y)*(v[i+1].x-v[i].x);
1177 		if(t<-FUZZ && s>=0) s = 1;
1178 		if(t>FUZZ && s<=0) s = -1;
1179 		if(-FUZZ<=t&&t<=FUZZ || t*s>0) {
1180 			s = 0;
1181 			break;
1182 		}
1183 	}
1184 	if(s==0)
1185 		error("improper clipping polygon");
1186 	for(i=0; i<nvert; i++) {	/*edge equation ax+by=c*/
1187 		e[i].a = s*(v[i+1].y - v[i].y);
1188 		e[i].b = s*(v[i].x - v[i+1].x);
1189 		e[i].c = s*(v[i].x*v[i+1].y - v[i].y*v[i+1].x);
1190 	}
1191 }
1192 
1193 int
1194 inpoly(double x, double y)
1195 {
1196 	int i;
1197 	for(i=0; i<nvert; i++) {
1198 		register struct edge *ei = &e[i];
1199 		double val = x*ei->a + y*ei->b - ei->c;
1200 		if(val>10*FUZZ)
1201 			return(0);
1202 	}
1203 	return 1;
1204 }
1205 
1206 void
1207 realcut(void)
1208 {
1209 	struct place g;
1210 	double lat;
1211 
1212 	if(cut != picut)	/* punt on unusual cuts */
1213 		return;
1214 	for(lat=window[0]; lat<=window[1]; lat+=grid[2]) {
1215 		g.wlon.l = PI;
1216 		sincos(&g.wlon);
1217 		g.nlat.l = lat*RAD;
1218 		sincos(&g.nlat);
1219 		if(!inwindow(&g)) {
1220 			break;
1221 }
1222 		invert(&g);
1223 		if(inlimits(&g)) {
1224 			return;
1225 }
1226 	}
1227 	longlines = shortlines = LONGLINES;
1228 	cut = nocut;		/* not necessary; small eff. gain */
1229 }