7 * uses lambda, beta, rad, motion
8 * sets alpha, delta, rp
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
22 * compute geocentric distance of object and
23 * compute light-time correction (i.i. planetary aberration)
26 xmp = rad*cos(beta)*cos(lambda);
27 ymp = rad*cos(beta)*sin(lambda);
29 rp = sqrt((xmp+xms)*(xmp+xms) +
32 lmb2 = lambda - .0057756e0*rp*motion;
34 xmp = rad*cos(beta)*cos(lmb2);
35 ymp = rad*cos(beta)*sin(lmb2);
39 * compute annual parallax from the position of the sun
45 rp = sqrt(xmp*xmp + ymp*ymp + zmp*zmp);
48 * compute annual (i.e. stellar) aberration
49 * from the orbital velocity of the earth
50 * (by an incorrect method)
58 * perform the nutation and so convert from the mean
59 * equator and equinox to the true
62 lmb2 = atan2(ymp, xmp);
63 beta2 = atan2(zmp, sqrt(xmp*xmp+ymp*ymp));
67 * change to equatorial coordinates
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));
74 alpha = atan2(ymp, xmp);
75 delta = atan2(zmp, sqrt(xmp*xmp+ymp*ymp));
77 hp = 8.794e0*radsec/rp;
79 if(rad > 0 && rad < 2.e5)
80 mag += 2.17*log(rad*rp);