Help me to port this NetHack function to Python please!

Question

Hi,

I am trying to write a Python function which returns the same moon phase value as in the game NetHack. This is found in hacklib.c.

I have tried to simply copy the corresponding function from the NetHack code but I don't believe I am getting the correct results.

The function which I have written is phase_of_the_moon().

The functions position() and phase(), I found on the net, and I am using them as an indication of the success of my function. They are very accurate and give results which approximately match the nethack.alt.org server (see http://alt.org/nethack/moon/pom.txt). What I am after however is an exact replication of the original NetHack function, idiosyncrasies intact.

I would expect my function and the 'control' function to give the same moon phase at least, but currently they do not and I'm not sure why!!

Any help would be appreciated,

Thanks.

Here is the NetHack code:

/*
 * moon period = 29.53058 days ~= 30, year = 365.2422 days
 * days moon phase advances on first day of year compared to preceding year
 *  = 365.2422 - 12*29.53058 ~= 11
 * years in Metonic cycle (time until same phases fall on the same days of
 *  the month) = 18.6 ~= 19
 * moon phase on first day of year (epact) ~= (11*(year%19) + 29) % 30
 *  (29 as initial condition)
 * current phase in days = first day phase + days elapsed in year
 * 6 moons ~= 177 days
 * 177 ~= 8 reported phases * 22
 * + 11/22 for rounding
 */
int
phase_of_the_moon()  /* 0-7, with 0: new, 4: full */
{
    register struct tm *lt = getlt();
    register int epact, diy, goldn;

    diy = lt->tm_yday;
    goldn = (lt->tm_year % 19) + 1;
    epact = (11 * goldn + 18) % 30;
    if ((epact == 25 && goldn > 11) || epact == 24)
        epact++;

    return( (((((diy + epact) * 6) + 11) % 177) / 22) & 7 );
}

Here is the getlt() function (also in hacklib.c):

static struct tm *
getlt()
{
    time_t date;

#if defined(BSD) && !defined(POSIX_TYPES)
    (void) time((long *)(&date));
#else
    (void) time(&date);
#endif
#if (defined(ULTRIX) && !(defined(ULTRIX_PROTO) || defined(NHSTDC))) || (defined(BSD) && !defined(POSIX_TYPES))
    return(localtime((long *)(&date)));
#else
    return(localtime(&date));
#endif
}

Here is my Python code:

from datetime import date

def phase_of_the_moon():
   lt = date.today()

   diy = (lt - date(lt.year, 1, 1)).days
   goldn = ((lt.year - 1900) % 19) + 1
   epact = (11 * goldn + 18) % 30;
   if ((epact == 25 and goldn > 11) or epact == 24):
      epact += 1
   return ( (((((diy + epact) * 6) + 11) % 177) / 22) & 7 )

import math, decimal, datetime
dec = decimal.Decimal

def position(now=None): 
   if now is None: 
      now = datetime.datetime.now()

   diff = now - datetime.datetime(2001, 1, 1)
   days = dec(diff.days) + (dec(diff.seconds) / dec(86400))
   lunations = dec("0.20439731") + (days * dec("0.03386319269"))

   return lunations % dec(1)

def phase(pos): 
   index = (pos * dec(8)) + dec("0.5")
   index = math.floor(index)
   return {
      0: "New Moon", 
      1: "Waxing Crescent", 
      2: "First Quarter", 
      3: "Waxing Gibbous", 
      4: "Full Moon", 
      5: "Waning Gibbous", 
      6: "Last Quarter", 
      7: "Waning Crescent"
   }[int(index) & 7]

def phase2(pos): 
   return {
      0: "New Moon", 
      1: "Waxing Crescent", 
      2: "First Quarter", 
      3: "Waxing Gibbous", 
      4: "Full Moon", 
      5: "Waning Gibbous", 
      6: "Last Quarter", 
      7: "Waning Crescent"
   }[int(pos)]

def main():
   ## Correct output
   pos = position()
   phasename = phase(pos)
   roundedpos = round(float(pos), 3)
   print "%s (%s)" % (phasename, roundedpos)

   ## My output
   print "%s (%s)" % (phase2(phase_of_the_moon()), phase_of_the_moon())

if __name__=="__main__": 
   main()

Answer 1

Edit: Turns out both of the "problems" I spotted here were based on a misunderstanding of the tm struct. I'll leave the answer intact for the sake of the discussion in the comments, but save your votes for someone who might actually be correct. ;-)

---

Caveat: I'm not terribly familiar with C time constructs; I'm mostly going off the field documentation supplied for strftime.

I see two "bugs" in your port. First, I believe tm_year is intended to be the year without century, not the year minus 1900, so, goldn should be ((lt.year % 100) % 19) + 1. Secondly, your calculation for diy is zero-based, whereas tm_yday appears (again, from the docs) to be one-based. However, I'm not certain about the latter, as fixing just the goldn line gives a correct result (at least for today), where as fixing both gives the wrong answer:

>>> def phase_of_the_moon():
    lt = date.today()

    diy = (lt - date(lt.year, 1, 1)).days
    goldn = ((lt.year % 100) % 19) + 1
    epact = (11 * goldn + 18) % 30
    if ((epact == 25 and goldn > 11) or epact == 24):
        epact += 1
    return ( (((((diy + epact) * 6) + 11) % 177) / 22) & 7 )

>>> phase_of_the_moon():
3

Again, this is mostly guesswork. Please be kind. :-)

Answer 2

Curiously, when I compile and run the nethack example I get "2" as the answer ("First Quarter" which is the same as your port)

#include <time.h>

static struct tm *
getlt()
{
        time_t date;
        (void) time(&date);
        return(localtime(&date));
}
/*
 * moon period = 29.53058 days ~= 30, year = 365.2422 days
 * days moon phase advances on first day of year compared to preceding year
 *  = 365.2422 - 12*29.53058 ~= 11
 * years in Metonic cycle (time until same phases fall on the same days of
 *  the month) = 18.6 ~= 19
 * moon phase on first day of year (epact) ~= (11*(year%19) + 29) % 30
 *  (29 as initial condition)
 * current phase in days = first day phase + days elapsed in year
 * 6 moons ~= 177 days
 * 177 ~= 8 reported phases * 22
 * + 11/22 for rounding
 */
int
phase_of_the_moon()     /* 0-7, with 0: new, 4: full */
{
    register struct tm *lt = getlt();
    register int epact, diy, goldn;

    diy = lt->tm_yday;
    goldn = (lt->tm_year % 19) + 1;
    epact = (11 * goldn + 18) % 30;
    if ((epact == 25 && goldn > 11) || epact == 24)
        epact++;

    return( (((((diy + epact) * 6) + 11) % 177) / 22) & 7 );
}

int main(int argc, char * argv[]) {
    printf ("phase of the moon %d\n\n", phase_of_the_moon());
}

output:

> a.out
phase of the moon 2

But that doesn't seem like the right answer, as today, weatherunderground.com and alt.org reports the phase of the moon as "Waxing Gibbous" (a.k.a 3).

I tried remove the "-1900" but that didn't result in the right answer either.

Answer 3

The code as written is largely untestable - and you need to make it testable. So, you need the C code to be:

int
phase_of_the_moon()     /* 0-7, with 0: new, 4: full */
{
    register struct tm *lt = getlt();
    return testable_potm(lt);
}

static int
testable_potm(const struct tm *lt)
{
    register int epact, diy, goldn;

    diy = lt->tm_yday;
    goldn = (lt->tm_year % 19) + 1;
    epact = (11 * goldn + 18) % 30;
    if ((epact == 25 && goldn > 11) || epact == 24)
        epact++;

    return( (((((diy + epact) * 6) + 11) % 177) / 22) & 7 );
}

Now you can run tests with multiple values of time. The alternative way to do this is to fake getlt() instead.

You then need parallel changes in your Python code. Then you create a file of time_t values which can be read by both Python and C, and then converted into an appropriate structure (via localtime() in C). Then you can see where things are deviating.

Answer 4

I like to think I know a thing or two about calendars, so let's see if I can clear a few things up.

The Catholic Church defines the date of Easter in terms of lunar phases (this is why the date jumps around from year to year). Because of this, it needs to be able to calculate the approximate moon phase, and its algorithm for doing so is explained here.

I have not done very detailed checking, but it appears that the NetHack algorithm is based heavily on the Church's algorithm. The NetHack algorithm seems to, like the Church's algorithm, pay attention only to the calendar date, ignoring time zones and the time of day.

The NetHack algorithm uses only the year and the day of the year. I can tell from examining the code that, to be Y2K compatible, that tm_year has to be the year minus 1900.

Answer 5

Following code is borrowed from this site, pasting it here for easy reference (and in case the other site goes down). Seems to do what you want.

# Determine the moon phase of a date given
# Python code by HAB

def moon_phase(month, day, year):
    ages = [18, 0, 11, 22, 3, 14, 25, 6, 17, 28, 9, 20, 1, 12, 23, 4, 15, 26, 7]
    offsets = [-1, 1, 0, 1, 2, 3, 4, 5, 7, 7, 9, 9]
    description = ["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)"]
    months = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]

    if day == 31:
        day = 1
    days_into_phase = ((ages[(year + 1) % 19] + ((day + offsets[month-1]) % 30) + (year < 1900)) % 30)
    index = int((days_into_phase + 2) * 16/59.0)
    if index > 7:
        index = 7
    status = description[index]

    # light should be 100% 15 days into phase
    light = int(2 * days_into_phase * 100/29)
    if light > 100:
        light = abs(light - 200);
    date = "%d%s%d" % (day, months[month-1], year)

    return date, status, light

# put in a date you want ...
month = 5
day = 14
year = 2006  # use yyyy format

date, status, light = moon_phase(month, day, year)
print "moon phase on %s is %s, light = %d%s" % (date, status, light, '%')

You can use the time module to get the current local time. Heres how I did it (paste below posted code to testrun):

import time
tm = time.localtime()
month = tm.tm_mon
day = tm.tm_mday
year = tm.tm_year
date, status, light = moon_phase(month, day, year)
print "moon phase on %s is %s, light = %d%s" % (date, status, light, '%')

Output:

moon phase on 22Dec2009 is waxing crescent (increasing to full), light = 34%

Moon stuff is fun. :)