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
option(char * s)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
conv(int k,struct coord * g)134 conv(int k, struct coord *g)
135 {
136 g->l = (0.0001/SCALERATIO)*k;
137 sincos(g);
138 }
139
140 int
main(int argc,char * argv[])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
fixproj(struct place * g,double * x,double * y)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
normproj(double lat,double lon,double * x,double * y)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
posproj(double lat,double lon,double * x,double * y)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
inwindow(struct place * geog)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
inlimits(struct place * g)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
getdata(char * mapfile)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
seeable(double lat0,double lon0)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
satellite(struct file * t)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
pnorm(double x)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
error(char * s)745 error(char *s)
746 {
747 fprintf(stderr,"map: \r\n%s\n",s);
748 exits("error");
749 }
750
751 void
filerror(char * s,char * f)752 filerror(char *s, char *f)
753 {
754 fprintf(stderr,"\r\n%s %s\n",s,f);
755 exits("error");
756 }
757
758 char *
mapindex(char * s)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
cpoint(int xi,int yi,int conn)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
plotpt(struct place * g,int conn)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
doproj(struct place * g,int * kx,int * ky)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
duple(struct place * g,double cutlon)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
revlon(struct place * g,double cutlon)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
picut(struct place * g,struct place * og,double * cutlon)893 picut(struct place *g, struct place *og, double *cutlon)
894 {
895 *cutlon = PI;
896 return(ckcut(g,og,PI));
897 }
898
899 int
nocut(struct place * g,struct place * og,double * cutlon)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
ckcut(struct place * g1,struct place * g2,double lon)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
diddle(struct place * g,double lon,double d)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
reduce(double lon)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
dogrid(double lat0,double lat1,double lon0,double lon1)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
gridpt(double lat,double lon,int conn)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
dobounds(double lolat,double hilat,double lolon,double hilon,int win)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
dolimb(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
radbds(double * w,double * rw)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
windlim(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
getshort(FILE * f)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
fmin(double x,double y)1102 fmin(double x, double y)
1103 {
1104 return(x<y?x:y);
1105 }
1106
1107 double
fmax(double x,double y)1108 fmax(double x, double y)
1109 {
1110 return(x>y?x:y);
1111 }
1112
1113 void
clamp(double * px,double v)1114 clamp(double *px, double v)
1115 {
1116 *px = (v<0?fmax:fmin)(*px,v);
1117 }
1118
1119 void
pathnames(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
clipinit(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
inpoly(double x,double y)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
realcut()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 }
1221