1 cd5bae78 2004-04-21 devnull #include "astro.h"
 2 cd5bae78 2004-04-21 devnull 
 3 cd5bae78 2004-04-21 devnull void
 4 cd5bae78 2004-04-21 devnull geo(void)
 5 cd5bae78 2004-04-21 devnull {
 6 cd5bae78 2004-04-21 devnull 
 7 cd5bae78 2004-04-21 devnull /*
 8 cd5bae78 2004-04-21 devnull  *	uses alpha, delta, rp
 9 cd5bae78 2004-04-21 devnull  */
10 cd5bae78 2004-04-21 devnull 
11 cd5bae78 2004-04-21 devnull /*
12 cd5bae78 2004-04-21 devnull  *	sets ra, decl, lha, decl2, az, el
13 cd5bae78 2004-04-21 devnull  */
14 cd5bae78 2004-04-21 devnull 
15 cd5bae78 2004-04-21 devnull /*
16 cd5bae78 2004-04-21 devnull  *	geo converts geocentric equatorial coordinates
17 cd5bae78 2004-04-21 devnull  *	to topocentric equatorial and topocentric horizon
18 cd5bae78 2004-04-21 devnull  *	coordinates.
19 cd5bae78 2004-04-21 devnull  *	All are (usually) referred to the true equator.
20 cd5bae78 2004-04-21 devnull  */
21 cd5bae78 2004-04-21 devnull 
22 cd5bae78 2004-04-21 devnull 	double sel, saz, caz;
23 cd5bae78 2004-04-21 devnull 	double f;
24 cd5bae78 2004-04-21 devnull 	double sa, ca, sd;
25 cd5bae78 2004-04-21 devnull 
26 cd5bae78 2004-04-21 devnull /*
27 cd5bae78 2004-04-21 devnull  *	convert to local hour angle and declination
28 cd5bae78 2004-04-21 devnull  */
29 cd5bae78 2004-04-21 devnull 
30 cd5bae78 2004-04-21 devnull 	lha = gst - alpha - wlong;
31 cd5bae78 2004-04-21 devnull 	decl = delta;
32 cd5bae78 2004-04-21 devnull 
33 cd5bae78 2004-04-21 devnull /*
34 cd5bae78 2004-04-21 devnull  *	compute diurnal parallax (requires geocentric latitude)
35 cd5bae78 2004-04-21 devnull  */
36 cd5bae78 2004-04-21 devnull 
37 cd5bae78 2004-04-21 devnull 	sa = cos(decl)*sin(lha);
38 cd5bae78 2004-04-21 devnull 	ca = cos(decl)*cos(lha) - erad*cos(glat)*sin(hp);
39 cd5bae78 2004-04-21 devnull 	sd = sin(decl)           - erad*sin(glat)*sin(hp);
40 cd5bae78 2004-04-21 devnull 
41 cd5bae78 2004-04-21 devnull 	lha = atan2(sa, ca);
42 cd5bae78 2004-04-21 devnull 	decl2 = atan2(sd, sqrt(sa*sa+ca*ca));
43 cd5bae78 2004-04-21 devnull 	f = sqrt(sa*sa+ca*ca+sd*sd);
44 cd5bae78 2004-04-21 devnull 	semi2 = semi/f;
45 cd5bae78 2004-04-21 devnull 	ra = gst - lha - wlong;
46 cd5bae78 2004-04-21 devnull 	ra = pinorm(ra);
47 cd5bae78 2004-04-21 devnull 
48 cd5bae78 2004-04-21 devnull /*
49 cd5bae78 2004-04-21 devnull  *	convert to horizon coordinates
50 cd5bae78 2004-04-21 devnull  */
51 cd5bae78 2004-04-21 devnull 
52 cd5bae78 2004-04-21 devnull 	sel = sin(nlat)*sin(decl2) + cos(nlat)*cos(decl2)*cos(lha);
53 cd5bae78 2004-04-21 devnull 	el = atan2(sel, pyth(sel));
54 cd5bae78 2004-04-21 devnull 	saz = sin(lha)*cos(decl2);
55 cd5bae78 2004-04-21 devnull 	caz = cos(nlat)*sin(decl2) - sin(nlat)*cos(decl2)*cos(lha);
56 cd5bae78 2004-04-21 devnull 	az = pi + atan2(saz, -caz);
57 cd5bae78 2004-04-21 devnull 
58 cd5bae78 2004-04-21 devnull 	az /= radian;
59 cd5bae78 2004-04-21 devnull 	el /= radian;
60 cd5bae78 2004-04-21 devnull }