123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140 |
- /*
- * This file is part of the UCB release of Plan 9. It is subject to the license
- * terms in the LICENSE file found in the top-level directory of this
- * distribution and at http://akaros.cs.berkeley.edu/files/Plan9License. No
- * part of the UCB release of Plan 9, including this file, may be copied,
- * modified, propagated, or distributed except according to the terms contained
- * in the LICENSE file.
- */
- #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);
- void
- dorot(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;
- }
- void
- rotate(double a, double b, double *x, double *y)
- {
- dorot(a,b,x,y,normalize);
- }
- void
- rinvert(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;
- }
- void
- doroute(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");
- }
- }
|