1 files changed, 198 insertions, 0 deletions
diff --git a/libmateweather/weather-moon.c b/libmateweather/weather-moon.c
new file mode 100644
index 0000000..99f0a25
--- /dev/null
+++ b/libmateweather/weather-moon.c
@@ -0,0 +1,198 @@
+/* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+/* weather-moon.c - Lunar calculations for mateweather
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation; either version 2 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, see
+ * <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * Formulas from:
+ * "Practical Astronomy With Your Calculator" (3e), Peter Duffett-Smith
+ * Cambridge University Press 1988
+ */
+
+#ifdef HAVE_CONFIG_H
+# include <config.h>
+#endif
+
+#ifdef __FreeBSD__
+#include <sys/types.h>
+#endif
+
+#include <math.h>
+#include <time.h>
+#include <string.h>
+#include <glib.h>
+
+#define MATEWEATHER_I_KNOW_THIS_IS_UNSTABLE
+#include "weather-priv.h"
+
+/*
+ * Elements of the Moon's matecorba, epoch 2000 Jan 1.5
+ * http://ssd.jpl.nasa.gov/?sat_elem#earth
+ * The page only lists most values to 2 decimal places
+ */
+
+#define LUNAR_MEAN_LONGITUDE 218.316
+#define LUNAR_PERIGEE_MEAN_LONG 318.15
+#define LUNAR_NODE_MEAN_LONG 125.08
+#define LUNAR_PROGRESSION 13.176358
+#define LUNAR_INCLINATION DEGREES_TO_RADIANS(5.145396)
+
+/**
+ * calc_moon:
+ * @info: WeatherInfo containing time_t of interest. The
+ * values moonphase, moonlatitude and moonValid are updated
+ * on success.
+ *
+ * Returns: gboolean indicating success or failure.
+ * moonphase is expressed as degrees where '0' is a new moon,
+ * '90' is first quarter, etc.
+ */
+
+gboolean
+calc_moon (WeatherInfo *info)
+{
+ time_t t;
+ gdouble ra_h;
+ gdouble decl_r;
+ gdouble ndays, sunMeanAnom_d;
+ gdouble moonLong_d;
+ gdouble moonMeanAnom_d, moonMeanAnom_r;
+ gdouble sunEclipLong_r;
+ gdouble ascNodeMeanLong_d;
+ gdouble corrLong_d, eviction_d;
+ gdouble sinSunMeanAnom;
+ gdouble Ae, A3, Ec, A4, lN_r;
+ gdouble lambda_r, beta_r;
+
+ /*
+ * The comments refer to the enumerated steps to calculate the
+ * position of the moon (section 65 of above reference)
+ */
+ t = info->update;
+ ndays = EPOCH_TO_J2000(t) / 86400.;
+ sunMeanAnom_d = fmod (MEAN_ECLIPTIC_LONGITUDE (ndays) - PERIGEE_LONGITUDE (ndays),
+ 360.);
+ sunEclipLong_r = sunEclipLongitude (t);
+ moonLong_d = fmod (LUNAR_MEAN_LONGITUDE + (ndays * LUNAR_PROGRESSION),
+ 360.);
+ /* 5: moon's mean anomaly */
+ moonMeanAnom_d = fmod ((moonLong_d - (0.1114041 * ndays)
+ - (LUNAR_PERIGEE_MEAN_LONG + LUNAR_NODE_MEAN_LONG)),
+ 360.);
+ /* 6: ascending node mean longitude */
+ ascNodeMeanLong_d = fmod (LUNAR_NODE_MEAN_LONG - (0.0529539 * ndays),
+ 360.);
+ eviction_d = 1.2739 /* 7: eviction */
+ * sin (DEGREES_TO_RADIANS (2.0 * (moonLong_d - RADIANS_TO_DEGREES (sunEclipLong_r))
+ - moonMeanAnom_d));
+ sinSunMeanAnom = sin (DEGREES_TO_RADIANS (sunMeanAnom_d));
+ Ae = 0.1858 * sinSunMeanAnom;
+ A3 = 0.37 * sinSunMeanAnom; /* 8: annual equation */
+ moonMeanAnom_d += eviction_d - Ae - A3; /* 9: "third correction" */
+ moonMeanAnom_r = DEGREES_TO_RADIANS (moonMeanAnom_d);
+ Ec = 6.2886 * sin (moonMeanAnom_r); /* 10: equation of center */
+ A4 = 0.214 * sin (2.0 * moonMeanAnom_r); /* 11: "yet another correction" */
+
+ /* Steps 12-14 give the true longitude after correcting for variation */
+ moonLong_d += eviction_d + Ec - Ae + A4
+ + (0.6583 * sin (2.0 * (moonMeanAnom_r - sunEclipLong_r)));
+
+ /* 15: corrected longitude of node */
+ corrLong_d = ascNodeMeanLong_d - 0.16 * sinSunMeanAnom;
+
+ /*
+ * Calculate ecliptic latitude (16-19) and longitude (20) of the moon,
+ * then convert to right ascension and declination.
+ */
+ lN_r = DEGREES_TO_RADIANS (moonLong_d - corrLong_d); /* l''-N' */
+ lambda_r = DEGREES_TO_RADIANS(corrLong_d)
+ + atan2 (sin (lN_r) * cos (LUNAR_INCLINATION), cos (lN_r));
+ beta_r = asin (sin (lN_r) * sin (LUNAR_INCLINATION));
+ ecl2equ (t, lambda_r, beta_r, &ra_h, &decl_r);
+
+ /*
+ * The phase is the angle from the sun's longitude to the moon's 
+ */
+ info->moonphase =
+ fmod (15.*ra_h - RADIANS_TO_DEGREES (sunEclipLongitude (info->update)),
+ 360.);
+ if (info->moonphase < 0)
+ info->moonphase += 360;
+ info->moonlatitude = RADIANS_TO_DEGREES (decl_r);
+ info->moonValid = TRUE;
+
+ return TRUE;
+}
+
+
+/**
+ * calc_moon_phases:
+ * @info: WeatherInfo containing the time_t of interest
+ * @phases: An array of four time_t values that will hold the returned values.
+ * The values are estimates of the time of the next new, quarter, full and
+ * three-quarter moons.
+ *
+ * Returns: gboolean indicating success or failure
+ */
+
+gboolean
+calc_moon_phases (WeatherInfo *info, time_t *phases)
+{
+ WeatherInfo temp;
+ time_t *ptime;
+ int idx;
+ gdouble advance;
+ int iter;
+ time_t dtime;
+
+ g_return_val_if_fail (info != NULL &&
+ (info->moonValid || calc_moon (info)),
+ FALSE);
+
+ ptime = phases;
+ memset(&temp, 0, sizeof(WeatherInfo));
+
+ for (idx = 0; idx < 4; idx++) {
+ temp.update = info->update;
+ temp.moonphase = info->moonphase;
+
+ /*
+ * First estimate on how far the moon needs to advance
+ * to get to the required phase
+ */
+ advance = (idx * 90.) - info->moonphase;
+ if (advance < 0.)
+ advance += 360.;
+
+ for (iter = 0; iter < 10; iter++) {
+ /* Convert angle change (degrees) to dtime (seconds) */
+ dtime = advance / LUNAR_PROGRESSION * 86400.;
+ if ((dtime > -10) && (dtime < 10))
+ break;
+ temp.update += dtime;
+ (void)calc_moon (&temp);
+
+ if (idx == 0 && temp.moonphase > 180.) {
+ advance = 360. - temp.moonphase;
+ } else {
+ advance = (idx * 90.) - temp.moonphase;
+ }
+ }
+ *ptime++ = temp.update;
+ }
+
+ return TRUE;
+}