op public repos

Blob

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