1 files changed, 353 insertions, 0 deletions
diff --git a/libmateweather/weather-sun.c b/libmateweather/weather-sun.c
new file mode 100644
index 0000000..523bf0c
--- /dev/null
+++ b/libmateweather/weather-sun.c
@@ -0,0 +1,353 @@
+/* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+/* weather-sun.c - Astronomy 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
+ * Unless otherwise noted, comments referencing "steps" are related to
+ * the algorithm presented in section 49 of above
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <time.h>
+#include <glib.h>
+
+#define MATEWEATHER_I_KNOW_THIS_IS_UNSTABLE
+#include "weather-priv.h"
+
+#define ECCENTRICITY(d) (0.01671123 - (d)/36525.*0.00004392)
+
+/* 
+ * Ecliptic longitude of the sun at specified time (UT)
+ * The algoithm is described in section 47 of Duffett-Smith
+ * Return value is in radians
+ */
+gdouble
+sunEclipLongitude(time_t t)
+{
+ gdouble ndays, meanAnom, eccenAnom, delta, e, longitude;
+
+ /*
+ * Start with an estimate based on a fixed daily rate
+ */
+ ndays = EPOCH_TO_J2000(t) / 86400.;
+ meanAnom = DEGREES_TO_RADIANS(MEAN_ECLIPTIC_LONGITUDE(ndays)
+ - PERIGEE_LONGITUDE(ndays));
+ 
+ /*
+ * Approximate solution of Kepler's equation:
+ * Find E which satisfies E - e sin(E) = M (mean anomaly)
+ */ 
+ eccenAnom = meanAnom;
+ e = ECCENTRICITY(ndays);
+
+ while (1e-12 < fabs( delta = eccenAnom - e * sin(eccenAnom) - meanAnom))
+ {
+ eccenAnom -= delta / (1.- e * cos(eccenAnom));
+ }
+
+ /*
+ * Earth's longitude on the ecliptic
+ */
+ longitude = fmod( DEGREES_TO_RADIANS (PERIGEE_LONGITUDE(ndays))
+ + 2. * atan (sqrt ((1.+e)/(1.-e))
+ * tan (eccenAnom / 2.)),
+ 2. * M_PI);
+ if (longitude < 0.) {
+ longitude += 2 * M_PI;
+ }
+ return longitude;
+}
+
+static gdouble
+ecliptic_obliquity (gdouble time)
+{
+ gdouble jc = EPOCH_TO_J2000 (time) / (36525. * 86400.);
+ gdouble eclip_secs = (84381.448
+ - (46.84024 * jc)
+ - (59.e-5 * jc * jc)
+ + (1.813e-3 * jc * jc * jc));
+ return DEGREES_TO_RADIANS(eclip_secs / 3600.);
+}
+
+/*
+ * Convert ecliptic longitude and latitude (radians) to equitorial
+ * coordinates, expressed as right ascension (hours) and
+ * declination (radians)
+ */
+void
+ecl2equ (gdouble time,
+ gdouble eclipLon, gdouble eclipLat,
+ gdouble *ra, gdouble *decl)
+{
+ gdouble mEclipObliq = ecliptic_obliquity(time);
+
+ if (ra) {
+ *ra = RADIANS_TO_HOURS (atan2 ((sin (eclipLon) * cos (mEclipObliq)
+ - tan (eclipLat) * sin(mEclipObliq)),
+ cos (eclipLon)));
+ if (*ra < 0.)
+ *ra += 24.;
+ }
+ if (decl) {
+ *decl = asin (( sin (eclipLat) * cos (mEclipObliq))
+ + cos (eclipLat) * sin (mEclipObliq) * sin(eclipLon));
+ }
+}
+
+/*
+ * Calculate rising and setting times for an object
+ * based on it equitorial coordinates (section 33 & 15)
+ * Returned "rise" and "set" values are sideral times in hours
+ */
+static void
+gstObsv (gdouble ra, gdouble decl,
+ gdouble obsLat, gdouble obsLon,
+ gdouble *rise, gdouble *set)
+{
+ double a = acos (-tan (obsLat) * tan (decl));
+ double b;
+
+ if (isnan (a) != 0) {
+ *set = *rise = a;
+ return;
+ }
+ a = RADIANS_TO_HOURS (a);
+ b = 24. - a + ra;
+ a += ra;
+ a -= RADIANS_TO_HOURS (obsLon);
+ b -= RADIANS_TO_HOURS (obsLon);
+ if ((a = fmod (a, 24.)) < 0)
+ a += 24.;
+ if ((b = fmod (b, 24.)) < 0)
+ b += 24.;
+
+ *set = a;
+ *rise = b;
+}
+
+
+static gdouble
+t0 (time_t date)
+{
+ gdouble t = ((gdouble)(EPOCH_TO_J2000 (date) / 86400)) / 36525.0;
+ gdouble t0 = fmod (6.697374558 + 2400.051366 * t + 2.5862e-5 * t * t, 24.);
+ if (t0 < 0.)
+ t0 += 24.;
+ return t0;
+}
+
+
+static gboolean
+calc_sun2 (WeatherInfo *info, time_t t)
+{
+ gdouble obsLat = info->location->latitude;
+ gdouble obsLon = info->location->longitude;
+ time_t gm_midn;
+ time_t lcl_midn;
+ gdouble gm_hoff, lambda;
+ gdouble ra1, ra2;
+ gdouble decl1, decl2;
+ gdouble decl_midn, decl_noon;
+ gdouble rise1, rise2;
+ gdouble set1, set2;
+ gdouble tt, t00;
+ gdouble x, u, dt;
+
+ /* Approximate preceding local midnight at observer's longitude */
+ obsLat = info->location->latitude;
+ obsLon = info->location->longitude;
+ gm_midn = t - (t % 86400);
+ gm_hoff = floor ((RADIANS_TO_DEGREES (obsLon) + 7.5) / 15.);
+ lcl_midn = gm_midn - 3600. * gm_hoff;
+ if (t - lcl_midn >= 86400)
+ lcl_midn += 86400;
+ else if (lcl_midn > t)
+ lcl_midn -= 86400;
+
+ lambda = sunEclipLongitude (lcl_midn);
+ 
+ /*
+ * Calculate equitorial coordinates of sun at previous
+ * and next local midnights
+ */
+ ecl2equ (lcl_midn, lambda, 0., &ra1, &decl1);
+ ecl2equ (lcl_midn + 86400.,
+ lambda + DEGREES_TO_RADIANS(SOL_PROGRESSION), 0.,
+ &ra2, &decl2);
+ 
+ /*
+ * If the observer is within the Arctic or Antarctic Circles then
+ * the sun may be above or below the horizon for the full day.
+ */
+ decl_midn = MIN(decl1,decl2);
+ decl_noon = (decl1+decl2)/2.;
+ info->midnightSun =
+ (obsLat > (M_PI/2.-decl_midn)) || (obsLat < (-M_PI/2.-decl_midn));
+ info->polarNight =
+ (obsLat > (M_PI/2.+decl_noon)) || (obsLat < (-M_PI/2.+decl_noon));
+ if (info->midnightSun || info->polarNight) {
+ info->sunriseValid = info->sunsetValid = FALSE;
+ return FALSE;
+ }
+
+ /*
+ * Convert to rise and set times based positions for the preceding
+ * and following local midnights.
+ */
+ gstObsv (ra1, decl1, obsLat, obsLon - (gm_hoff * M_PI / 12.), &rise1, &set1);
+ gstObsv (ra2, decl2, obsLat, obsLon - (gm_hoff * M_PI / 12.), &rise2, &set2);
+
+ /* TODO: include calculations for regions near the poles. */
+ if (isnan(rise1) || isnan(rise2)) {
+ info->sunriseValid = info->sunsetValid = FALSE;
+ return FALSE;
+ }
+
+ if (rise2 < rise1) {
+ rise2 += 24.;
+ }
+ if (set2 < set1) {
+ set2 += 24.;
+ }
+
+ tt = t0(lcl_midn);
+ t00 = tt - (gm_hoff + RADIANS_TO_HOURS(obsLon)) * 1.002737909;
+
+ if (t00 < 0.)
+ t00 += 24.;
+
+ if (rise1 < t00) {
+ rise1 += 24.;
+ rise2 += 24.;
+ }
+ if (set1 < t00) {
+ set1 += 24.;
+ set2 += 24.;
+ }
+ 
+ /*
+ * Interpolate between the two to get a rise and set time
+ * based on the sun's position at local noon (step 8)
+ */
+ rise1 = (24.07 * rise1 - t00 * (rise2 - rise1)) / (24.07 + rise1 - rise2);
+ set1 = (24.07 * set1 - t00 * (set2 - set1)) / (24.07 + set1 - set2);
+
+ /*
+ * Calculate an adjustment value to account for parallax,
+ * refraction and the Sun's finite diameter (steps 9,10)
+ */
+ decl2 = (decl1 + decl2) / 2.;
+ x = DEGREES_TO_RADIANS(0.830725);
+ u = acos ( sin(obsLat) / cos(decl2) );
+ dt = RADIANS_TO_HOURS ( asin ( sin(x) / sin(u) ) / cos(decl2) );
+ 
+ /*
+ * Subtract the correction value from sunrise and add to sunset,
+ * then (step 11) convert sideral times to UT
+ */
+ rise1 = (rise1 - dt - tt) * 0.9972695661;
+ if (rise1 < 0.)
+ rise1 += 24;
+ else if (rise1 >= 24.)
+ rise1 -= 24.;
+ info->sunriseValid = ((rise1 >= 0.) && (rise1 < 24.));
+ info->sunrise = (rise1 * 3600.) + lcl_midn;
+
+ set1 = (set1 + dt - tt) * 0.9972695661;
+ if (set1 < 0.)
+ set1 += 24;
+ else if (set1 >= 24.)
+ set1 -= 24.;
+ info->sunsetValid = ((set1 >= 0.) && (set1 < 24.));
+ info->sunset = (set1 * 3600.) + lcl_midn;
+
+ return (info->sunriseValid || info->sunsetValid);
+}
+
+
+/**
+ * calc_sun_time:
+ * @info: #WeatherInfo structure containing the observer's latitude
+ * and longitude in radians, fills in the sunrise and sunset times.
+ * @t: time_t
+ *
+ * Returns: gboolean indicating if the results are valid.
+ */
+gboolean
+calc_sun_time (WeatherInfo *info, time_t t)
+{
+ return info->location->latlon_valid && calc_sun2 (info, t);
+}
+
+/**
+ * calc_sun:
+ * @info: #WeatherInfo structure containing the observer's latitude
+ * and longitude in radians, fills in the sunrise and sunset times.
+ *
+ * Returns: gboolean indicating if the results are valid.
+ */
+gboolean
+calc_sun (WeatherInfo *info)
+{
+ return calc_sun_time(info, time(NULL));
+}
+
+
+/**
+ * weather_info_next_sun_event:
+ * @info: #WeatherInfo structure
+ *
+ * Returns: the interval, in seconds, until the next "sun event":
+ * - local midnight, when rise and set times are recomputed
+ * - next sunrise, when icon changes to daytime version
+ * - next sunset, when icon changes to nighttime version
+ */
+gint
+weather_info_next_sun_event (WeatherInfo *info)
+{
+ time_t now = time (NULL);
+ struct tm ltm;
+ time_t nxtEvent;
+
+ g_return_val_if_fail (info != NULL, -1);
+
+ if (!calc_sun (info))
+ return -1;
+
+ /* Determine when the next local midnight occurs */
+ (void) localtime_r (&now, &ltm);
+ ltm.tm_sec = 0;
+ ltm.tm_min = 0;
+ ltm.tm_hour = 0;
+ ltm.tm_mday++;
+ nxtEvent = mktime (&ltm);
+
+ if (info->sunsetValid &&
+ (info->sunset > now) && (info->sunset < nxtEvent))
+ nxtEvent = info->sunset;
+ if (info->sunriseValid &&
+ (info->sunrise > now) && (info->sunrise < nxtEvent))
+ nxtEvent = info->sunrise;
+ return (gint)(nxtEvent - now);
+}