1 #include <u.h>
2 #include <libc.h>
3 #include <stdio.h>
4 #include "map.h"
5 
6 static double cirmod(double);
7 
8 static struct place pole;	/* map pole is tilted to here */
9 static struct coord twist;	/* then twisted this much */
10 static struct place ipole;	/* inverse transfrom */
11 static struct coord itwist;
12 
13 void
14 orient(double lat, double lon, double theta)
15 {
16 	lat = cirmod(lat);
17 	if(lat>90.) {
18 		lat = 180. - lat;
19 		lon -= 180.;
20 		theta -= 180.;
21 	} else if(lat < -90.) {
22 		lat = -180. - lat;
23 		lon -= 180.;
24 		theta -= 180;
25 	}
26 	latlon(lat,lon,&pole);
27 	deg2rad(theta, &twist);
28 	latlon(lat,180.-theta,&ipole);
29 	deg2rad(180.-lon, &itwist);
30 }
31 
32 void
33 latlon(double lat, double lon, struct place *p)
34 {
35 	lat = cirmod(lat);
36 	if(lat>90.) {
37 		lat = 180. - lat;
38 		lon -= 180.;
39 	} else if(lat < -90.) {
40 		lat = -180. - lat;
41 		lon -= 180.;
42 	}
43 	deg2rad(lat,&p->nlat);
44 	deg2rad(lon,&p->wlon);
45 }
46 
47 void
48 deg2rad(double theta, struct coord *coord)
49 {
50 	theta = cirmod(theta);
51 	coord->l = theta*RAD;
52 	if(theta==90) {
53 		coord->s = 1;
54 		coord->c = 0;
55 	} else if(theta== -90) {
56 		coord->s = -1;
57 		coord->c = 0;
58 	} else
59 		sincos(coord);
60 }
61 
62 static double
63 cirmod(double theta)
64 {
65 	while(theta >= 180.)
66 		theta -= 360;
67 	while(theta<-180.)
68 		theta += 360.;
69 	return(theta);
70 }
71 
72 void
73 sincos(struct coord *coord)
74 {
75 	coord->s = sin(coord->l);
76 	coord->c = cos(coord->l);
77 }
78 
79 void
80 normalize(struct place *gg)
81 {
82 	norm(gg,&pole,&twist);
83 }
84 
85 void
86 invert(struct place *g)
87 {
88 	norm(g,&ipole,&itwist);
89 }
90 
91 void
92 norm(struct place *gg, struct place *pp, struct coord *tw)
93 {
94 	register struct place *g;	/*geographic coords */
95 	register struct place *p;	/* new pole in old coords*/
96 	struct place m;			/* standard map coords*/
97 	g = gg;
98 	p = pp;
99 	if(p->nlat.s == 1.) {
100 		if(p->wlon.l+tw->l == 0.)
101 			return;
102 		g->wlon.l -= p->wlon.l+tw->l;
103 	} else {
104 		if(p->wlon.l != 0) {
105 			g->wlon.l -= p->wlon.l;
106 			sincos(&g->wlon);
107 		}
108 		m.nlat.s = p->nlat.s * g->nlat.s
109 			+ p->nlat.c * g->nlat.c * g->wlon.c;
110 		m.nlat.c = sqrt(1. - m.nlat.s * m.nlat.s);
111 		m.nlat.l = atan2(m.nlat.s, m.nlat.c);
112 		m.wlon.s = g->nlat.c * g->wlon.s;
113 		m.wlon.c = p->nlat.c * g->nlat.s
114 			- p->nlat.s * g->nlat.c * g->wlon.c;
115 		m.wlon.l = atan2(m.wlon.s, - m.wlon.c)
116 			- tw->l;
117 		*g = m;
118 	}
119 	sincos(&g->wlon);
120 	if(g->wlon.l>PI)
121 		g->wlon.l -= 2*PI;
122 	else if(g->wlon.l<-PI)
123 		g->wlon.l += 2*PI;
124 }
125 
126 void
127 printp(struct place *g)
128 {
129 printf("%.3f %.3f %.3f %.3f %.3f %.3f\n",
130 g->nlat.l,g->nlat.s,g->nlat.c,g->wlon.l,g->wlon.s,g->wlon.c);
131 }
132 
133 void
134 copyplace(struct place *g1, struct place *g2)
135 {
136 	*g2 = *g1;
137 }