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
dorot(double a,double b,double * x,double * y,void (* f)(struct place *))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
rotate(double a,double b,double * x,double * y)39 rotate(double a, double b, double *x, double *y)
40 {
41 dorot(a,b,x,y,normalize);
42 }
43
44 void
rinvert(double a,double b,double * x,double * y)45 rinvert(double a, double b, double *x, double *y)
46 {
47 dorot(a,b,x,y,invert);
48 }
49
main(int argc,char ** argv)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
doroute(double dir,double an,double aw,double bn,double bw)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 }
132