Index: ChangeLog
===================================================================
--- ChangeLog	(revision 612)
+++ ChangeLog	(working copy)
@@ -1,3 +1,26 @@
+2009-02-14: Frank Solensky  <frank@solensky.org>
+
+	reviewed by: <delete if not using a buddy>
+
+	Bug #564598: add phase-of-moon calculations to library
+
+	* libgweather/weather-moon.c: new file
+	* libgweather/Makefile.am: build new source file
+	* libgweather/weather.h: add new typedefs and function.
+	* libgweather/weather-priv.h: add calc_moon(); move some macros
+	  out of weather-sun.c; calculation corrections
+	* libgweather/weather-sun.c (sunEclipLongitude) new function,
+	  (ecl2equ) change from static to not, gives calc_moon access
+	  (gstObsv) push section out to sunEclipLongitude;
+          calculation corrections
+	* libgweather/weather.c (request_done) add call to calc_moon(),
+          (_weather_info_fill) initialize new values,
+	  (weather_info_get_icon_name) include moon phase references,
+          (weather_info_get_value_update), (weather_info_get_value_sunrise),
+	    (weather_info_get_value_sunset) zero is an OK value to return,
+	  (weather_info_get_value_moonphase),
+            (weather_info_get_upcoming_moonphases): new functions
+
 2009-01-19  Vincent Untz  <vuntz@gnome.org>
 
 	* configure.in: post-release bump to 2.25.90
Index: libgweather/weather.c
===================================================================
--- libgweather/weather.c	(revision 612)
+++ libgweather/weather.c	(working copy)
@@ -45,6 +45,8 @@
 #include "gweather-win32.h"
 #endif
 
+#define MOON_PHASES 36
+
 static void _weather_internal_check (void);
 
 const char *
@@ -336,8 +338,10 @@
 
 void request_done (WeatherInfo *info, gboolean ok)
 {
-    if (ok)
+    if (ok) {
 	info->sunValid = info->valid && calc_sun (info);
+	info->moonValid = info->valid && calc_moon (info);
+    }
 
     if (!--info->requests_pending)
         info->finish_cb (info, info->cb_data);
@@ -520,8 +524,11 @@
     info->pressure = -1.0;
     info->visibility = -1.0;
     info->sunValid = FALSE;
+    info->moonValid = FALSE;
     info->sunrise = 0;
     info->sunset = 0;
+    info->moonphase = 0;
+    info->moonlatitude = 0;
     info->forecast = NULL;
     info->forecast_list = NULL;
     info->radar = NULL;
@@ -1117,9 +1124,14 @@
 weather_info_get_icon_name (WeatherInfo *info)
 {
     WeatherConditions cond;
-    WeatherSky sky;
-    time_t current_time;
-    gboolean daytime;
+    WeatherSky   sky;
+    time_t       current_time;
+    gboolean     daytime;
+    gchar*       icon;
+    static gchar icon_buffer[32];
+    WeatherMoonPhase moonPhase;
+    WeatherMoonLatitude moonLat;
+    gint             phase;
 
     g_return_val_if_fail (info != NULL, NULL);
 
@@ -1185,7 +1197,8 @@
 	if (daytime)
 	    return "weather-clear";
 	else
-	    return "weather-clear-night";
+	    icon = g_stpcpy(icon_buffer, "weather-clear-night");
+	break;
 
     case SKY_BROKEN:
     case SKY_SCATTERED:
@@ -1193,13 +1206,38 @@
 	if (daytime)
 	    return "weather-few-clouds";
 	else
-	    return "weather-few-clouds-night";
+	    icon = g_stpcpy(icon_buffer, "weather-few-clouds-night");
+	break;
 
     case SKY_OVERCAST:
 	return "weather-overcast";
+
+    default:  /* unrecognized */
+	return NULL;
     }
 
-    return NULL;
+    /*
+     * A phase-of-moon icon is to be returned.
+     * Determine which one based on the moon's location
+     */
+    if (info->moonValid && weather_info_get_value_moonphase(info, &moonPhase, &moonLat)) {
+	
+	phase = (gint)((moonPhase * MOON_PHASES / 360.) + 0.5);
+	if (phase == MOON_PHASES) {
+	    phase = 0;
+	} else if (phase > 0 &&
+		   weather_info_get_location(info)->latitude < moonLat) {
+	    /*
+	     * Locations south of the moon's latitude will see the moon in the northern sky.
+	     * The moon waxes and wanes from left to right so we reference an icon running
+	     * in the opposite direction.
+	     */
+	    phase = MOON_PHASES - phase;
+	}
+	g_snprintf(icon, sizeof(icon_buffer) - strlen(icon_buffer), "_%03d",
+		   phase * 360 / MOON_PHASES);
+    }
+    return icon_buffer;
 }
 
 static gboolean
@@ -1461,7 +1499,7 @@
 
     *value = info->update;
 
-    return *value != 0;
+    return TRUE;
 }
 
 gboolean
@@ -1475,7 +1513,7 @@
 
     *value = info->sunrise;
 
-    return *value != 0;
+    return TRUE;
 }
 
 gboolean
@@ -1489,10 +1527,27 @@
 
     *value = info->sunset;
 
-    return *value != 0;
+    return TRUE;
 }
 
 gboolean
+weather_info_get_value_moonphase (WeatherInfo *info,
+				  WeatherMoonPhase *value,
+				  WeatherMoonLatitude *lat)
+{
+    g_return_val_if_fail (info != NULL, FALSE);
+    g_return_val_if_fail (value != NULL, FALSE);
+
+    if (!info->valid || !info->moonValid)
+	return FALSE;
+
+    *value = info->moonphase;
+    *lat = info->moonlatitude;
+
+    return TRUE;
+}
+
+gboolean
 weather_info_get_value_wind (WeatherInfo *info, SpeedUnit unit, gdouble *speed, WeatherWindDirection *direction)
 {
     gboolean res = FALSE;
@@ -1537,6 +1592,24 @@
     return distance_value (info->visibility, unit, value, info->distance_unit);
 }
 
+/**
+ * weather_info_get_upcoming_moonphases:
+ * @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
+weather_info_get_upcoming_moonphases (WeatherInfo *info, time_t *phases)
+{
+    g_return_val_if_fail (info != NULL, FALSE);
+    g_return_val_if_fail (phases != NULL, FALSE);
+
+    return calc_moon_phases(info, phases);
+}
+
 static void
 _weather_internal_check (void)
 {
Index: libgweather/weather-sun.c
===================================================================
--- libgweather/weather-sun.c	(revision 612)
+++ libgweather/weather-sun.c	(working copy)
@@ -20,10 +20,8 @@
  * Formulas from:
  * "Practical Astronomy With Your Calculator" (3e), Peter Duffett-Smith
  * Cambridge University Press 1988
- *
- * Planetary Mean Orbit parameters from http://ssd.jpl.nasa.gov/elem_planets.html,
- * converting longitudes from heliocentric to geocentric coordinates
- * epoch J2000 (2000 Jan 1 12:00:00 UTC)
+ * Unless otherwise noted, comments referencing "steps" are related to
+ * the algorithm presented in section 49 of above
  */
 
 #ifdef HAVE_CONFIG_H
@@ -41,19 +39,56 @@
 #include "gweather-win32.h"
 #endif
 
-#define EPOCH_TO_J2000(t)       (t - 946728000)
-#define MEAN_ECLIPTIC_LONGITUDE 280.46435
-#define PERIGEE_LONGITUDE       282.94719
-#define ECCENTRICITY            0.01671002
-#define SOL_PROGRESSION         (360./365.242191)
+#define ECCENTRICITY(d)         (0.01671123 - (d)/36525.*0.00004392)
 
+/* 
+ * Ecliptic longitude of the sun at specified time (UT)
+ * The algorithm 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;
+}
+
 /*
- * Convert ecliptic longitude (radians) to equitorial coordinates,
- * expressed as right ascension (hours) and declination (radians)
- * Assume ecliptic latitude is 0.0
+ * Convert ecliptic longitude and latitude (radians) to equitorial
+ * coordinates, expressed as right ascension (hours) and
+ * declination (radians)
  */
-static void
-ecl2equ (gdouble time, gdouble eclipLon,
+void
+ecl2equ (gdouble time,
+	 gdouble eclipLon, gdouble eclipLat,
 	 gdouble *ra, gdouble *decl)
 {
     gdouble jc = EPOCH_TO_J2000 (time) / (36525. * 86400.);
@@ -61,16 +96,23 @@
 					      - .013004166 * jc
 					      - 1.666667e-7 * jc * jc
 					      + 5.027777e-7 * jc * jc * jc);
-    *ra = RADIANS_TO_HOURS (atan2 (sin (eclipLon) * cos (mEclipObliq),
-				   cos (eclipLon)));
-    if (*ra < 0.)
-        *ra += 24.;
-    *decl = asin (sin (mEclipObliq) * sin (eclipLon));
+    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 of an object
- * based on it equitorial coordinates
+ * 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,
@@ -116,7 +158,7 @@
 {
     time_t gm_midn;
     time_t lcl_midn;
-    gdouble gm_hoff, ndays, meanAnom, eccenAnom, delta, lambda;
+    gdouble gm_hoff, lambda;
     gdouble ra1, ra2, decl1, decl2;
     gdouble rise1, rise2, set1, set2;
     gdouble tt, t00;
@@ -131,39 +173,21 @@
     else if (lcl_midn > t)
         lcl_midn -= 86400;
 
+    lambda = sunEclipLongitude (lcl_midn);
+    
     /*
-     * Ecliptic longitude of the sun at midnight local time
-     * Start with an estimate based on a fixed daily rate
+     * Calculate equitorial coordinates of sun at previous
+     * and next local midnights
      */
-    ndays = EPOCH_TO_J2000 (lcl_midn) / 86400.;
-    meanAnom = DEGREES_TO_RADIANS (ndays * SOL_PROGRESSION
-				   + MEAN_ECLIPTIC_LONGITUDE - PERIGEE_LONGITUDE);
-
+    ecl2equ (lcl_midn, lambda, 0., &ra1, &decl1);
+    ecl2equ (lcl_midn + 86400.,
+	     lambda + DEGREES_TO_RADIANS(SOL_PROGRESSION), 0.,
+	     &ra2, &decl2);
+    
     /*
-     * Approximate solution of Kepler's equation:
-     * Find E which satisfies  E - e sin(E) = M (mean anomaly)
+     * Convert to rise and set times based positions for the preceding
+     * and following local midnights.
      */
-    eccenAnom = meanAnom;
-
-    while (1e-12 < fabs (delta = eccenAnom - ECCENTRICITY * sin (eccenAnom) - meanAnom))
-	eccenAnom -= delta / (1. - ECCENTRICITY * cos (eccenAnom));
-
-    /* Sun's longitude on the ecliptic */
-    lambda = fmod (DEGREES_TO_RADIANS (PERIGEE_LONGITUDE)
-		   + 2. * atan (sqrt ((1. + ECCENTRICITY) / (1. - ECCENTRICITY))
-				* tan (eccenAnom / 2.)),
-		   2. * M_PI);
-
-    /* Calculate equitorial coordinates of sun at previous and next local midnights */
-
-    ecl2equ (lcl_midn, lambda, &ra1, &decl1);
-    ecl2equ (lcl_midn + 86400., lambda + DEGREES_TO_RADIANS (SOL_PROGRESSION), &ra2, &decl2);
-
-    /* Convert to rise and set times assuming the earth were to stay in its position
-     * in its orbit around the sun.  This will soon be followed by interpolating
-     * between the two
-     */
-
     gstObsv (ra1, decl1, obsLat, obsLon - (gm_hoff * M_PI / 12.), &rise1, &set1);
     gstObsv (ra2, decl2, obsLat, obsLon - (gm_hoff * M_PI / 12.), &rise2, &set2);
 
@@ -192,16 +216,25 @@
     }
 
     /*
-     * Interpolate between the two
+     * 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;
     set1  = (set1 + dt - tt) * 0.9972695661;
     if (rise1 < 0.)
Index: libgweather/weather.h
===================================================================
--- libgweather/weather.h	(revision 612)
+++ libgweather/weather.h	(working copy)
@@ -259,6 +259,8 @@
 };
 
 typedef enum _WeatherConditionQualifier WeatherConditionQualifier;
+typedef gdouble WeatherMoonPhase;
+typedef gdouble WeatherMoonLatitude;
 
 gboolean weather_info_get_value_update		(WeatherInfo *info, time_t *value);
 gboolean weather_info_get_value_sky		(WeatherInfo *info, WeatherSky *sky);
@@ -273,6 +275,8 @@
 gboolean weather_info_get_value_visibility	(WeatherInfo *info, DistanceUnit unit, gdouble *value);
 gboolean weather_info_get_value_sunrise		(WeatherInfo *info, time_t *value);
 gboolean weather_info_get_value_sunset 		(WeatherInfo *info, time_t *value);
+gboolean weather_info_get_value_moonphase       (WeatherInfo *info, WeatherMoonPhase *value, WeatherMoonLatitude *lat);
+gboolean weather_info_get_upcoming_moonphases   (WeatherInfo *info, time_t *phases);
 
 G_END_DECLS
 
Index: libgweather/weather-priv.h
===================================================================
--- libgweather/weather-priv.h	(revision 612)
+++ libgweather/weather-priv.h	(working copy)
@@ -73,6 +73,7 @@
     gboolean valid;
     gboolean network_error;
     gboolean sunValid;
+    gboolean moonValid;
     gboolean tempMinMaxValid;
     WeatherLocation *location;
     WeatherUpdate update;
@@ -88,6 +89,8 @@
     WeatherVisibility visibility;
     WeatherUpdate sunrise;
     WeatherUpdate sunset;
+    WeatherMoonPhase moonphase;
+    WeatherMoonLatitude moonlatitude;
     gchar *forecast;
     GSList *forecast_list; /* list of WeatherInfo* for the forecast, NULL if not available */
     gchar *radar_buffer;
@@ -139,10 +142,19 @@
 #define VISIBILITY_SM_TO_KM(sm)		((sm) * 1.609344)
 #define VISIBILITY_SM_TO_M(sm)		(VISIBILITY_SM_TO_KM (sm) * 1000)
 
-#define DEGREES_TO_RADIANS(deg)		((fmod (deg,360.) / 180.) * M_PI)
+#define DEGREES_TO_RADIANS(deg)		((fmod ((deg),360.) / 180.) * M_PI)
 #define RADIANS_TO_DEGREES(rad)		((rad) * 180. / M_PI)
 #define RADIANS_TO_HOURS(rad)		((rad) * 12. / M_PI)
 
+/*
+ * Planetary Mean Orbit and their progressions from J2000 are based on the
+ * values in http://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf
+ */
+#define EPOCH_TO_J2000(t)          ((gdouble)(t)-946727935.816)
+#define MEAN_ECLIPTIC_LONGITUDE(d) (280.46457166 + (d)/36525.*35999.37244981)
+#define SOL_PROGRESSION            (360./365.242191)
+#define PERIGEE_LONGITUDE(d)       (282.93768193 + (d)/36525.*0.32327364)
+
 void		metar_start_open	(WeatherInfo *info);
 void		iwin_start_open		(WeatherInfo *info);
 void		metoffice_start_open	(WeatherInfo *info);
@@ -156,9 +168,16 @@
 void		request_done		(WeatherInfo *info,
 					 gboolean     ok);
 
+void		ecl2equ			(gdouble t,
+					 gdouble eclipLon,
+					 gdouble eclipLat,
+					 gdouble *ra,
+					 gdouble *decl);
+gdouble 	sunEclipLongitude	(time_t t);
 gboolean	calc_sun		(WeatherInfo *info);
+gboolean	calc_moon		(WeatherInfo *info);
+gboolean	calc_moon_phases	(WeatherInfo *info, time_t *phases);
 
 void		free_forecast_list	(WeatherInfo *info);
 
 #endif /* __WEATHER_PRIV_H_ */
-
Index: libgweather/Makefile.am
===================================================================
--- libgweather/Makefile.am	(revision 612)
+++ libgweather/Makefile.am	(working copy)
@@ -18,7 +18,7 @@
 	weather.c weather.h weather-priv.h \
 	weather-metar.c weather-iwin.c weather-met.c \
 	weather-bom.c weather-wx.c \
-	weather-sun.c \
+	weather-sun.c weather-moon.c \
 	gweather-enum-types.c \
 	gweather-prefs.c gweather-prefs.h \
 	gweather-gconf.c gweather-gconf.h \
