![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to map.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / map / export|
Return to map.c CVS log | Up to [Research Unix] / researchv10no / cmd / map / export |
1.1 root 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 = "/usr/dict"; /* 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: static int xflag = 0; /* 1 if option -x occurred */
86: int vflag = 1; /* -1 if option -v occurred */
87: static double position[3]; /* option -p */
88: static double center[3] = {0., 0., 0.}; /* option -c */
89: static struct coord crot; /* option -c */
90: static double grid[3] = { 10., 10., RESOL }; /* option -g */
91: static double dlat, dlon; /* resolution for tracing grid in lat and lon */
92: static double scaling; /* to compute final integer output */
93: static struct file *track; /* options -t and -u */
94: static int ntrack; /* number of tracks present */
95: static char *symbolfile; /* option -y */
96:
97: void clamp(double *px, double v);
98: void clipinit(void);
99: double diddle(struct place *, double, double);
100: double diddle(struct place *, double, double);
101: void dobounds(double, double, double, double, int);
102: void dogrid(double, double, double, double);
103: int duple(struct place *, double);
104: double fmax(double, double);
105: double fmin(double, double);
106: void getdata(char *);
107: int gridpt(double, double, int);
108: int inpoly(double, double);
109: int inwindow(struct place *);
110: void pathnames(void);
111: int pnorm(double);
112: void radbds(double *w, double *rw);
113: void revlon(struct place *, double);
114: void satellite(struct file *);
115: int seeable(double, double);
116: void windlim(void);
117: void realcut(void);
118:
119: int
120: option(char *s)
121: {
122:
123: if(s[0]=='-' && (s[1]<'0'||s[1]>'9'))
124: return(s[1]!='.'&&s[1]!=0);
125: else
126: return(0);
127: }
128:
129: void
130: conv(int k, struct coord *g)
131: {
132: g->l = (0.0001/SCALERATIO)*k;
133: sincos(g);
134: }
135:
136: int
137: main(int argc, char *argv[])
138: {
139: int i,k;
140: char *s, *t, *style;
141: double x, y;
142: double lat, lon;
143: double *wlim;
144: double dd;
145: if(sizeof(short)!=2)
146: abort(); /* getshort() won't work */
147: s = getenv("MAP");
148: if(s)
149: file[0].name = s;
150: s = getenv("MAPDIR");
151: if(s)
152: mapdir = s;
153: if(argc<=1)
154: error("usage: map projection params options");
155: for(k=0;index[k].name;k++) {
156: s = index[k].name;
157: t = argv[1];
158: while(*s == *t){
159: if(*s==0) goto found;
160: s++;
161: t++;
162: }
163: }
164: fprintf(stderr,"projections:\n");
165: for(i=0;index[i].name;i++) {
166: fprintf(stderr,"%s",index[i].name);
167: for(k=0; k<index[i].npar; k++)
168: fprintf(stderr," p%d", k);
169: fprintf(stderr,"\n");
170: }
171: exit(1);
172: found:
173: argv += 2;
174: argc -= 2;
175: cut = index[k].cut;
176: limb = index[k].limb;
177: poles = index[k].poles;
178: for(i=0;i<index[k].npar;i++) {
179: if(i>=argc||option(argv[i])) {
180: fprintf(stderr,"%s needs %d params\n",index[k].name,index[k].npar);
181: exit(1);
182: }
183: params[i] = atof(argv[i]);
184: }
185: argv += i;
186: argc -= i;
187: while(argc>0&&option(argv[0])) {
188: argc--;
189: argv++;
190: switch(argv[-1][1]) {
191: case 'm':
192: if(file == &dfltfile) {
193: file = 0;
194: nfile = 0;
195: }
196: while(argc && !option(*argv)) {
197: file = realloc(file,(nfile+1)*sizeof(*file));
198: file[nfile].name = *argv;
199: file[nfile].color = currcolor;
200: file[nfile].style = SOLID;
201: nfile++;
202: argv++;
203: argc--;
204: }
205: break;
206: case 'b':
207: bflag = 0;
208: for(nvert=0;nvert<NVERT&&argc>=2;nvert++) {
209: if(option(*argv))
210: break;
211: v[nvert].x = atof(*argv++);
212: argc--;
213: if(option(*argv))
214: break;
215: v[nvert].y = atof(*argv++);
216: argc--;
217: }
218: if(nvert>=NVERT)
219: error("too many clipping vertices");
220: break;
221: case 'g':
222: gridcolor = currcolor;
223: for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
224: grid[i] = atof(argv[i]);
225: switch(i) {
226: case 0:
227: grid[0] = grid[1] = 0.;
228: break;
229: case 1:
230: grid[1] = grid[0];
231: }
232: argc -= i;
233: argv += i;
234: break;
235: case 't':
236: style = SOLID;
237: goto casetu;
238: case 'u':
239: style = DOTDASH;
240: casetu:
241: while(argc && !option(*argv)) {
242: track = realloc(track,(ntrack+1)*sizeof(*track));
243: track[ntrack].name = *argv;
244: track[ntrack].color = currcolor;
245: track[ntrack].style = style;
246: ntrack++;
247: argv++;
248: argc--;
249: }
250: break;
251: case 'r':
252: rflag++;
253: break;
254: case 's':
255: switch(argv[-1][2]) {
256: case '1':
257: s1flag++;
258: break;
259: case 0: /* compatibility */
260: case '2':
261: s2flag++;
262: }
263: break;
264: case 'o':
265: for(i=0;i<3&&i<argc&&!option(argv[i]);i++)
266: orientation[i] = atof(argv[i]);
267: oriented++;
268: argv += i;
269: argc -= i;
270: break;
271: case 'l':
272: bordcolor = currcolor;
273: for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
274: limits[i] = atof(argv[i]);
275: argv += i;
276: argc -= i;
277: break;
278: case 'k':
279: kflag++;
280: for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
281: klimits[i] = atof(argv[i]);
282: argv += i;
283: argc -= i;
284: break;
285: case 'd':
286: if(argc>0&&!option(argv[0])) {
287: delta = atoi(argv[0]);
288: argv++;
289: argc--;
290: }
291: break;
292: case 'w':
293: bordcolor = currcolor;
294: windowed++;
295: for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
296: window[i] = atof(argv[i]);
297: argv += i;
298: argc -= i;
299: break;
300: case 'c':
301: for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
302: center[i] = atof(argv[i]);
303: argc -= i;
304: argv += i;
305: break;
306: case 'p':
307: for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
308: position[i] = atof(argv[i]);
309: argc -= i;
310: argv += i;
311: if(i!=3||position[2]<=0)
312: error("incomplete positioning");
313: break;
314: case 'y':
315: if(argc>0&&!option(argv[0])) {
316: symbolfile = argv[0];
317: argc--;
318: argv++;
319: }
320: break;
321: case 'v':
322: if(index[k].limb == 0)
323: error("-v does not apply here");
324: vflag = -1;
325: break;
326: case 'x':
327: xflag = 1;
328: break;
329: case 'C':
330: if(argc && !option(*argv)) {
331: currcolor = *argv;
332: argc--;
333: argv++;
334: }
335: break;
336: }
337: }
338: if(argc>0)
339: error("error in arguments");
340: pathnames();
341: clamp(&limits[0],-90.);
342: clamp(&limits[1],90.);
343: clamp(&klimits[0],-90.);
344: clamp(&klimits[1],90.);
345: clamp(&window[0],-90.);
346: clamp(&window[1],90.);
347: radbds(limits,rlimits);
348: limcase = limits[2]<-180.?0:
349: limits[3]>180.?2:
350: 1;
351: if(
352: window[0]>=window[1]||
353: window[2]>=window[3]||
354: window[0]>90.||
355: window[1]<-90.||
356: window[2]>180.||
357: window[3]<-180.)
358: error("unreasonable window");
359: windlim();
360: radbds(window,rwindow);
361: upright = orientation[0]==90? 1: orientation[0]==-90? -1: 0;
362: if(index[k].spheroid && !upright)
363: error("can't tilt the spheroid");
364: if(limits[2]>limits[3])
365: limits[3] += 360;
366: if(!oriented)
367: orientation[2] = (limits[2]+limits[3])/2;
368: orient(orientation[0],orientation[1],orientation[2]);
369: projection = (*index[k].prog)(params[0],params[1]);
370: if(projection == 0)
371: error("unreasonable projection parameters");
372: clipinit();
373: grid[0] = fabs(grid[0]);
374: grid[1] = fabs(grid[1]);
375: if(!kflag)
376: for(i=0;i<4;i++)
377: klimits[i] = limits[i];
378: if(klimits[2]>klimits[3])
379: klimits[3] += 360;
380: lolat = limits[0];
381: hilat = limits[1];
382: lolon = limits[2];
383: hilon = limits[3];
384: if(lolon>=hilon||lolat>=hilat||lolat<-90.||hilat>90.)
385: error("unreasonable limits");
386: wlim = kflag? klimits: window;
387: dlat = fmin(hilat-lolat,wlim[1]-wlim[0])/16;
388: dlon = fmin(hilon-lolon,wlim[3]-wlim[2])/32;
389: dd = fmax(dlat,dlon);
390: while(grid[2]>fmin(dlat,dlon)/2)
391: grid[2] /= 2;
392: realcut();
393: if(nvert<=0) {
394: for(lat=klimits[0];lat<klimits[1]+dd-FUZZ;lat+=dd) {
395: if(lat>klimits[1])
396: lat = klimits[1];
397: for(lon=klimits[2];lon<klimits[3]+dd-FUZZ;lon+=dd) {
398: i = (kflag?posproj:normproj)
399: (lat,lon+(lon<klimits[3]?FUZZ:-FUZZ),
400: &x,&y);
401: if(i*vflag <= 0)
402: continue;
403: if(x<xmin) xmin = x;
404: if(x>xmax) xmax = x;
405: if(y<ymin) ymin = y;
406: if(y>ymax) ymax = y;
407: }
408: }
409: } else {
410: for(i=0; i<nvert; i++) {
411: x = v[i].x;
412: y = v[i].y;
413: if(x<xmin) xmin = x;
414: if(x>xmax) xmax = x;
415: if(y<ymin) ymin = y;
416: if(y>ymax) ymax = y;
417: }
418: }
419: xrange = xmax - xmin;
420: yrange = ymax - ymin;
421: if(xrange<=0||yrange<=0)
422: error("map seems to be empty");
423: scaling = 2; /*plotting area from -1 to 1*/
424: if(position[2]!=0) {
425: if(posproj(position[0]-.5,position[1],&xcent,&ycent)==0||
426: posproj(position[0]+.5,position[1],&x,&y)==0)
427: error("unreasonable position");
428: scaling /= (position[2]*hypot(x-xcent,y-ycent));
429: if(posproj(position[0],position[1],&xcent,&ycent)==0)
430: error("unreasonable position");
431: } else {
432: scaling /= (xrange>yrange?xrange:yrange);
433: xcent = (xmin+xmax)/2;
434: ycent = (ymin+ymax)/2;
435: }
436: xoff = center[0]/scaling;
437: yoff = center[1]/scaling;
438: crot.l = center[2]*RAD;
439: sincos(&crot);
440: scaling *= HALFWIDTH*0.9;
441: if(symbolfile)
442: getsyms(symbolfile);
443: if(!s2flag) {
444: openpl();
445: erase();
446: }
447: range(left,bottom,right,top);
448: color(gridcolor);
449: pen(DOTTED);
450: if(grid[0]>0.)
451: for(lat=ceil(lolat/grid[0])*grid[0];
452: lat<=hilat;lat+=grid[0])
453: dogrid(lat,lat,lolon,hilon);
454: if(grid[1]>0.)
455: for(lon=ceil(lolon/grid[1])*grid[1];
456: lon<=hilon;lon+=grid[1])
457: dogrid(lolat,hilat,lon,lon);
458: color(bordcolor);
459: pen(SOLID);
460: if(bflag) {
461: dolimb();
462: dobounds(lolat,hilat,lolon,hilon,0);
463: dobounds(window[0],window[1],window[2],window[3],1);
464: }
465: lolat = floor(limits[0]/10)*10;
466: hilat = ceil(limits[1]/10)*10;
467: lolon = floor(limits[2]/10)*10;
468: hilon = ceil(limits[3]/10)*10;
469: if(lolon>hilon)
470: hilon += 360.;
471: /*do tracks first so as not to lose the standard input*/
472: for(i=0;i<ntrack;i++) {
473: longlines = LONGLINES;
474: satellite(&track[i]);
475: longlines = shortlines;
476: }
477: for(i=0;i<nfile;i++) {
478: pen(file[i].style);
479: color(file[i].color);
480: getdata(file[i].name);
481: }
482: move(right,bottom);
483: if(!s1flag)
484: closepl();
485: return 0;
486: }
487:
488: /* Out of perverseness (really to recover from a dubious,
489: but documented, convention) the returns from projection
490: functions (-1 unplottable, 0 wrong sheet, 1 good) are
491: recoded into -1 wrong sheet, 0 unplottable, 1 good. */
492:
493: int
494: fixproj(struct place *g, double *x, double *y)
495: {
496: int i = (*projection)(g,x,y);
497: return i<0? 0: i==0? -1: 1;
498: }
499:
500: int
501: normproj(double lat, double lon, double *x, double *y)
502: {
503: int i;
504: struct place geog;
505: latlon(lat,lon,&geog);
506: /*
507: printp(&geog);
508: */
509: normalize(&geog);
510: if(!inwindow(&geog))
511: return(-1);
512: i = fixproj(&geog,x,y);
513: if(rflag)
514: *x = -*x;
515: /*
516: printp(&geog);
517: fprintf(stderr,"%d %.3f %.3f\n",i,*x,*y);
518: */
519: return(i);
520: }
521:
522: int
523: posproj(double lat, double lon, double *x, double *y)
524: {
525: int i;
526: struct place geog;
527: latlon(lat,lon,&geog);
528: normalize(&geog);
529: i = fixproj(&geog,x,y);
530: if(rflag)
531: *x = -*x;
532: return(i);
533: }
534:
535: int
536: inwindow(struct place *geog)
537: {
538: if(geog->nlat.l<rwindow[0]-FUZZ||
539: geog->nlat.l>rwindow[1]+FUZZ||
540: geog->wlon.l<rwindow[2]-FUZZ||
541: geog->wlon.l>rwindow[3]+FUZZ)
542: return(0);
543: else return(1);
544: }
545:
546: int
547: inlimits(struct place *g)
548: {
549: if(rlimits[0]-FUZZ>g->nlat.l||
550: rlimits[1]+FUZZ<g->nlat.l)
551: return(0);
552: switch(limcase) {
553: case 0:
554: if(rlimits[2]+TWOPI-FUZZ>g->wlon.l&&
555: rlimits[3]+FUZZ<g->wlon.l)
556: return(0);
557: break;
558: case 1:
559: if(rlimits[2]-FUZZ>g->wlon.l||
560: rlimits[3]+FUZZ<g->wlon.l)
561: return(0);
562: break;
563: case 2:
564: if(rlimits[2]>g->wlon.l&&
565: rlimits[3]-TWOPI+FUZZ<g->wlon.l)
566: return(0);
567: break;
568: }
569: return(1);
570: }
571:
572:
573: long patch[18][36];
574:
575: void
576: getdata(char *mapfile)
577: {
578: char *indexfile;
579: int kx,ky,c;
580: int k;
581: long b;
582: long *p;
583: int ip, jp;
584: int n;
585: struct place g;
586: int i, j;
587: double lat, lon;
588: int conn;
589: FILE *ifile, *xfile;
590:
591: indexfile = mapindex(mapfile);
592: xfile = fopen(indexfile,"r");
593: if(xfile==NULL)
594: filerror("can't find map index", indexfile);
595: free(indexfile);
596: for(i=0,p=patch[0];i<18*36;i++,p++)
597: *p = 1;
598: while(!feof(xfile) && fscanf(xfile,"%d%d%ld",&i,&j,&b)==3)
599: patch[i+9][j+18] = b;
600: fclose(xfile);
601: ifile = fopen(mapfile,"r");
602: if(ifile==NULL)
603: filerror("can't find map data", mapfile);
604: for(lat=lolat;lat<hilat;lat+=10.)
605: for(lon=lolon;lon<hilon;lon+=10.) {
606: if(!seeable(lat,lon))
607: continue;
608: i = pnorm(lat);
609: j = pnorm(lon);
610: if((b=patch[i+9][j+18])&1)
611: continue;
612: fseek(ifile,b,0);
613: while((ip=getc(ifile))>=0&&(jp=getc(ifile))>=0){
614: if(ip!=(i&0377)||jp!=(j&0377))
615: break;
616: n = getshort(ifile);
617: conn = 0;
618: if(n > 0) { /* absolute coordinates */
619: for(k=0;k<n;k++){
620: kx = SCALERATIO*getshort(ifile);
621: ky = SCALERATIO*getshort(ifile);
622: if (((k%delta) != 0) && (k != (n-1)))
623: continue;
624: conv(kx,&g.nlat);
625: conv(ky,&g.wlon);
626: conn = plotpt(&g,conn);
627: }
628: } else { /* differential, scaled by SCALERATI0 */
629: n = -n;
630: kx = SCALERATIO*getshort(ifile);
631: ky = SCALERATIO*getshort(ifile);
632: for(k=0; k<n; k++) {
633: c = getc(ifile);
634: if(c&0200) c|= ~0177;
635: kx += c;
636: c = getc(ifile);
637: if(c&0200) c|= ~0177;
638: ky += c;
639: if(k%delta!=0&&k!=n-1)
640: continue;
641: conv(kx,&g.nlat);
642: conv(ky,&g.wlon);
643: conn = plotpt(&g,conn);
644: }
645: }
646: if(k==1) {
647: conv(kx,&g.nlat);
648: conv(ky,&g.wlon);
649: conn = plotpt(&g,conn);
650: }
651: }
652: }
653: fclose(ifile);
654: }
655:
656: int
657: seeable(double lat0, double lon0)
658: {
659: double x, y;
660: double lat, lon;
661: for(lat=lat0;lat<=lat0+10;lat+=2*grid[2])
662: for(lon=lon0;lon<=lon0+10;lon+=2*grid[2])
663: if(normproj(lat,lon,&x,&y)*vflag>0)
664: return(1);
665: return(0);
666: }
667:
668: void
669: satellite(struct file *t)
670: {
671: char sym[50];
672: char lbl[50];
673: double scale;
674: register conn;
675: double lat,lon;
676: struct place place;
677: static FILE *ifile = stdin;
678: if(t->name[0]!='-'||t->name[1]!=0) {
679: fclose(ifile);
680: if((ifile=fopen(t->name,"r"))==NULL)
681: filerror("can't find track", t->name);
682: }
683: color(t->color);
684: pen(t->style);
685: for(;;) {
686: conn = 0;
687: while(!feof(ifile) && fscanf(ifile,"%lf%lf",&lat,&lon)==2){
688: latlon(lat,lon,&place);
689: if(fscanf(ifile,"%1s",lbl) == 1) {
690: if(strchr("+-.0123456789",*lbl)==0)
691: break;
692: ungetc(*lbl,ifile);
693: }
694: conn = plotpt(&place,conn);
695: }
696: if(feof(ifile))
697: return;
698: fscanf(ifile,"%[^\n]",lbl+1);
699: switch(*lbl) {
700: case '"':
701: if(plotpt(&place,conn))
702: text(lbl+1);
703: break;
704: case ':':
705: case '!':
706: if(sscanf(lbl+1,"%s %lf",sym,&scale) <= 1)
707: scale = 1;
708: if(plotpt(&place,conn?conn:-1)) {
709: int r = *lbl=='!'?0:rflag?-1:1;
710: if(putsym(&place,sym,scale,r) == 0)
711: text(lbl);
712: }
713: break;
714: default:
715: if(plotpt(&place,conn))
716: text(lbl);
717: break;
718: }
719: }
720: }
721:
722: int
723: pnorm(double x)
724: {
725: int i;
726: i = x/10.;
727: i %= 36;
728: if(i>=18) return(i-36);
729: if(i<-18) return(i+36);
730: return(i);
731: }
732:
733: void
734: error(char *s)
735: {
736: fprintf(stderr,"map: \r\n%s\n",s);
737: exit(1);
738: }
739:
740: void
741: filerror(char *s, char *f)
742: {
743: fprintf(stderr,"\r\n%s %s\n",s,f);
744: exit(1);
745: }
746:
747: char *
748: mapindex(char *s)
749: {
750: char *t = malloc(strlen(s)+3);
751: strcpy(t,s);
752: strcat(t,".x");
753: return t;
754: }
755:
756: #define NOPT 32767
757: static ox = NOPT, oy = NOPT;
758:
759: int
760: cpoint(int xi, int yi, int conn)
761: {
762: int dx = abs(ox-xi);
763: int dy = abs(oy-yi);
764: if(!xflag && (xi<left||xi>=right || yi<bottom||yi>=top)) {
765: ox = oy = NOPT;
766: return 0;
767: }
768: if(conn == -1) /* isolated plotting symbol */
769: ;
770: else if(!conn)
771: move(xi,yi);
772: else {
773: if(dx+dy>longlines) {
774: ox = oy = NOPT; /* don't leap across cuts */
775: return 0;
776: }
777: if(dx || dy)
778: vec(xi,yi);
779: }
780: ox = xi, oy = yi;
781: return dx+dy<=2? 2: 1; /* 2=very near; see dogrid */
782: }
783:
784:
785: struct place oldg;
786:
787: int
788: plotpt(struct place *g, int conn)
789: {
790: int kx,ky;
791: int ret;
792: double cutlon;
793: if(!inlimits(g)) {
794: return(0);
795: }
796: normalize(g);
797: if(!inwindow(g)) {
798: return(0);
799: }
800: switch((*cut)(g,&oldg,&cutlon)) {
801: case 2:
802: if(conn) {
803: ret = duple(g,cutlon)|duple(g,cutlon);
804: oldg = *g;
805: return(ret);
806: }
807: case 0:
808: conn = 0;
809: default: /* prevent diags about bad return value */
810: case 1:
811: oldg = *g;
812: ret = doproj(g,&kx,&ky);
813: if(ret==0 || !onlimb && ret*vflag<=0)
814: return(0);
815: ret = cpoint(kx,ky,conn);
816: return ret;
817: }
818: }
819:
820: int
821: doproj(struct place *g, int *kx, int *ky)
822: {
823: int i;
824: double x,y,x1,y1;
825: /*fprintf(stderr,"dopr1 %f %f \n",g->nlat.l,g->wlon.l);*/
826: i = fixproj(g,&x,&y);
827: if(i == 0)
828: return(0);
829: if(rflag)
830: x = -x;
831: /*fprintf(stderr,"dopr2 %f %f\n",x,y);*/
832: if(!inpoly(x,y)) {
833: return 0;
834: }
835: x1 = x - xcent;
836: y1 = y - ycent;
837: x = (x1*crot.c - y1*crot.s + xoff)*scaling;
838: y = (x1*crot.s + y1*crot.c + yoff)*scaling;
839: *kx = x + (x>0?.5:-.5);
840: *ky = y + (y>0?.5:-.5);
841: return(i);
842: }
843:
844: int
845: duple(struct place *g, double cutlon)
846: {
847: int kx,ky;
848: int okx,oky;
849: struct place ig;
850: revlon(g,cutlon);
851: revlon(&oldg,cutlon);
852: ig = *g;
853: invert(&ig);
854: if(!inlimits(&ig))
855: return(0);
856: if(doproj(g,&kx,&ky)*vflag<=0 ||
857: doproj(&oldg,&okx,&oky)*vflag<=0)
858: return(0);
859: cpoint(okx,oky,0);
860: cpoint(kx,ky,1);
861: return(1);
862: }
863:
864: void
865: revlon(struct place *g, double cutlon)
866: {
867: g->wlon.l = reduce(cutlon-reduce(g->wlon.l-cutlon));
868: sincos(&g->wlon);
869: }
870:
871:
872: /* recognize problems of cuts
873: * move a point across cut to side of its predecessor
874: * if its very close to the cut
875: * return(0) if cut interrupts the line
876: * return(1) if line is to be drawn normally
877: * return(2) if line is so close to cut as to
878: * be properly drawn on both sheets
879: */
880:
881: int
882: picut(struct place *g, struct place *og, double *cutlon)
883: {
884: *cutlon = PI;
885: return(ckcut(g,og,PI));
886: }
887:
888: int
889: nocut(struct place *g, struct place *og, double *cutlon)
890: {
891: #pragma ref g
892: #pragma ref og
893: #pragma ref cutlon
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: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.