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