#include <u.h>#include <libc.h>#include "map.h" /* Given two lat-lon pairs, find an orientation for the   -o option of "map" that will place those two points   on the equator of a standard projection, equally spaced   about the prime meridian.    -w and -l options are suggested also.    Option -t prints out a series of   coordinates that follows the (great circle) track   in the original coordinate system,   followed by ".   This data is just right for map -t.    Option -i inverts the map top-to-bottom.*/struct place pole;struct coord twist;int track;int inv = -1; extern void doroute(double, double, double, double, double); voiddorot(double a, double b, double *x, double *y, void (*f)(struct place *)){	struct place g;	deg2rad(a,&g.nlat);	deg2rad(b,&g.wlon);	(*f)(&g);	*x = g.nlat.l/RAD;	*y = g.wlon.l/RAD;} voidrotate(double a, double b, double *x, double *y){	dorot(a,b,x,y,normalize);} voidrinvert(double a, double b, double *x, double *y){	dorot(a,b,x,y,invert);} main(int argc, char **argv){#pragma ref argv	double an,aw,bn,bw;	ARGBEGIN {	case 't':		track = 1;		break; 	case 'i':		inv = 1;		break; 	default:		exits("route: bad option");	} ARGEND;	if (argc<4) {		print("use route [-t] [-i] lat lon lat lon\n");		exits("arg count");	}	an = atof(argv[0]);	aw = atof(argv[1]);	bn = atof(argv[2]);	bw = atof(argv[3]);	doroute(inv*90.,an,aw,bn,bw);	return 0;} voiddoroute(double dir, double an, double aw, double bn, double bw){	double an1,aw1,bn1,bw1,pn,pw;	double theta;	double cn,cw,cn1,cw1;	int i,n;	orient(an,aw,0.);	rotate(bn,bw,&bn1,&bw1);/*	printf("b %f %f\n",bn1,bw1);*/	orient(an,aw,bw1);	rinvert(0.,dir,&pn,&pw);/*	printf("p %f %f\n",pn,pw);*/	orient(pn,pw,0.);	rotate(an,aw,&an1,&aw1);	rotate(bn,bw,&bn1,&bw1);	theta = (aw1+bw1)/2;/*	printf("a %f %f \n",an1,aw1);*/	orient(pn,pw,theta);	rotate(an,aw,&an1,&aw1);	rotate(bn,bw,&bn1,&bw1);	if(fabs(aw1-bw1)>180)		if(theta<0.) theta+=180;		else theta -= 180;	orient(pn,pw,theta);	rotate(an,aw,&an1,&aw1);	rotate(bn,bw,&bn1,&bw1);	if(!track) {		double dlat, dlon, t;		/* printf("A %.4f %.4f\n",an1,aw1); */		/* printf("B %.4f %.4f\n",bn1,bw1); */		cw1 = fabs(bw1-aw1);	/* angular difference for map margins */		/* while (aw<0.0)			aw += 360.;		while (bw<0.0)			bw += 360.; */		dlon = fabs(aw-bw);		if (dlon>180)			dlon = 360-dlon;		dlat = fabs(an-bn);		printf("-o %.4f %.4f %.4f -w %.2f %.2f %.2f %.2f \n",		  pn,pw,theta, -0.3*cw1, .3*cw1, -.6*cw1, .6*cw1);			} else {		cn1 = 0;		n = 1 + fabs(bw1-aw1)/.2;		for(i=0;i<=n;i++) {			cw1 = aw1 + i*(bw1-aw1)/n;			rinvert(cn1,cw1,&cn,&cw);			printf("%f %f\n",cn,cw);		}		printf("\"\n");	}}