Does anyone know an algorithm to either calculate the moon phase or age on a given date or find the dates for new/full moons in a given year?

Googling tells me the answer is in some Astronomy book, but I don't really want to buy a whole book when I only need a single page.

Update:

I should have qualified my statement about googling a little better. I did find solutions that only worked over some subset of time (like the 1900's); and the trig based solutions that would be more computationally expensive than I'd like.

S Lott in his Python book has several algorithms for calculating Easter on a given year, most are less than ten lines of code and some work for all days in the Gregorian calendar. Finding the full moon in March is a key piece of finding Easter so I figured there should be an algorithm that doesn't require trig and works for all dates in the Gregorian calendar.

If you're like me, you try to be a careful programmer. So it makes you nervous when you see random code scattered across the internet that purports to solve a complex astronomical problem, but doesn't explain why the solution is correct.

You believe that there must be authoritative sources such as books which contain careful, and complete, solutions. For instance:

Meeus, Jean. Astronomical Algorithms. Richmond: Willmann-Bell, 1991. ISBN 0-943396-35-2.

Duffett-Smith, Peter. Practical Astronomy With Your Calculator. 3rd ed. Cambridge: Cambridge University Press, 1981. ISBN 0-521-28411-2.

You place your trust in widely-used, well-tested, open source libraries which can have their errors corrected (unlike static web pages). Here then, is a Python solution to your question based on the PyEphem library, using the Phases of the Moon interface.

#!/usr/bin/python
import datetime
import ephem
from typing import List, Tupledef get_phase_on_day(year: int, month: int, day: int):"""Returns a floating-point number from 0-1. where 0=new, 0.5=full, 1=new"""#Ephem stores its date numbers as floating points, which the following uses#to conveniently extract the percent time between one new moon and the next#This corresponds (somewhat roughly) to the phase of the moon.#Use Year, Month, Day as argumentsdate = ephem.Date(datetime.date(year,month,day))nnm = ephem.next_new_moon(date)pnm = ephem.previous_new_moon(date)lunation = (date-pnm)/(nnm-pnm)#Note that there is a ephem.Moon().phase() command, but this returns the#percentage of the moon which is illuminated. This is not really what we want.return lunationdef get_moons_in_year(year: int) -> List[Tuple[ephem.Date, str]]:"""Returns a list of the full and new moons in a year. The list contains tuples
of either the form (DATE,'full') or the form (DATE,'new')"""moons=[]date=ephem.Date(datetime.date(year,1,1))while date.datetime().year==year:date=ephem.next_full_moon(date)moons.append( (date,'full') )date=ephem.Date(datetime.date(year,1,1))while date.datetime().year==year:date=ephem.next_new_moon(date)moons.append( (date,'new') )#Note that previous_first_quarter_moon() and previous_last_quarter_moon()#are also methodsmoons.sort(key=lambda x: x[0])return moonsprint(get_phase_on_day(2013,1,1))print(get_moons_in_year(2013))

This returns

0.632652265318[(2013/1/11 19:43:37, 'new'), (2013/1/27 04:38:22, 'full'), (2013/2/10 07:20:06, 'new'), (2013/2/25 20:26:03, 'full'), (2013/3/11 19:51:00, 'new'), (2013/3/27 09:27:18, 'full'), (2013/4/10 09:35:17, 'new'), (2013/4/25 19:57:06, 'full'), (2013/5/10 00:28:22, 'new'), (2013/5/25 04:24:55, 'full'), (2013/6/8 15:56:19, 'new'), (2013/6/23 11:32:15, 'full'), (2013/7/8 07:14:16, 'new'), (2013/7/22 18:15:31, 'full'), (2013/8/6 21:50:40, 'new'), (2013/8/21 01:44:35, 'full'), (2013/9/5 11:36:07, 'new'), (2013/9/19 11:12:49, 'full'), (2013/10/5 00:34:31, 'new'), (2013/10/18 23:37:39, 'full'), (2013/11/3 12:49:57, 'new'), (2013/11/17 15:15:44, 'full'), (2013/12/3 00:22:22, 'new'), (2013/12/17 09:28:05, 'full'), (2014/1/1 11:14:10, 'new'), (2014/1/16 04:52:10, 'full')]