sunrise_sunset.rb

require 'date'

# Ruby implementation of algorithm for approximating
# sunrise and sunset for a given date, latitude, and longitude.
# Ported from instructions at http://users.electromagnetic.net/bu/astro/sunrise-set.php 
# I think this is accurate to within about 10 minutes.
class SunriseSunset
  def sunset
    @sunset ||= begin
      jd = (julian_sunset + 0.5).to_i
      seconds = (julian_sunset - jd + 0.5) * 24 * 60 * 60
      date = Date.jd(jd)
      (Time.utc(date.year, date.month, date.day) + seconds).localtime
    end
  end

  def sunrise
    @sunrise ||= begin
      jd = (julian_sunrise + 0.5).to_i
      seconds = (julian_sunrise - jd + 0.5) * 24 * 60 * 60
      date = Date.jd(jd)
      (Time.utc(date.year, date.month, date.day) + seconds).localtime
    end
  end

  private

  def initialize(lat, lng, date = Date.today)
    @lat, @lng, @date = lat, lng, date
  end

  def julian_cycle
    @julian_cycle ||= ((@date.jd - 2451545 - 0.0009) - (lng_w/360)).round
  end

  def approximate_solar_noon
    @approximate_solar_noon = 2451545 + 0.0009 + (lng_w/360) + julian_cycle
  end

  def mean_solar_anomaly
    @mean_solar_anomaly = (357.5291 + 0.98560028 * (approximate_solar_noon - 2451545)) % 360 
  end

  def equation_of_center
    @equation_of_center ||= 1.9148*degree_sin(mean_solar_anomaly) + 0.0200*degree_sin(2*mean_solar_anomaly) + 0.0003*degree_sin(3*mean_solar_anomaly)
  end

  def ecliptical_longitude
    @ecliptical_longitude ||= (mean_solar_anomaly + 102.9372 + equation_of_center + 180) % 360
  end

  def solar_noon
    @solar_noon ||= approximate_solar_noon + 0.0053*degree_sin(mean_solar_anomaly) - 0.0069*degree_sin(2*ecliptical_longitude)
  end

  def declination
    @declination ||= degree_asin(degree_sin(ecliptical_longitude) * degree_sin(23.45))
  end

  def hour_angle
    @hour_angle ||= degree_acos((degree_sin(-0.83) - degree_sin(lat_n) * degree_sin(declination))/(degree_cos(lat_n) * degree_cos(declination)))
  end

  def j_star_star
    @j_star_star ||= 2451545 + 0.0009 + ((hour_angle + lng_w)/360) + julian_cycle
  end

  def julian_sunset
    @julian_sunset ||= j_star_star + 0.0053 * Math.sin(mean_solar_anomaly) - 0.0069 * Math.sin(2*ecliptical_longitude)
  end

  def julian_sunrise
    @julian_sunrise ||= 2*solar_noon - julian_sunset
  end

  def lng_w
    -@lng
  end

  def lat_n
    @lat
  end

  def degree_sin(degrees)
    Math.sin(degrees / 360 * 2 * Math::PI)
  end

  def degree_asin(sin)
    Math.asin(sin) / (2 * Math::PI) * 360
  end

  def degree_cos(degrees)
    Math.cos(degrees / 360 * 2 * Math::PI)
  end

  def degree_acos(cos)
    Math.acos(cos) / (2 * Math::PI) * 360
  end
end