route.c (2716B)
1 #include <u.h> 2 #include <libc.h> 3 #include "map.h" 4 5 /* Given two lat-lon pairs, find an orientation for the 6 -o option of "map" that will place those two points 7 on the equator of a standard projection, equally spaced 8 about the prime meridian. 9 10 -w and -l options are suggested also. 11 12 Option -t prints out a series of 13 coordinates that follows the (great circle) track 14 in the original coordinate system, 15 followed by ". 16 This data is just right for map -t. 17 18 Option -i inverts the map top-to-bottom. 19 */ 20 struct place pole; 21 struct coord twist; 22 int track; 23 int inv = -1; 24 25 extern void doroute(double, double, double, double, double); 26 27 void 28 dorot(double a, double b, double *x, double *y, void (*f)(struct place *)) 29 { 30 struct place g; 31 deg2rad(a,&g.nlat); 32 deg2rad(b,&g.wlon); 33 (*f)(&g); 34 *x = g.nlat.l/RAD; 35 *y = g.wlon.l/RAD; 36 } 37 38 void 39 rotate(double a, double b, double *x, double *y) 40 { 41 dorot(a,b,x,y,normalize); 42 } 43 44 void 45 rinvert(double a, double b, double *x, double *y) 46 { 47 dorot(a,b,x,y,invert); 48 } 49 50 main(int argc, char **argv) 51 { 52 #pragma ref argv 53 double an,aw,bn,bw; 54 ARGBEGIN { 55 case 't': 56 track = 1; 57 break; 58 59 case 'i': 60 inv = 1; 61 break; 62 63 default: 64 exits("route: bad option"); 65 } ARGEND; 66 if (argc<4) { 67 print("use route [-t] [-i] lat lon lat lon\n"); 68 exits("arg count"); 69 } 70 an = atof(argv[0]); 71 aw = atof(argv[1]); 72 bn = atof(argv[2]); 73 bw = atof(argv[3]); 74 doroute(inv*90.,an,aw,bn,bw); 75 return 0; 76 } 77 78 void 79 doroute(double dir, double an, double aw, double bn, double bw) 80 { 81 double an1,aw1,bn1,bw1,pn,pw; 82 double theta; 83 double cn,cw,cn1,cw1; 84 int i,n; 85 orient(an,aw,0.); 86 rotate(bn,bw,&bn1,&bw1); 87 /* printf("b %f %f\n",bn1,bw1);*/ 88 orient(an,aw,bw1); 89 rinvert(0.,dir,&pn,&pw); 90 /* printf("p %f %f\n",pn,pw);*/ 91 orient(pn,pw,0.); 92 rotate(an,aw,&an1,&aw1); 93 rotate(bn,bw,&bn1,&bw1); 94 theta = (aw1+bw1)/2; 95 /* printf("a %f %f \n",an1,aw1);*/ 96 orient(pn,pw,theta); 97 rotate(an,aw,&an1,&aw1); 98 rotate(bn,bw,&bn1,&bw1); 99 if(fabs(aw1-bw1)>180) 100 if(theta<0.) theta+=180; 101 else theta -= 180; 102 orient(pn,pw,theta); 103 rotate(an,aw,&an1,&aw1); 104 rotate(bn,bw,&bn1,&bw1); 105 if(!track) { 106 double dlat, dlon, t; 107 /* printf("A %.4f %.4f\n",an1,aw1); */ 108 /* printf("B %.4f %.4f\n",bn1,bw1); */ 109 cw1 = fabs(bw1-aw1); /* angular difference for map margins */ 110 /* while (aw<0.0) 111 aw += 360.; 112 while (bw<0.0) 113 bw += 360.; */ 114 dlon = fabs(aw-bw); 115 if (dlon>180) 116 dlon = 360-dlon; 117 dlat = fabs(an-bn); 118 printf("-o %.4f %.4f %.4f -w %.2f %.2f %.2f %.2f \n", 119 pn,pw,theta, -0.3*cw1, .3*cw1, -.6*cw1, .6*cw1); 120 121 } else { 122 cn1 = 0; 123 n = 1 + fabs(bw1-aw1)/.2; 124 for(i=0;i<=n;i++) { 125 cw1 = aw1 + i*(bw1-aw1)/n; 126 rinvert(cn1,cw1,&cn,&cw); 127 printf("%f %f\n",cn,cw); 128 } 129 printf("\"\n"); 130 } 131 }