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