Index: ChangeLog
===================================================================
--- ChangeLog	(revision 553)
+++ ChangeLog	(working copy)
@@ -1,3 +1,17 @@
+2008-12-24  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:
+	* libgweather/weather-sun.c (sunEclipLongitude), (ecl2equ),
+	(calc_sun2):
+	* libgweather/weather-priv.h:
+	* libgweather/weather.c (request_done), (_weather_info_fill),
+	(weather_info_get_icon_name), (weather_info_get_value_moonphase):
+	* libgweather/weather.h:
+
 2008-12-19  Milan Crha  <mcrha@redhat.com>
 
 	** Part of fix for bug #564652
Index: libgweather/weather.c
===================================================================
--- libgweather/weather.c	(revision 553)
+++ libgweather/weather.c	(working copy)
@@ -336,9 +336,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 +521,10 @@
     info->pressure = -1.0;
     info->visibility = -1.0;
     info->sunValid = FALSE;
+    info->moonValid = FALSE;
     info->sunrise = 0;
     info->sunset = 0;
+    info->moonphase = 0;
     info->forecast = NULL;
     info->forecast_list = NULL;
     info->radar = NULL;
@@ -1114,13 +1117,16 @@
 }
 
 const gchar *
-weather_info_get_icon_name (WeatherInfo *info)
+weather_info_get_icon_name (WeatherInfo *info, gint moon_phases)
 {
     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];
+    gdouble      moonphase;
+    gint         phase;
     g_return_val_if_fail (info != NULL, NULL);
 
     if (!info->valid)
@@ -1184,22 +1190,51 @@
     case SKY_CLEAR:
 	if (daytime)
 	    return "weather-clear";
+	else if (moon_phases <= 0)
+	    return "weather-clear-night";
 	else
-	    return "weather-clear-night";
+	    icon = g_stpcpy(icon_buffer, "weather-clear-night");
+	break;
 
     case SKY_BROKEN:
     case SKY_SCATTERED:
     case SKY_FEW:
 	if (daytime)
 	    return "weather-few-clouds";
+	else if (moon_phases <= 0)
+	    return "weather-few-clouds-night";
 	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 caller's indication of how many are available
+     */
+    if (weather_info_get_value_moonphase(info, &moonphase)) {
+	phase = (gint)((moonphase * moon_phases / 360.) + 0.5);
+	if (phase == moon_phases) {
+	    phase = 0;
+	} else if (phase > 0 &&
+		   weather_info_get_location(info)->latitude < 0) {
+	    /*
+	     * Locations in the southern hemisphere see the moon in the northern sky.
+	     * The moon phase advances to the right.
+	     * (consider flipping when south of moon's current declination rather than equator)
+	     */
+	    phase = moon_phases - phase;
+	}
+	g_snprintf(icon, sizeof(icon_buffer) - strlen(icon_buffer), "_%03d",
+		   phase * 360 / moon_phases);
+    }
+    return icon_buffer;
 }
 
 static gboolean
@@ -1493,6 +1528,20 @@
 }
 
 gboolean
+weather_info_get_value_moonphase (WeatherInfo *info, gdouble *value)
+{
+    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;
+
+    return *value != 0;
+}
+
+gboolean
 weather_info_get_value_wind (WeatherInfo *info, SpeedUnit unit, gdouble *speed, WeatherWindDirection *direction)
 {
     gboolean res = FALSE;
Index: libgweather/weather-sun.c
===================================================================
--- libgweather/weather-sun.c	(revision 553)
+++ 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,36 +39,81 @@
 #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 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(ndays * SOL_PROGRESSION
+				  + 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));
+    }
+
+    /*
+     * Sun'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
  */
-static void
-ecl2equ (gdouble time, gdouble eclipLon,
-	 gdouble *ra, gdouble *decl)
+void ecl2equ (gdouble time,
+	      gdouble eclipLon, gdouble eclipLat,
+	      gdouble *ra, gdouble *decl)
 {
     gdouble jc = EPOCH_TO_J2000 (time) / (36525. * 86400.);
     gdouble mEclipObliq = DEGREES_TO_RADIANS (23.439291667
 					      - .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 +159,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,56 +174,42 @@
     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;
+    gstObsv (ra1, decl1, obsLat, obsLon - (gm_hoff * M_PI / 12.),
+	     &rise1, &set1);
+    gstObsv (ra2, decl2, obsLat, obsLon - (gm_hoff * M_PI / 12.),
+	     &rise2, &set2);
 
-    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);
-
     /* TODO: include calculations for regions near the poles. */
-    if (isnan (rise1) || isnan (rise2))
+    if (isnan(rise1) || isnan(rise2))
         return FALSE;
 
-    if (rise2 < rise1)
+    if (rise2 < rise1) {
         rise2 += 24.;
-    if (set2 < set1)
+    }
+    if (set2 < set1) {
         set2 += 24.;
+    }
 
-    tt = t0 (lcl_midn);
-    t00 = tt - (gm_hoff + RADIANS_TO_HOURS (obsLon)) * 1.002737909;
+    tt = t0(lcl_midn);
+    t00 = tt - (gm_hoff + RADIANS_TO_HOURS(obsLon)) * 1.002737909;
 
     if (t00 < 0.)
-	t00 += 24.;
+        t00 += 24.;
 
     if (rise1 < t00) {
         rise1 += 24.;
@@ -190,18 +219,27 @@
         set1  += 24.;
         set2  += 24.;
     }
-
+   
     /*
-     * 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));
-
+    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 553)
+++ libgweather/weather.h	(working copy)
@@ -171,7 +171,7 @@
 const gchar *		weather_info_get_temp_summary	(WeatherInfo *info);
 gchar *			weather_info_get_weather_summary(WeatherInfo *info);
 
-const gchar *		weather_info_get_icon_name	(WeatherInfo *info);
+const gchar *		weather_info_get_icon_name	(WeatherInfo *info, gint moon_phases);
 gint			weather_info_next_sun_event	(WeatherInfo *info);
 
 /* values retrieving functions */
@@ -273,6 +273,7 @@
 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, gdouble* moon_phases);
 
 G_END_DECLS
 
Index: libgweather/weather-priv.h
===================================================================
--- libgweather/weather-priv.h	(revision 553)
+++ libgweather/weather-priv.h	(working copy)
@@ -61,6 +61,7 @@
 typedef gdouble WeatherPressure;
 typedef gdouble WeatherVisibility;
 typedef time_t WeatherUpdate;
+typedef gdouble WeatherMoonPhase;
 
 struct _WeatherInfo {
     WeatherForecastType forecast_type;
@@ -73,6 +74,7 @@
     gboolean valid;
     gboolean network_error;
     gboolean sunValid;
+    gboolean moonValid;
     gboolean tempMinMaxValid;
     WeatherLocation *location;
     WeatherUpdate update;
@@ -88,6 +90,7 @@
     WeatherVisibility visibility;
     WeatherUpdate sunrise;
     WeatherUpdate sunset;
+    WeatherMoonPhase moonphase;
     gchar *forecast;
     GSList *forecast_list; /* list of WeatherInfo* for the forecast, NULL if not available */
     gchar *radar_buffer;
@@ -139,10 +142,20 @@
 #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
+ * converting longitudes from heliocentric to geocentric coordinates (+180)
+ */
+#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,7 +169,14 @@
 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);
 
 void		free_forecast_list	(WeatherInfo *info);
 
Index: libgweather/Makefile.am
===================================================================
--- libgweather/Makefile.am	(revision 553)
+++ 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 \
