Index: weather.c
===================================================================
--- weather.c	(revision 10141)
+++ weather.c	(working copy)
@@ -304,6 +304,7 @@
     gnome_vfs_async_close(handle, close_cb, info);
 
     info->sunValid = info->valid && calc_sun(info);
+    info->moonValid = info->valid && calc_moon(info);
     return;
 }
 
@@ -510,8 +511,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->radar = NULL;
     info->radar_url = prefs->radar && prefs->radar_custom_url ?
@@ -1007,6 +1010,8 @@
 static GdkPixbuf **weather_pixbufs_mini = NULL;
 static GdkPixbuf **weather_pixbufs = NULL;
 
+#define N_MOONPHASES 8
+
 enum
 {
 	PIX_UNKNOWN,
@@ -1018,9 +1023,9 @@
 	PIX_SNOW,
 	PIX_FOG,
 	PIX_MOON,
-	PIX_MOONCLOUD,
+	PIX_MOONCLOUD = PIX_MOON + N_MOONPHASES,
 
-	NUM_PIX
+	NUM_PIX = PIX_MOONCLOUD + N_MOONPHASES
 };
 
 #if 0
@@ -1053,10 +1058,25 @@
       "stock_weather-storm",
       "stock_weather-snow",
       "stock_weather-fog",
-      "stock_weather-night-clear",
-      "stock_weather-night-few-clouds"
+      "stock_weather-new-clear",
+      "stock_weather-waxing-crescent-clear",
+      "stock_weather-first-quarter-clear",
+      "stock_weather-waxing-gibbous-clear",
+      "stock_weather-full-clear",
+      "stock_weather-waning-gibbous-clear",
+      "stock_weather-third-quarter-clear",
+      "stock_weather-waning-crescent-clear",
+      "stock_weather-new-few-clouds",
+      "stock_weather-waxing-crescent-few-clouds",
+      "stock_weather-first-quarter-few-clouds",
+      "stock_weather-waxing-gibbous-few-clouds",
+      "stock_weather-full-few-clouds",
+      "stock_weather-waning-gibbous-few-clouds",
+      "stock_weather-third-quarter-few-clouds",
+      "stock_weather-waning-crescent-few-clouds",
     };
     int idx;
+    GdkPixbuf *wpmMiniFewClouds, *wpmMiniClear, *wpmFewClouds, *wpmClear;
     
     if (initialized)
        return;
@@ -1067,9 +1087,36 @@
     weather_pixbufs_mini = g_new(GdkPixbuf *, NUM_PIX);
     weather_pixbufs = g_new(GdkPixbuf *, NUM_PIX);
 
+    // Look for the backup nighttime images in case new artwork
+    // for the different moon phases is still unavailable
+    wpmMiniFewClouds = gtk_icon_theme_load_icon (icon_theme, "stock-weather-night-few-clouds",
+						 16, 0, NULL);
+    wpmMiniClear = gtk_icon_theme_load_icon (icon_theme, "stock-weather-night-clear",
+					     16, 0, NULL);
+    wpmFewClouds = gtk_icon_theme_load_icon (icon_theme, "stock-weather-night-few-clouds",
+					     48, 0, NULL);
+    wpmClear = gtk_icon_theme_load_icon (icon_theme, "stock-weather-night-clear",
+					     48, 0, NULL);
+
+    // Set up the table
     for (idx = 0; idx < NUM_PIX; idx++) {
         weather_pixbufs_mini[idx] = gtk_icon_theme_load_icon (icon_theme, icons[idx], 16, 0, NULL);
+	if (weather_pixbufs_mini[idx] == NULL) {
+	    if (idx >= PIX_MOONCLOUD) {
+	        weather_pixbufs_mini[idx] = wpmMiniFewClouds;
+	    } else if (idx >= PIX_MOON) {
+	        weather_pixbufs_mini[idx] = wpmMiniClear;
+	    }
+	}
+
         weather_pixbufs[idx] = gtk_icon_theme_load_icon (icon_theme, icons[idx], 48, 0, NULL);
+	if (weather_pixbufs[idx] == NULL) {
+	    if (idx >= PIX_MOONCLOUD) {
+	        weather_pixbufs[idx] = wpmFewClouds;
+	    } else if (idx >= PIX_MOON) {
+	        weather_pixbufs[idx] = wpmClear;
+	    }
+	}
     }
 }
 
@@ -1183,5 +1230,12 @@
         }
     }
 
+    /*
+     * If a moon image was selected above, determine which phase image to use
+     */
+    if (idx == PIX_MOON || idx == PIX_MOONCLOUD) {
+        idx += (gint)((info->moonphase * (gdouble)N_MOONPHASES / 360.) + 0.5);
+    }
+ 
     *pixbuf = pixbufs[idx];
 }
Index: weather-sun.c
===================================================================
--- weather-sun.c	(revision 10141)
+++ weather-sun.c	(working copy)
@@ -6,10 +6,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
@@ -26,35 +24,83 @@
 
 #include "weather-priv.h"
 
-#define EPOCH_TO_J2000(t)       (t-946728000)
-#define MEAN_ECLIPTIC_LONGITUDE 280.46435
-#define PERIGEE_LONGITUDE       282.94719
-#define ECCENTRICITY            0.01671002
-#define MEAN_ECLIPTIC_OBLIQUITY 
-#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,
 		     gdouble obsLat, gdouble obsLon,
@@ -72,10 +118,12 @@
     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.;
+    if ((a = fmod(a, 24.)) < 0) {
+        a += 24.;
+    }
+    if ((b = fmod(b, 24.)) < 0) {
+        b += 24.;
+    }
 
     *set = a;
     *rise = b;
@@ -98,7 +146,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;
@@ -112,45 +160,27 @@
         lcl_midn += 86400;
     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))
         return FALSE;
@@ -162,6 +192,10 @@
         set2 += 24.;
     }    
 
+    /*
+     * Calculate Greenwich sideral time at 0h UT (step 5)
+     * and at observer's longitude (steps 6,7)
+     */
     tt = t0(lcl_midn);
     t00 = tt - (gm_hoff + RADIANS_TO_HOURS(obsLon)) * 1.002737909;
 
@@ -178,17 +212,27 @@
     }
     
     /*
-     * 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);
 
     decl2 = (decl1 + decl2) / 2.;
+
+    /*
+     * Calculate an adjustment value to account for parallax,
+     * refraction and the Sun's finite diameter (steps 9,10)
+     */
     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: ChangeLog
===================================================================
--- ChangeLog	(revision 10141)
+++ ChangeLog	(working copy)
@@ -1,3 +1,18 @@
+2007-01-24  Frank Solensky  <frank@solensky.org>
+
+	reviewed by: <delete if not using a buddy>
+
+	* Makefile.am:
+	* weather-priv.h:
+	* weather-sun.c: (sunEclipLongitude), (ecl2equ), (gstObsv),
+	(calc_sun2):
+	* weather.c: (request_done), (_weather_info_fill), (init_pixbufs),
+	(_weather_info_get_pixbuf):
+	* weather-moon.c: new file
+	Include calculations for selecting different moon phase icons,
+	fall back to original moon icons if artwork is not installed
+	Closes #170277
+
 2007-01-05  Kjartan Maraas  <kmaraas@gnome.org>
 
 	* weather-met.c: s/malloc.h/stdlib.h/
Index: weather-priv.h
===================================================================
--- weather-priv.h	(revision 10141)
+++ weather-priv.h	(working copy)
@@ -107,6 +107,7 @@
 typedef gdouble WeatherPressure;
 typedef gdouble WeatherVisibility;
 typedef time_t WeatherUpdate;
+typedef gdouble WeatherMoonPhase;
 
 struct _WeatherInfo {
     WeatherForecastType forecast_type;
@@ -118,6 +119,7 @@
 
     gboolean valid;
     gboolean sunValid;
+    gboolean moonValid;
     WeatherLocation *location;
     WeatherUpdate update;
     WeatherSky sky;
@@ -130,6 +132,7 @@
     WeatherVisibility visibility;
     WeatherUpdate sunrise;
     WeatherUpdate sunset;
+    WeatherMoonPhase moonphase;
     gchar *forecast;
     gchar *metar_buffer;
     gchar *iwin_buffer;
@@ -187,10 +190,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);
@@ -205,7 +218,14 @@
 					 WeatherInfo *info);
 void		requests_done_check	(WeatherInfo *info);
 
+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);
 
 #endif /* __WEATHER_PRIV_H_ */
 
Index: Makefile.am
===================================================================
--- Makefile.am	(revision 10141)
+++ Makefile.am	(working copy)
@@ -12,7 +12,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-prefs.c gweather-prefs.h \
 	gweather-gconf.c gweather-gconf.h
 
