Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 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
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
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
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
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
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
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
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
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 
}