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 } 1221