diff --git a/bin/weeutil/Moon.py b/bin/weeutil/Moon.py
index b4ed80cb..b876e8f9 100644
--- a/bin/weeutil/Moon.py
+++ b/bin/weeutil/Moon.py
@@ -1,42 +1,46 @@
 #
-#    Copyright (c) 2009-2015 Tom Keffer <tkeffer@gmail.com>
+#    Copyright (c) 2009-2018 Tom Keffer <tkeffer@gmail.com>
 #
 #    See the file LICENSE.txt for your full rights.
 #
-"""Determine the phase of the moon phase given a date.
-    CF: http://en.wikipedia.org/wiki/Lunar_phase
-"""
+"""Determine the phase of the moon phase given a date."""
 
-import datetime
+import time
 import math
-import decimal
 
-dec = decimal.Decimal
-
-moon_phases = ["new (totally dark)", 
-               "waxing crescent (increasing to full)", 
-               "in its first quarter (increasing to full)", 
-               "waxing gibbous (increasing to full)", 
-               "full (full light)", 
-               "waning gibbous (decreasing from full)", 
-               "in its last quarter (decreasing from full)", 
+moon_phases = ["new (totally dark)",
+               "waxing crescent (increasing to full)",
+               "in its first quarter (increasing to full)",
+               "waxing gibbous (increasing to full)",
+               "full (full light)",
+               "waning gibbous (decreasing from full)",
+               "in its last quarter (decreasing from full)",
                "waning crescent (decreasing from full)"]
 
+# First new moon of 2018: 16-Jan-2018 at 1017UTC
+new_moon_2018 = 1516097820
 
-def moon_phase(year, month, day):
+
+def moon_phase(year, month, day, hour=12):
     """Calculates the phase of the moon, given a year, month, day.
-    
+
     returns: a tuple. First value is an index into an array
     of moon phases, such as Moon.moon_phases above. Second
     value is the percent fullness of the moon.
     """
-    time_dt = datetime.datetime(year, month, day)
-    diff = time_dt - datetime.datetime(2001, 1, 1)
-    
-    days = dec(diff.days) + (dec(diff.seconds) / dec(86400))
-    lunations = dec('0.20439731') + (days / dec('29.530589'))
+
+    # Convert to UTC
+    time_ts = time.mktime((year, month, day, hour, 0, 0, 0, 0, -1))
+
+    # How many days since the first moon of 2018
+    delta_days = (time_ts - new_moon_2018) / 86400
+    # Number of lunations
+    lunations = delta_days / 29.530589
+
+    # The fraction of the lunar cycle
     position = float(lunations) % 1.0
-    fullness = int(100.0*(1.0 - math.cos(2.0 * math.pi * position))/2.0 + 0.5)
+    # The percent illumination, rounded to the nearest integer
+    fullness = int(100.0 * (1.0 - math.cos(2.0 * math.pi * position)) / 2.0 + 0.5)
     index = int((position * 8) + 0.5) & 7
 
-    return (index, fullness)
+    return index, fullness
diff --git a/docs/changes.txt b/docs/changes.txt
index f4e5dab8..17d181cb 100644
--- a/docs/changes.txt
+++ b/docs/changes.txt
@@ -34,6 +34,10 @@ Guard against a negative value for 'interval' in WMR200 stations.
 Fixed bug that caused older, "type A" Vantage Pro2, to crash.
 Fixes issue #343.
 
+Changed the formula used to calculate percentage illumination of the moon to
+something more accurate around 2018. This formula is only used if pyephem is not
+installed.
+
 
 3.8.2 08/15/2018