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