Subversion Repositories planix.SVN

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 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 
}