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