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.
10 -w and -l options are suggested also.
12 Option -t prints out a series of
13 coordinates that follows the (great circle) track
14 in the original coordinate system,
16 This data is just right for map -t.
18 Option -i inverts the map top-to-bottom.
25 extern void doroute(double, double, double, double, double);
28 dorot(double a, double b, double *x, double *y, void (*f)(struct place *))
39 rotate(double a, double b, double *x, double *y)
41 dorot(a,b,x,y,normalize);
45 rinvert(double a, double b, double *x, double *y)
47 dorot(a,b,x,y,invert);
50 main(int argc, char **argv)
64 exits("route: bad option");
67 print("use route [-t] [-i] lat lon lat lon\n");
74 doroute(inv*90.,an,aw,bn,bw);
79 doroute(double dir, double an, double aw, double bn, double bw)
81 double an1,aw1,bn1,bw1,pn,pw;
86 rotate(bn,bw,&bn1,&bw1);
87 /* printf("b %f %f\n",bn1,bw1);*/
89 rinvert(0.,dir,&pn,&pw);
90 /* printf("p %f %f\n",pn,pw);*/
92 rotate(an,aw,&an1,&aw1);
93 rotate(bn,bw,&bn1,&bw1);
95 /* printf("a %f %f \n",an1,aw1);*/
97 rotate(an,aw,&an1,&aw1);
98 rotate(bn,bw,&bn1,&bw1);
100 if(theta<0.) theta+=180;
103 rotate(an,aw,&an1,&aw1);
104 rotate(bn,bw,&bn1,&bw1);
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 */
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);
123 n = 1 + fabs(bw1-aw1)/.2;
125 cw1 = aw1 + i*(bw1-aw1)/n;
126 rinvert(cn1,cw1,&cn,&cw);
127 printf("%f %f\n",cn,cw);