summaryrefslogtreecommitdiff
path: root/sys/src/cmd/astro/helio.c
diff options
context:
space:
mode:
authorTaru Karttunen <taruti@taruti.net>2011-03-30 15:46:40 +0300
committerTaru Karttunen <taruti@taruti.net>2011-03-30 15:46:40 +0300
commite5888a1ffdae813d7575f5fb02275c6bb07e5199 (patch)
treed8d51eac403f07814b9e936eed0c9a79195e2450 /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-xsys/src/cmd/astro/helio.c81
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);
+}