Subversion Repositories planix.SVN

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 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 oriented;	/* nonzero if -o option occurred */
49 
static 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 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 
	register conn;
683 
	double lat,lon;
684 
	struct place place;
685 
	static FILE *ifile = stdin;
686 
	if(t->name[0]!='-'||t->name[1]!=0) {
687 
		fclose(ifile);
688 
		if((ifile=fopen(t->name,"r"))==NULL)
689 
			filerror("can't find track", t->name);
690 
	}
691 
	comment("track",t->name);
692 
	colorx(t->color);
693 
	pen(t->style);
694 
	for(;;) {
695 
		conn = 0;
696 
		while(!feof(ifile) && fscanf(ifile,"%lf%lf",&lat,&lon)==2){
697 
			latlon(lat,lon,&place);
698 
			if(fscanf(ifile,"%1s",lbl) == 1) {
699 
				if(strchr("+-.0123456789",*lbl)==0)
700 
					break;
701 
				ungetc(*lbl,ifile);
702 
			}
703 
			conn = plotpt(&place,conn);
704 
		}
705 
		if(feof(ifile))
706 
			return;
707 
		fscanf(ifile,"%[^\n]",lbl+1);
708 
		switch(*lbl) {
709 
		case '"':
710 
			if(plotpt(&place,conn))
711 
				text(lbl+1);
712 
			break;
713 
		case ':':
714 
		case '!':
715 
			if(sscanf(lbl+1,"%s %lf",sym,&scale) <= 1)
716 
				scale = 1;
717 
			if(plotpt(&place,conn?conn:-1)) {
718 
				int r = *lbl=='!'?0:rflag?-1:1;
719 
				pen(SOLID);
720 
				if(putsym(&place,sym,scale,r) == 0)
721 
					text(lbl);
722 
				pen(t->style);
723 
			}
724 
			break;
725 
		default:
726 
			if(plotpt(&place,conn))
727 
				text(lbl);
728 
			break;
729 
		}
730 
	}
731 
}
732 
 
733 
int
734 
pnorm(double x)
735 
{
736 
	int i;
737 
	i = x/10.;
738 
	i %= 36;
739 
	if(i>=18) return(i-36);
740 
	if(i<-18) return(i+36);
741 
	return(i);
742 
}
743 
 
744 
void
745 
error(char *s)
746 
{
747 
	fprintf(stderr,"map: \r\n%s\n",s);
748 
	exits("error");
749 
}
750 
 
751 
void
752 
filerror(char *s, char *f)
753 
{
754 
	fprintf(stderr,"\r\n%s %s\n",s,f);
755 
	exits("error");
756 
}
757 
 
758 
char *
759 
mapindex(char *s)
760 
{
761 
	char *t = malloc(strlen(s)+3);
762 
	strcpy(t,s);
763 
	strcat(t,".x");
764 
	return t;
765 
}
766 
 
767 
#define NOPT 32767
768 
static ox = NOPT, oy = NOPT;
769 
 
770 
int
771 
cpoint(int xi, int yi, int conn)
772 
{
773 
	int dx = abs(ox-xi);
774 
	int dy = abs(oy-yi);
775 
	if(!xflag && (xi<left||xi>=right || yi<bottom||yi>=top)) {
776 
		ox = oy = NOPT;
777 
		return 0;
778 
	}
779 
	if(conn == -1)		/* isolated plotting symbol */
780 
		{}
781 
	else if(!conn)
782 
		point(xi,yi);
783 
	else {
784 
		if(dx+dy>longlines) {
785 
			ox = oy = NOPT;	/* don't leap across cuts */
786 
			return 0;
787 
		}
788 
		if(dx || dy)
789 
			vec(xi,yi);
790 
	}
791 
	ox = xi, oy = yi;
792 
	return dx+dy<=2? 2: 1;	/* 2=very near; see dogrid */
793 
}
794 
 
795 
 
796 
struct place oldg;
797 
 
798 
int
799 
plotpt(struct place *g, int conn)
800 
{
801 
	int kx,ky;
802 
	int ret;
803 
	double cutlon;
804 
	if(!inlimits(g)) {
805 
		return(0);
806 
}
807 
	normalize(g);
808 
	if(!inwindow(g)) {
809 
		return(0);
810 
}
811 
	switch((*cut)(g,&oldg,&cutlon)) {
812 
	case 2:
813 
		if(conn) {
814 
			ret = duple(g,cutlon)|duple(g,cutlon);
815 
			oldg = *g;
816 
			return(ret);
817 
		}
818 
	case 0:
819 
		conn = 0;
820 
	default:	/* prevent diags about bad return value */
821 
	case 1:
822 
		oldg = *g;
823 
		ret = doproj(g,&kx,&ky);
824 
		if(ret==0 || !onlimb && ret*vflag<=0)
825 
			return(0);
826 
		ret = cpoint(kx,ky,conn);
827 
		return ret;
828 
	}
829 
}
830 
 
831 
int
832 
doproj(struct place *g, int *kx, int *ky)
833 
{
834 
	int i;
835 
	double x,y,x1,y1;
836 
/*fprintf(stderr,"dopr1 %f %f \n",g->nlat.l,g->wlon.l);*/
837 
	i = fixproj(g,&x,&y);
838 
	if(i == 0)
839 
		return(0);
840 
	if(rflag)
841 
		x = -x;
842 
/*fprintf(stderr,"dopr2 %f %f\n",x,y);*/
843 
	if(!inpoly(x,y)) {
844 
		return 0;
845 
}
846 
	x1 = x - xcent;
847 
	y1 = y - ycent;
848 
	x = (x1*crot.c - y1*crot.s + xoff)*scaling;
849 
	y = (x1*crot.s + y1*crot.c + yoff)*scaling;
850 
	*kx = x + (x>0?.5:-.5);
851 
	*ky = y + (y>0?.5:-.5);
852 
	return(i);
853 
}
854 
 
855 
int
856 
duple(struct place *g, double cutlon)
857 
{
858 
	int kx,ky;
859 
	int okx,oky;
860 
	struct place ig;
861 
	revlon(g,cutlon);
862 
	revlon(&oldg,cutlon);
863 
	ig = *g;
864 
	invert(&ig);
865 
	if(!inlimits(&ig))
866 
		return(0);
867 
	if(doproj(g,&kx,&ky)*vflag<=0 ||
868 
	   doproj(&oldg,&okx,&oky)*vflag<=0)
869 
		return(0);
870 
	cpoint(okx,oky,0);
871 
	cpoint(kx,ky,1);
872 
	return(1);
873 
}
874 
 
875 
void
876 
revlon(struct place *g, double cutlon)
877 
{
878 
	g->wlon.l = reduce(cutlon-reduce(g->wlon.l-cutlon));
879 
	sincos(&g->wlon);
880 
}
881 
 
882 
 
883 
/*	recognize problems of cuts
884 
 *	move a point across cut to side of its predecessor
885 
 *	if its very close to the cut
886 
 *	return(0) if cut interrupts the line
887 
 *	return(1) if line is to be drawn normally
888 
 *	return(2) if line is so close to cut as to
889 
 *	be properly drawn on both sheets
890 
*/
891 
 
892 
int
893 
picut(struct place *g, struct place *og, double *cutlon)
894 
{
895 
	*cutlon = PI;
896 
	return(ckcut(g,og,PI));
897 
}
898 
 
899 
int
900 
nocut(struct place *g, struct place *og, double *cutlon)
901 
{
902 
	USED(g, og, cutlon);
903 
/*
904 
#pragma	ref g
905 
#pragma	ref og
906 
#pragma	ref cutlon
907 
*/
908 
	return(1);
909 
}
910 
 
911 
int
912 
ckcut(struct place *g1, struct place *g2, double lon)
913 
{
914 
	double d1, d2;
915 
	double f1, f2;
916 
	int kx,ky;
917 
	d1 = reduce(g1->wlon.l -lon);
918 
	d2 = reduce(g2->wlon.l -lon);
919 
	if((f1=fabs(d1))<FUZZ)
920 
		d1 = diddle(g1,lon,d2);
921 
	if((f2=fabs(d2))<FUZZ) {
922 
		d2 = diddle(g2,lon,d1);
923 
		if(doproj(g2,&kx,&ky)*vflag>0)
924 
			cpoint(kx,ky,0);
925 
	}
926 
	if(f1<FUZZ&&f2<FUZZ)
927 
		return(2);
928 
	if(f1>PI*TWO_THRD||f2>PI*TWO_THRD)
929 
		return(1);
930 
	return(d1*d2>=0);
931 
}
932 
 
933 
double
934 
diddle(struct place *g, double lon, double d)
935 
{
936 
	double d1;
937 
	d1 = FUZZ/2;
938 
	if(d<0)
939 
		d1 = -d1;
940 
	g->wlon.l = reduce(lon+d1);
941 
	sincos(&g->wlon);
942 
	return(d1);
943 
}
944 
 
945 
double
946 
reduce(double lon)
947 
{
948 
	if(lon>PI)
949 
		lon -= 2*PI;
950 
	else if(lon<-PI)
951 
		lon += 2*PI;
952 
	return(lon);
953 
}
954 
 
955 
 
956 
double tetrapt = 35.26438968;	/* atan(1/sqrt(2)) */
957 
 
958 
void
959 
dogrid(double lat0, double lat1, double lon0, double lon1)
960 
{
961 
	double slat,slon,tlat,tlon;
962 
	register int conn, oconn;
963 
	slat = tlat = slon = tlon = 0;
964 
	if(lat1>lat0)
965 
		slat = tlat = fmin(grid[2],dlat);
966 
	else
967 
		slon = tlon = fmin(grid[2],dlon);;
968 
	conn = oconn = 0;
969 
	while(lat0<=lat1&&lon0<=lon1) {
970 
		conn = gridpt(lat0,lon0,conn);
971 
		if(projection==Xguyou&&slat>0) {
972 
			if(lat0<-45&&lat0+slat>-45)
973 
				conn = gridpt(-45.,lon0,conn);
974 
			else if(lat0<45&&lat0+slat>45)
975 
				conn = gridpt(45.,lon0,conn);
976 
		} else if(projection==Xtetra&&slat>0) {
977 
			if(lat0<-tetrapt&&lat0+slat>-tetrapt) {
978 
				gridpt(-tetrapt-.001,lon0,conn);
979 
				conn = gridpt(-tetrapt+.001,lon0,0);
980 
			}
981 
			else if(lat0<tetrapt&&lat0+slat>tetrapt) {
982 
				gridpt(tetrapt-.001,lon0,conn);
983 
				conn = gridpt(tetrapt+.001,lon0,0);
984 
			}
985 
		}
986 
		if(conn==0 && oconn!=0) {
987 
			if(slat+slon>.05) {
988 
				lat0 -= slat;	/* steps too big */
989 
				lon0 -= slon;	/* or near bdry */
990 
				slat /= 2;
991 
				slon /= 2;
992 
				conn = oconn = gridpt(lat0,lon0,conn);
993 
			} else
994 
				oconn = 0;
995 
		} else {
996 
			if(conn==2) {
997 
				slat = tlat;
998 
				slon = tlon;
999 
				conn = 1;
1000 
			}
1001 
			oconn = conn;
1002 
	 	}
1003 
		lat0 += slat;
1004 
		lon0 += slon;
1005 
	}
1006 
	gridpt(lat1,lon1,conn);
1007 
}
1008 
 
1009 
static gridinv;		/* nonzero when doing window bounds */
1010 
 
1011 
int
1012 
gridpt(double lat, double lon, int conn)
1013 
{
1014 
	struct place g;
1015 
/*fprintf(stderr,"%f %f\n",lat,lon);*/
1016 
	latlon(lat,lon,&g);
1017 
	if(gridinv)
1018 
		invert(&g);
1019 
	return(plotpt(&g,conn));
1020 
}
1021 
 
1022 
/* win=0 ordinary grid lines, win=1 window lines */
1023 
 
1024 
void
1025 
dobounds(double lolat, double hilat, double lolon, double hilon, int win)
1026 
{
1027 
	gridinv = win;
1028 
	if(lolat>-90 || win && (poles&1)!=0)
1029 
		dogrid(lolat+FUZZ,lolat+FUZZ,lolon,hilon);
1030 
	if(hilat<90 || win && (poles&2)!=0)
1031 
		dogrid(hilat-FUZZ,hilat-FUZZ,lolon,hilon);
1032 
	if(hilon-lolon<360 || win && cut==picut) {
1033 
		dogrid(lolat,hilat,lolon+FUZZ,lolon+FUZZ);
1034 
		dogrid(lolat,hilat,hilon-FUZZ,hilon-FUZZ);
1035 
	}
1036 
	gridinv = 0;
1037 
}
1038 
 
1039 
static void
1040 
dolimb(void)
1041 
{
1042 
	double lat, lon;
1043 
	double res = fmin(dlat, dlon)/4;
1044 
	int conn = 0;
1045 
	int newconn;
1046 
	if(limb == 0)
1047 
		return;
1048 
	onlimb = gridinv = 1;
1049 
	for(;;) {
1050 
		newconn = (*limb)(&lat, &lon, res);
1051 
		if(newconn == -1)
1052 
			break;
1053 
		conn = gridpt(lat, lon, conn*newconn);
1054 
	}
1055 
	onlimb = gridinv = 0;
1056 
}
1057 
 
1058 
 
1059 
void
1060 
radbds(double *w, double *rw)
1061 
{
1062 
	int i;
1063 
	for(i=0;i<4;i++)
1064 
		rw[i] = w[i]*RAD;
1065 
	rw[0] -= FUZZ;
1066 
	rw[1] += FUZZ;
1067 
	rw[2] -= FUZZ;
1068 
	rw[3] += FUZZ;
1069 
}
1070 
 
1071 
void
1072 
windlim(void)
1073 
{
1074 
	double center = orientation[0];
1075 
	double colat;
1076 
	if(center>90)
1077 
		center = 180 - center;
1078 
	if(center<-90)
1079 
		center = -180 - center;
1080 
	if(fabs(center)>90)
1081 
		error("unreasonable orientation");
1082 
	colat = 90 - window[0];
1083 
	if(center-colat>limits[0])
1084 
		limits[0] = center - colat;
1085 
	if(center+colat<limits[1])
1086 
		limits[1] = center + colat;
1087 
}
1088 
 
1089 
 
1090 
short
1091 
getshort(FILE *f)
1092 
{
1093 
	int c, r;
1094 
	c = getc(f);
1095 
	r = (c | getc(f)<<8);
1096 
	if (r&0x8000)
1097 
		r |= ~0xFFFF;	/* in case short > 16 bits */
1098 
	return r;
1099 
}
1100 
 
1101 
double
1102 
fmin(double x, double y)
1103 
{
1104 
	return(x<y?x:y);
1105 
}
1106 
 
1107 
double
1108 
fmax(double x, double y)
1109 
{
1110 
	return(x>y?x:y);
1111 
}
1112 
 
1113 
void
1114 
clamp(double *px, double v)
1115 
{
1116 
	*px = (v<0?fmax:fmin)(*px,v);
1117 
}
1118 
 
1119 
void
1120 
pathnames(void)
1121 
{
1122 
	int i;
1123 
	char *t, *indexfile, *name;
1124 
	FILE *f, *fx;
1125 
	for(i=0; i<nfile; i++) {
1126 
		name = file[i].name;
1127 
		if(*name=='/')
1128 
			continue;
1129 
		indexfile = mapindex(name);
1130 
			/* ansi equiv of unix access() call */
1131 
		f = fopen(name, "r");
1132 
		fx = fopen(indexfile, "r");
1133 
		if(f) fclose(f);
1134 
		if(fx) fclose(fx);
1135 
		free(indexfile);
1136 
		if(f && fx)
1137 
			continue;
1138 
		t = malloc(strlen(name)+strlen(mapdir)+2);
1139 
		strcpy(t,mapdir);
1140 
		strcat(t,"/");
1141 
		strcat(t,name);
1142 
		file[i].name = t;
1143 
	}
1144 
}
1145 
 
1146 
void
1147 
clipinit(void)
1148 
{
1149 
	register i;
1150 
	double s,t;
1151 
	if(nvert<=0)
1152 
		return;
1153 
	for(i=0; i<nvert; i++) {	/*convert latlon to xy*/
1154 
		if(normproj(v[i].x,v[i].y,&v[i].x,&v[i].y)==0)
1155 
			error("invisible clipping vertex");
1156 
	}
1157 
	if(nvert==2) {			/*rectangle with diag specified*/
1158 
		nvert = 4;
1159 
		v[2] = v[1];
1160 
		v[1].x=v[0].x, v[1].y=v[2].y, v[3].x=v[2].x, v[3].y=v[0].y;
1161 
	}
1162 
	v[nvert] = v[0];
1163 
	v[nvert+1] = v[1];
1164 
	s = 0;
1165 
	for(i=1; i<=nvert; i++) {	/*test for convexity*/
1166 
		t = (v[i-1].x-v[i].x)*(v[i+1].y-v[i].y) -
1167 
		    (v[i-1].y-v[i].y)*(v[i+1].x-v[i].x);
1168 
		if(t<-FUZZ && s>=0) s = 1;
1169 
		if(t>FUZZ && s<=0) s = -1;
1170 
		if(-FUZZ<=t&&t<=FUZZ || t*s>0) {
1171 
			s = 0;
1172 
			break;
1173 
		}
1174 
	}
1175 
	if(s==0)
1176 
		error("improper clipping polygon");
1177 
	for(i=0; i<nvert; i++) {	/*edge equation ax+by=c*/
1178 
		e[i].a = s*(v[i+1].y - v[i].y);
1179 
		e[i].b = s*(v[i].x - v[i+1].x);
1180 
		e[i].c = s*(v[i].x*v[i+1].y - v[i].y*v[i+1].x);
1181 
	}
1182 
}
1183 
 
1184 
int
1185 
inpoly(double x, double y)
1186 
{
1187 
	register i;
1188 
	for(i=0; i<nvert; i++) {
1189 
		register struct edge *ei = &e[i];
1190 
		double val = x*ei->a + y*ei->b - ei->c;
1191 
		if(val>10*FUZZ)
1192 
			return(0);
1193 
	}
1194 
	return 1;
1195 
}
1196 
 
1197 
void
1198 
realcut()
1199 
{
1200 
	struct place g;
1201 
	double lat;
1202 
 
1203 
	if(cut != picut)	/* punt on unusual cuts */
1204 
		return;
1205 
	for(lat=window[0]; lat<=window[1]; lat+=grid[2]) {
1206 
		g.wlon.l = PI;
1207 
		sincos(&g.wlon);
1208 
		g.nlat.l = lat*RAD;
1209 
		sincos(&g.nlat);
1210 
		if(!inwindow(&g)) {
1211 
			break;
1212 
}
1213 
		invert(&g);
1214 
		if(inlimits(&g)) {
1215 
			return;
1216 
}
1217 
	}
1218 
	longlines = shortlines = LONGLINES;
1219 
	cut = nocut;		/* not necessary; small eff. gain */
1220 
}