Phases Of The Moon

January 22, 2010

We calculate the moon’s age in days since the last new moon using the formula given previously, using the julian function from the Standard Prelude to number the days. Normalize is a convenience function that extracts the fractional part of a real number:

(define (days year month day)
  (define (normalize x)
    (let ((x (- x (floor x))))
      (if (negative? x) (+ x 1) x)))
  (let ((new 2451550.1) (moon 29.530588853)
        (j (julian year month day)))
    (* (normalize (/ (- j new) moon)) moon)))

We could initialize a vector and calculate the phase as an index into the vector, but the calculation involves real numbers, so we choose a safer approach:

(define (phase year month day)
  (let ((d (days year month day)))
    (cond ((< d  1.84566) "New")
          ((< d  5.53699) "Waxing crescent")
          ((< d  9.22831) "First quarter")
          ((< d 12.91963) "Waxing gibbous")
          ((< d 16.61096) "Full")
          ((< d 20.30228) "Waning gibbous")
          ((< d 23.99361) "Last quarter")
          ((< d 27.68493) "Waning crescent")
          (else           "New"))))

Here are some examples:

> (phase 2000 1 6)
"New"
> (days 2010 1 22)
7.106982227906547
> (phase 2010 1 22)
"First quarter"

You can run the program at http://programmingpraxis.codepad.org/23rZcxsl.

About these ads

Pages: 1 2

8 Responses to “Phases Of The Moon”

  1. Remco Niemeijer said

    My Haskell solution (see http://bonsaicode.wordpress.com/2010/01/22/programming-praxis-phases-of-the-moon/ for a version with comments):

    import Data.Time
    import Data.Time.Calendar.Julian
    import Data.Fixed
    
    moonPhase :: Integer -> Int -> Int -> String
    moonPhase y m d = phase . flip mod' 29.530588853 . fromIntegral $
        diffDays (fromJulian y m d) (fromJulian 2000 1 6) where
        phase x | x <  1.84566 = "New"
                | x <  5.53699 = "Waxing crescent"
                | x <  9.22831 = "First quarter"
                | x < 12.91963 = "Waxing gibbous"
                | x < 16.61096 = "Full"
                | x < 20.30228 = "Waning gibbous"
                | x < 23.99361 = "Last quarter"
                | x < 27.68493 = "Waning crescent"
                | otherwise    = "New"
    
  2. novatech said

    ruby version

    def julian (year, month, day)
    a = (14-month)/12
    y = year+4800-a
    m = (12*a)-3+month
    return day + (153*m+2)/5 + (365*y) + y/4 - y/100 + y/400 - 32045
    end
    def phase (year,month,day)
    p=(julian(year,month,day)-julian(2000,1,6))%29.530588853
    if p<1.84566: return "New"
    elsif p<5.53699: return "Waxing crescent"
    elsif p<9.22831: return "First quarter"
    elsif p<12.91963: return "Waxing gibbous"
    elsif p<16.61096: return "Full"
    elsif p<20.30228: return "Waning gibbous"
    elsif p<23.99361: return "Last quarter"
    elsif p<27.68493: return "Waning crescent"
    else return "New"
    end
    end
    
    print "#{phase(2020,1,23)}\n"
    print "#{phase(1999,1,6)}\n"
    print "#{phase(2010,2,10)}\n"
    print "#{phase(1987,5,10)}\n"
    
    
  3. iyo said

    There’s my python version. I’m surprised python’s math.fmod(num, n) doesn’t give results into <0;n) but it is simple remainder – results are from interval (-n;n). Inverse indexing of list is also possible.

    Listing of the code:


    from datetime import date
    from math import fmod

    def phaseOfMoon(day):
    period = 29.530588853
    referenceDate = date(2000, 1, 6)
    phases = ["new", "waxing crescent", "first quarter", "waxing gibbous",
    "full", "waning gibbous", "last quarter", "waning crescent"]

    daysDelta = (day - referenceDate).days
    moonAge = fmod(daysDelta, period)

    phaseNum = int( moonAge / (period / len(phases)) )

    return phases[phaseNum]

    if __name__ == "__main__":
    samples = [date(2000, 1, 1), date(2000, 1, 6), date(2000, 2, 8), date.today()]
    for sample in samples:
    print sample, "=", phaseOfMoon(sample)

  4. I didn’t have any time for actually coding this up, but John Conway has a way of computing the phase of the moon that can be done in your head as part of his (formerly two vollume, now republished by AK Peterson as 4 volume) series Winning Ways. The “nice” feature of this is that it doesn’t actually require any conversion to julian dates. I think that most of these proposed solutions are slightly off: some don’t perform proper rounding, and most seem to ignore the fact that the first new moon of January 2000 was January 6,

  5. Sorry, hit the wrong button.

    … was January 6, 18:14 UTC.

  6. Maurits said

    A quibble – the phases of the moon are only four:

    1) waxing crescent
    2) waxing gibbous
    3) waning gibbous
    4) waning crescent

    The other four states (new, first quarter, full, third quarter) are more in the nature of events than phases… that is to say, the moon is only new for an instant.

  7. Mike said
    import datetime
    from math import fmod
    
    def moonphase( date=None ):
        '''return phase of the moon on a given date.
    
        date is a datetime.date object.  Defaults to today.
        '''
        
        if date is None:
            date = datetime.date.today()
    
        known_new_moon = datetime.date(2000,1,6)
        lunar_cycle = 29.530588853   # days per lunation
        phase_length = lunar_cycle/8 # days per phase
    
        days = (date - known_new_moon).days
        days = fmod(days + phase_length/2, lunar_cycle)
    
        phase = int( days * ( 8 / lunar_cycle ) )
    
        phasetext = ("new", "waxing crescent",
                     "first quarter", "waxing gibbous",
                     "full", "waning gibbous",
                     "last quater", "waning crescent")
    
        return phase, phasetext[ phase ]
    
  8. David said

    Erlang version

    -module(lunar).
    -export([phase/3]).
    
    -define(LUNAR_CYCLE, 29.530588853).
    
    jday(Year, Month, Day) ->
        A = (14 - Month) div 12,
        Y = Year + 4800 - A,
        M = Month + 12*A - 3,
        Day + (153*M + 2) div 5 + 365*Y + Y div 4 - Y div 100 + Y div 400 - 32045.
    
    fmod(A, B) -> A - B*trunc(A/B).
    
    phase(X) when X <  1.84566 -> new;
    phase(X) when X <  5.53699 -> waxing_crescent;
    phase(X) when X <  9.22831 -> first_quarter;
    phase(X) when X < 12.91963 -> waxing_gibbous;
    phase(X) when X < 16.61096 -> full;
    phase(X) when X < 20.30228 -> waning_gibbous;
    phase(X) when X < 23.99361 -> last_quarter;
    phase(X) when X < 27.68493 -> waning_crescent;
    phase(_) -> new.
    
    phase(Year, Month, Day) ->
        Days = jday(Year, Month, Day) - jday(2000, 1, 6),
        Offset = fmod(Days, ?LUNAR_CYCLE),
        phase(Offset).
    
    8> lunar:phase(2014,10,11).
    waning_gibbous
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 628 other followers

%d bloggers like this: