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