summaryrefslogtreecommitdiffstats
path: root/kworldwatch/astro.c
diff options
context:
space:
mode:
authortoma <toma@283d02a7-25f6-0310-bc7c-ecb5cbfe19da>2009-11-25 17:56:58 +0000
committertoma <toma@283d02a7-25f6-0310-bc7c-ecb5cbfe19da>2009-11-25 17:56:58 +0000
commitae2a03c2941bf92573f89b88ef73f8aa842bea0a (patch)
tree3566563f3fb6ac3cb3496669d8f233062d3091bc /kworldwatch/astro.c
downloadtdetoys-ae2a03c2941bf92573f89b88ef73f8aa842bea0a.tar.gz
tdetoys-ae2a03c2941bf92573f89b88ef73f8aa842bea0a.zip
Copy the KDE 3.5 branch to branches/trinity for new KDE 3.5 features.
BUG:215923 git-svn-id: svn://anonsvn.kde.org/home/kde/branches/trinity/kdetoys@1054174 283d02a7-25f6-0310-bc7c-ecb5cbfe19da
Diffstat (limited to 'kworldwatch/astro.c')
-rw-r--r--kworldwatch/astro.c166
1 files changed, 166 insertions, 0 deletions
diff --git a/kworldwatch/astro.c b/kworldwatch/astro.c
new file mode 100644
index 0000000..3423119
--- /dev/null
+++ b/kworldwatch/astro.c
@@ -0,0 +1,166 @@
+/*
+ * Sun clock - astronomical routines.
+ */
+
+#include "sunclock.h"
+
+long jdate(struct tm *);
+double jtime(struct tm *);
+double kepler(double m, double ecc);
+void sunpos(double jd, int apparent, double *ra, double *dec, double *rv, double *slong);
+double gmst(double jd);
+
+/* JDATE -- Convert internal GMT date and time to Julian day
+ and fraction. */
+
+long
+jdate(t)
+struct tm *t;
+{
+ long c, m, y;
+
+ y = t->tm_year + 1900;
+ m = t->tm_mon + 1;
+ if (m > 2)
+ m = m - 3;
+ else {
+ m = m + 9;
+ y--;
+ }
+ c = y / 100L; /* Compute century */
+ y -= 100L * c;
+ return t->tm_mday + (c * 146097L) / 4 + (y * 1461L) / 4 +
+ (m * 153L + 2) / 5 + 1721119L;
+}
+
+/* JTIME -- Convert internal GMT date and time to astronomical
+ Julian time (i.e. Julian date plus day fraction,
+ expressed as a double). */
+
+double
+jtime(t)
+struct tm *t;
+{
+ return (jdate(t) - 0.5) +
+ (((long) t->tm_sec) +
+ 60L * (t->tm_min + 60L * t->tm_hour)) / 86400.0;
+}
+
+/* KEPLER -- Solve the equation of Kepler. */
+
+double
+kepler(m, ecc)
+double m, ecc;
+{
+ double e, delta;
+#define EPSILON 1E-6
+
+ e = m = dtr(m);
+ do {
+ delta = e - ecc * sin(e) - m;
+ e -= delta / (1 - ecc * cos(e));
+ } while (abs(delta) > EPSILON);
+ return e;
+}
+
+/* SUNPOS -- Calculate position of the Sun. JD is the Julian date
+ of the instant for which the position is desired and
+ APPARENT should be nonzero if the apparent position
+ (corrected for nutation and aberration) is desired.
+ The Sun's co-ordinates are returned in RA and DEC,
+ both specified in degrees (divide RA by 15 to obtain
+ hours). The radius vector to the Sun in astronomical
+ units is returned in RV and the Sun's longitude (true
+ or apparent, as desired) is returned as degrees in
+ SLONG. */
+
+void
+sunpos(jd, apparent, ra, dec, rv, slong)
+double jd;
+int apparent;
+double *ra, *dec, *rv, *slong;
+{
+ double t, t2, t3, l, m, e, ea, v, theta, omega,
+ eps;
+
+ /* Time, in Julian centuries of 36525 ephemeris days,
+ measured from the epoch 1900 January 0.5 ET. */
+
+ t = (jd - 2415020.0) / 36525.0;
+ t2 = t * t;
+ t3 = t2 * t;
+
+ /* Geometric mean longitude of the Sun, referred to the
+ mean equinox of the date. */
+
+ l = fixangle(279.69668 + 36000.76892 * t + 0.0003025 * t2);
+
+ /* Sun's mean anomaly. */
+
+ m = fixangle(358.47583 + 35999.04975*t - 0.000150*t2 - 0.0000033*t3);
+
+ /* Eccentricity of the Earth's orbit. */
+
+ e = 0.01675104 - 0.0000418 * t - 0.000000126 * t2;
+
+ /* Eccentric anomaly. */
+
+ ea = kepler(m, e);
+
+ /* True anomaly */
+
+ v = fixangle(2 * rtd(atan(sqrt((1 + e) / (1 - e)) * tan(ea / 2))));
+
+ /* Sun's true longitude. */
+
+ theta = l + v - m;
+
+ /* Obliquity of the ecliptic. */
+
+ eps = 23.452294 - 0.0130125 * t - 0.00000164 * t2 + 0.000000503 * t3;
+
+ /* Corrections for Sun's apparent longitude, if desired. */
+
+ if (apparent) {
+ omega = fixangle(259.18 - 1934.142 * t);
+ theta = theta - 0.00569 - 0.00479 * sin(dtr(omega));
+ eps += 0.00256 * cos(dtr(omega));
+ }
+
+ /* Return Sun's longitude and radius vector */
+
+ *slong = theta;
+ *rv = (1.0000002 * (1 - e * e)) / (1 + e * cos(dtr(v)));
+
+ /* Determine solar co-ordinates. */
+
+ *ra =
+ fixangle(rtd(atan2(cos(dtr(eps)) * sin(dtr(theta)), cos(dtr(theta)))));
+ *dec = rtd(asin(sin(dtr(eps)) * sin(dtr(theta))));
+}
+
+/* GMST -- Calculate Greenwich Mean Siderial Time for a given
+ instant expressed as a Julian date and fraction. */
+
+double
+gmst(jd)
+double jd;
+{
+ double t, theta0;
+
+
+ /* Time, in Julian centuries of 36525 ephemeris days,
+ measured from the epoch 1900 January 0.5 ET. */
+
+ t = ((floor(jd + 0.5) - 0.5) - 2415020.0) / 36525.0;
+
+ theta0 = 6.6460656 + 2400.051262 * t + 0.00002581 * t * t;
+
+ t = (jd + 0.5) - (floor(jd + 0.5));
+
+ theta0 += (t * 24.0) * 1.002737908;
+
+ theta0 = (theta0 - 24.0 * (floor(theta0 / 24.0)));
+
+ return theta0;
+}