1*3e12c5d1SDavid du Colombier #include "astro.h"
2*3e12c5d1SDavid du Colombier
3*3e12c5d1SDavid du Colombier void
helio(void)4*3e12c5d1SDavid du Colombier helio(void)
5*3e12c5d1SDavid du Colombier {
6*3e12c5d1SDavid du Colombier /*
7*3e12c5d1SDavid du Colombier * uses lambda, beta, rad, motion
8*3e12c5d1SDavid du Colombier * sets alpha, delta, rp
9*3e12c5d1SDavid du Colombier */
10*3e12c5d1SDavid du Colombier
11*3e12c5d1SDavid du Colombier /*
12*3e12c5d1SDavid du Colombier * helio converts from ecliptic heliocentric coordinates
13*3e12c5d1SDavid du Colombier * referred to the mean equinox of date
14*3e12c5d1SDavid du Colombier * to equatorial geocentric coordinates referred to
15*3e12c5d1SDavid du Colombier * the true equator and equinox
16*3e12c5d1SDavid du Colombier */
17*3e12c5d1SDavid du Colombier
18*3e12c5d1SDavid du Colombier double xmp, ymp, zmp;
19*3e12c5d1SDavid du Colombier double beta2;
20*3e12c5d1SDavid du Colombier
21*3e12c5d1SDavid du Colombier /*
22*3e12c5d1SDavid du Colombier * compute geocentric distance of object and
23*3e12c5d1SDavid du Colombier * compute light-time correction (i.i. planetary aberration)
24*3e12c5d1SDavid du Colombier */
25*3e12c5d1SDavid du Colombier
26*3e12c5d1SDavid du Colombier xmp = rad*cos(beta)*cos(lambda);
27*3e12c5d1SDavid du Colombier ymp = rad*cos(beta)*sin(lambda);
28*3e12c5d1SDavid du Colombier zmp = rad*sin(beta);
29*3e12c5d1SDavid du Colombier rp = sqrt((xmp+xms)*(xmp+xms) +
30*3e12c5d1SDavid du Colombier (ymp+yms)*(ymp+yms) +
31*3e12c5d1SDavid du Colombier (zmp+zms)*(zmp+zms));
32*3e12c5d1SDavid du Colombier lmb2 = lambda - .0057756e0*rp*motion;
33*3e12c5d1SDavid du Colombier
34*3e12c5d1SDavid du Colombier xmp = rad*cos(beta)*cos(lmb2);
35*3e12c5d1SDavid du Colombier ymp = rad*cos(beta)*sin(lmb2);
36*3e12c5d1SDavid du Colombier zmp = rad*sin(beta);
37*3e12c5d1SDavid du Colombier
38*3e12c5d1SDavid du Colombier /*
39*3e12c5d1SDavid du Colombier * compute annual parallax from the position of the sun
40*3e12c5d1SDavid du Colombier */
41*3e12c5d1SDavid du Colombier
42*3e12c5d1SDavid du Colombier xmp += xms;
43*3e12c5d1SDavid du Colombier ymp += yms;
44*3e12c5d1SDavid du Colombier zmp += zms;
45*3e12c5d1SDavid du Colombier rp = sqrt(xmp*xmp + ymp*ymp + zmp*zmp);
46*3e12c5d1SDavid du Colombier
47*3e12c5d1SDavid du Colombier /*
48*3e12c5d1SDavid du Colombier * compute annual (i.e. stellar) aberration
49*3e12c5d1SDavid du Colombier * from the orbital velocity of the earth
50*3e12c5d1SDavid du Colombier * (by an incorrect method)
51*3e12c5d1SDavid du Colombier */
52*3e12c5d1SDavid du Colombier
53*3e12c5d1SDavid du Colombier xmp -= xdot*rp;
54*3e12c5d1SDavid du Colombier ymp -= ydot*rp;
55*3e12c5d1SDavid du Colombier zmp -= zdot*rp;
56*3e12c5d1SDavid du Colombier
57*3e12c5d1SDavid du Colombier /*
58*3e12c5d1SDavid du Colombier * perform the nutation and so convert from the mean
59*3e12c5d1SDavid du Colombier * equator and equinox to the true
60*3e12c5d1SDavid du Colombier */
61*3e12c5d1SDavid du Colombier
62*3e12c5d1SDavid du Colombier lmb2 = atan2(ymp, xmp);
63*3e12c5d1SDavid du Colombier beta2 = atan2(zmp, sqrt(xmp*xmp+ymp*ymp));
64*3e12c5d1SDavid du Colombier lmb2 += phi;
65*3e12c5d1SDavid du Colombier
66*3e12c5d1SDavid du Colombier /*
67*3e12c5d1SDavid du Colombier * change to equatorial coordinates
68*3e12c5d1SDavid du Colombier */
69*3e12c5d1SDavid du Colombier
70*3e12c5d1SDavid du Colombier xmp = rp*cos(lmb2)*cos(beta2);
71*3e12c5d1SDavid du Colombier ymp = rp*(sin(lmb2)*cos(beta2)*cos(tobliq) - sin(tobliq)*sin(beta2));
72*3e12c5d1SDavid du Colombier zmp = rp*(sin(lmb2)*cos(beta2)*sin(tobliq) + cos(tobliq)*sin(beta2));
73*3e12c5d1SDavid du Colombier
74*3e12c5d1SDavid du Colombier alpha = atan2(ymp, xmp);
75*3e12c5d1SDavid du Colombier delta = atan2(zmp, sqrt(xmp*xmp+ymp*ymp));
76*3e12c5d1SDavid du Colombier
77*3e12c5d1SDavid du Colombier hp = 8.794e0*radsec/rp;
78*3e12c5d1SDavid du Colombier semi /= rp;
79*3e12c5d1SDavid du Colombier if(rad > 0 && rad < 2.e5)
80*3e12c5d1SDavid du Colombier mag += 2.17*log(rad*rp);
81*3e12c5d1SDavid du Colombier }
82