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
orient(double lat,double lon,double theta)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
latlon(double lat,double lon,struct place * p)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
deg2rad(double theta,struct coord * coord)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
cirmod(double theta)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
sincos(struct coord * coord)73 sincos(struct coord *coord)
74 {
75 coord->s = sin(coord->l);
76 coord->c = cos(coord->l);
77 }
78
79 void
normalize(struct place * gg)80 normalize(struct place *gg)
81 {
82 norm(gg,&pole,&twist);
83 }
84
85 void
invert(struct place * g)86 invert(struct place *g)
87 {
88 norm(g,&ipole,&itwist);
89 }
90
91 void
norm(struct place * gg,struct place * pp,struct coord * tw)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 double
tan(double x)127 tan(double x)
128 {
129 return(sin(x)/cos(x));
130 }
131
132 void
printp(struct place * g)133 printp(struct place *g)
134 {
135 printf("%.3f %.3f %.3f %.3f %.3f %.3f\n",
136 g->nlat.l,g->nlat.s,g->nlat.c,g->wlon.l,g->wlon.s,g->wlon.c);
137 }
138
139 void
copyplace(struct place * g1,struct place * g2)140 copyplace(struct place *g1, struct place *g2)
141 {
142 *g2 = *g1;
143 }
144