#!/usr/bin/env ruby
require "pp"
module MoonPhase
  MOON_RAD = (Math::PI/180.0)
  MOON_SMALL_FLOAT = 1e-12

  #
  # http://www.voidware.com/moon_phase.htm
  # Math and numbers from site above
  #
  def MoonPhase.moon_phase(time)
      #/*
      #  Calculates more accurately than Moon_phase , the phase of the moon at
      #  the given epoch.
      #  number might not be 0 or 1 because of rounding issues
      #  returns the moon phase as a real number (0-1)
      #  */

      julian = julian_days_float(time)-2444238.5;

      time = moon_position(julian) - sun_position(julian)
      time += 360 if time < 0

      ((1.0 - Math.cos((time)*MOON_RAD))/2)
  end

  def moon_phase
    MoonPhase.moon_phase(self)
  end

  private

  #Returns the number of julian days for the specified day.
  def MoonPhase.julian_days_float(time)
      year = time.year
      month = time.month
      day = time.day
      b = 0
      if (month < 3)
          year -= 1
          month += 12
      end

      if (year > 1582 || (year == 1582 && month > 10) || (year == 1582 && month == 10 && day > 15))
          a = year / 100;
          b = 2 - a + a/4;
      end

      c = 365.25*year;
      e = 30.6001*(month+1);
      b + c + e + day + 1720994.5;
  end

  def MoonPhase.sun_position(julian_days)

      n = 360 / 365.2422 * julian_days
      n -= (n/360).to_i * 360.0
      x = n - 3.762863
      x += 360 if (x<0)
      x = x * MOON_RAD
      e = x
      begin
          dl = e - 0.016718 * Math.sin(e) - x
          e -=  dl / (1 - 0.016718 * Math.cos(e))
      end while (dl.abs >= MOON_SMALL_FLOAT)
      v = 360 / Math::PI * Math.atan(1.01686011182 * Math.tan(e/2))
      l = v + 282.596403
      l = l - (l / 360).to_i * 360.0
  end

  def MoonPhase.moon_position(julian, ls = nil)
      ls = sun_position(julian) unless ls

      ms = (0.985647332099 * julian) - 3.762863
      ms += 360.0 if ms < 0

      l = (13.176396 * julian) + 64.975464
      l -= ((l/360).to_i * 360.0)
      l += 360.0 if (l < 0)

      mm = l - 0.1114041 * julian - 349.383063
      mm -= (mm/360).to_i * 360.0

      #n is not used i dont think.
      n = 151.950429 - 0.0529539 * julian
      n -= (n/360).to_i * 360.0

      ev = 1.2739 * Math.sin( (2 * (l - ls) - mm) * MOON_RAD)
      sms = Math.sin(ms * MOON_RAD)
      ae = 0.1858 * sms
      mm += ev - ae - 0.37 * sms
      ec = 6.2886 * Math.sin(mm * MOON_RAD)
      l += ev + ec- ae + 0.214 * Math.sin(2 * mm * MOON_RAD)
      0.6583 * Math.sin( 2 * (l - ls) * MOON_RAD) + l;
  end

end

#pp MoonPhase.moon_phase((Time.now+(-3*24*60*60)).utc)