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