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