diff options
author | Taru Karttunen <taruti@taruti.net> | 2011-03-30 15:46:40 +0300 |
---|---|---|
committer | Taru Karttunen <taruti@taruti.net> | 2011-03-30 15:46:40 +0300 |
commit | e5888a1ffdae813d7575f5fb02275c6bb07e5199 (patch) | |
tree | d8d51eac403f07814b9e936eed0c9a79195e2450 /sys/src/cmd/astro/helio.c |
Import sources from 2011-03-30 iso image
Diffstat (limited to 'sys/src/cmd/astro/helio.c')
-rwxr-xr-x | sys/src/cmd/astro/helio.c | 81 |
1 files changed, 81 insertions, 0 deletions
diff --git a/sys/src/cmd/astro/helio.c b/sys/src/cmd/astro/helio.c new file mode 100755 index 000000000..cb4157863 --- /dev/null +++ b/sys/src/cmd/astro/helio.c @@ -0,0 +1,81 @@ +#include "astro.h" + +void +helio(void) +{ +/* + * uses lambda, beta, rad, motion + * sets alpha, delta, rp + */ + +/* + * helio converts from ecliptic heliocentric coordinates + * referred to the mean equinox of date + * to equatorial geocentric coordinates referred to + * the true equator and equinox + */ + + double xmp, ymp, zmp; + double beta2; + +/* + * compute geocentric distance of object and + * compute light-time correction (i.i. planetary aberration) + */ + + xmp = rad*cos(beta)*cos(lambda); + ymp = rad*cos(beta)*sin(lambda); + zmp = rad*sin(beta); + rp = sqrt((xmp+xms)*(xmp+xms) + + (ymp+yms)*(ymp+yms) + + (zmp+zms)*(zmp+zms)); + lmb2 = lambda - .0057756e0*rp*motion; + + xmp = rad*cos(beta)*cos(lmb2); + ymp = rad*cos(beta)*sin(lmb2); + zmp = rad*sin(beta); + +/* + * compute annual parallax from the position of the sun + */ + + xmp += xms; + ymp += yms; + zmp += zms; + rp = sqrt(xmp*xmp + ymp*ymp + zmp*zmp); + +/* + * compute annual (i.e. stellar) aberration + * from the orbital velocity of the earth + * (by an incorrect method) + */ + + xmp -= xdot*rp; + ymp -= ydot*rp; + zmp -= zdot*rp; + +/* + * perform the nutation and so convert from the mean + * equator and equinox to the true + */ + + lmb2 = atan2(ymp, xmp); + beta2 = atan2(zmp, sqrt(xmp*xmp+ymp*ymp)); + lmb2 += phi; + +/* + * change to equatorial coordinates + */ + + xmp = rp*cos(lmb2)*cos(beta2); + ymp = rp*(sin(lmb2)*cos(beta2)*cos(tobliq) - sin(tobliq)*sin(beta2)); + zmp = rp*(sin(lmb2)*cos(beta2)*sin(tobliq) + cos(tobliq)*sin(beta2)); + + alpha = atan2(ymp, xmp); + delta = atan2(zmp, sqrt(xmp*xmp+ymp*ymp)); + + hp = 8.794e0*radsec/rp; + semi /= rp; + if(rad > 0 && rad < 2.e5) + mag += 2.17*log(rad*rp); +} |