sys/src/cmd/map/map.c

   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 }