1 #include <u.h>
2 #include <libc.h>
3 #include "map.h"
4
5 static struct coord p0; /* standard parallel */
6
7 int first;
8
9 static double
trigclamp(double x)10 trigclamp(double x)
11 {
12 return x>1? 1: x<-1? -1: x;
13 }
14
15 static struct coord az; /* azimuth of p0 seen from place */
16 static struct coord rad; /* angular dist from place to p0 */
17
18 static int
azimuth(struct place * place)19 azimuth(struct place *place)
20 {
21 if(place->nlat.c < FUZZ) {
22 az.l = PI/2 + place->nlat.l - place->wlon.l;
23 sincos(&az);
24 rad.l = fabs(place->nlat.l - p0.l);
25 if(rad.l > PI)
26 rad.l = 2*PI - rad.l;
27 sincos(&rad);
28 return 1;
29 }
30 rad.c = trigclamp(p0.s*place->nlat.s + /* law of cosines */
31 p0.c*place->nlat.c*place->wlon.c);
32 rad.s = sqrt(1 - rad.c*rad.c);
33 if(fabs(rad.s) < .001) {
34 az.s = 0;
35 az.c = 1;
36 } else {
37 az.s = trigclamp(p0.c*place->wlon.s/rad.s); /* sines */
38 az.c = trigclamp((p0.s - rad.c*place->nlat.s)
39 /(rad.s*place->nlat.c));
40 }
41 rad.l = atan2(rad.s, rad.c);
42 return 1;
43 }
44
45 static int
Xmecca(struct place * place,double * x,double * y)46 Xmecca(struct place *place, double *x, double *y)
47 {
48 if(!azimuth(place))
49 return 0;
50 *x = -place->wlon.l;
51 *y = fabs(az.s)<.02? -az.c*rad.s/p0.c: *x*az.c/az.s;
52 return fabs(*y)>2? -1:
53 rad.c<0? 0:
54 1;
55 }
56
57 proj
mecca(double par)58 mecca(double par)
59 {
60 first = 1;
61 if(fabs(par)>80.)
62 return(0);
63 deg2rad(par,&p0);
64 return(Xmecca);
65 }
66
67 static int
Xhoming(struct place * place,double * x,double * y)68 Xhoming(struct place *place, double *x, double *y)
69 {
70 if(!azimuth(place))
71 return 0;
72 *x = -rad.l*az.s;
73 *y = -rad.l*az.c;
74 return place->wlon.c<0? 0: 1;
75 }
76
77 proj
homing(double par)78 homing(double par)
79 {
80 first = 1;
81 if(fabs(par)>80.)
82 return(0);
83 deg2rad(par,&p0);
84 return(Xhoming);
85 }
86
87 int
hlimb(double * lat,double * lon,double res)88 hlimb(double *lat, double *lon, double res)
89 {
90 if(first) {
91 *lon = -90;
92 *lat = -90;
93 first = 0;
94 return 0;
95 }
96 *lat += res;
97 if(*lat <= 90)
98 return 1;
99 if(*lon == 90)
100 return -1;
101 *lon = 90;
102 *lat = -90;
103 return 0;
104 }
105
106 int
mlimb(double * lat,double * lon,double res)107 mlimb(double *lat, double *lon, double res)
108 {
109 int ret = !first;
110 if(fabs(p0.s) < .01)
111 return -1;
112 if(first) {
113 *lon = -180;
114 first = 0;
115 } else
116 *lon += res;
117 if(*lon > 180)
118 return -1;
119 *lat = atan(-cos(*lon*RAD)/p0.s*p0.c)/RAD;
120 return ret;
121 }
122