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