class SunCalc

# File lib/sun_calc.rb, line 67
def self.add_sun_time(angle, rise_name, set_name)
  @__sun_times__.push([angle, rise_name, set_name])
end

# File lib/sun_calc.rb, line 158
def self.moon_illumination(date = Time.now)
  d = to_days(date)
  s = sun_coords(d)
  m = moon_coords(d)
  sdist = 149_598_000 # Distance from Earth to Sun in kilometers
  phi = Math.acos(Math.sin(s[:dec]) * Math.sin(m[:dec]) +
             Math.cos(s[:dec]) * Math.cos(m[:dec]) * Math.cos(s[:ra] -
                                                              m[:ra]))
  inc = Math.atan2(sdist * Math.sin(phi), m[:dist] - sdist * Math.cos(phi))
  angle = Math.atan2(Math.cos(s[:dec]) * Math.sin(s[:ra] - m[:ra]),
                     Math.sin(s[:dec]) * Math.cos(m[:dec]) -
                       Math.cos(s[:dec]) * Math.sin(m[:dec]) *
                         Math.cos(s[:ra] - m[:ra]))
  { fraction: (1 + Math.cos(inc)) / 2,
    phase: 0.5 + 0.5 * inc * (angle < 0 ? -1 : 1) / Math::PI,
    angle: angle }
end

# File lib/sun_calc.rb, line 72
def self.moon_position(date, lat, lng)
  lw  = ONE_RADIAN * -lng
  phi = ONE_RADIAN * lat
  d   = to_days(date)
  c   = moon_coords(d)
  h_  = sidereal_time(d, lw) - c[:ra]
  h   = altitude(h_, phi, c[:dec])
  # Formula 14.1 from "Astronomical Algorithms" 2nd edition by Jean Meeus
  # (Willmann-Bell, Richmond) 1998.
  pa  = Math.atan2(Math.sin(h_),
                   Math.tan(phi) * Math.cos(c[:dec]) -
                     Math.sin(c[:dec]) * Math.cos(h_))
  h += astro_refraction(h) # altitude correction for refraction

  { azimuth: azimuth(h_, phi, c[:dec]),
    altitude: h,
    distance: c[:dist],
    parallacticAngle: pa }
end

# File lib/sun_calc.rb, line 96
def self.moon_times(date, lat, lng)
  t = Time.utc(date.year, date.month, date.day)
  hc = 0.133 * ONE_RADIAN
  h0 = SunCalc.moon_position(t, lat, lng)[:altitude] - hc
  ye = 0
  max = nil
  min = nil
  rise = nil
  set = nil
  # Iterate in 2-hour steps checking if a 3-point quadratic curve crosses zero
  # (which means rise or set). Assumes x values -1, 0, +1.
  (1...24).step(2).each do |i|
    h1 = SunCalc.moon_position(hours_later(t, i), lat, lng)[:altitude] - hc
    h2 = SunCalc.moon_position(
      hours_later(t, i + 1), lat, lng
    )[:altitude] - hc
    a = (h0 + h2) / 2 - h1
    b = (h2 - h0) / 2
    xe = -b / (2 * a)
    ye = (a * xe + b) * xe + h1
    d = b * b - 4 * a * h1
    roots = 0
    x1 = 0
    x2 = 0
    min = i + xe if xe.abs <= 1 && ye < 0
    max = i + xe if xe.abs <= 1 && ye > 0
    if d >= 0
      dx = Math.sqrt(d) / (a.abs * 2)
      x1 = xe - dx
      x2 = xe + dx
      roots += 1 if x1.abs <= 1
      roots += 1 if x2.abs <= 1
      x1 = x2 if x1 < -1
    end
    if roots == 1
      if h0 < 0
        rise = i + x1
      else
        set = i + x1
      end
    elsif roots == 2
      rise = i + (ye < 0 ? x2 : x1)
      set = i + (ye < 0 ? x1 : x2)
    end
    break if rise && set && min && max
    h0 = h2
  end
  {}.tap do |result|
    result[:nadir] = hours_later(t, min) if min
    result[:lunar_noon] = hours_later(t, max) if max
    result[:moonrise] = hours_later(t, rise) if rise
    result[:moonset] = hours_later(t, set) if set
    result[ye > 0 ? :always_up : :always_down] = true if !rise && !set
  end
end

# File lib/sun_calc.rb, line 29
def self.sun_position(date, lat, lng)
  lw  = ONE_RADIAN * -lng
  phi = ONE_RADIAN * lat
  d   = to_days(date)
  c   = sun_coords(d)
  h   = sidereal_time(d, lw) - c[:ra]
  { azimuth: azimuth(h, phi, c[:dec]),
    altitude: altitude(h, phi, c[:dec]) }
end

# File lib/sun_calc.rb, line 40
def self.sun_times(date, lat, lng)
  lw     = ONE_RADIAN * -lng
  phi    = ONE_RADIAN * lat
  d      = to_days(date)
  n      = julian_cycle(d, lw)
  ds     = approx_transit(0, lw, n)
  m      = solar_mean_anomaly(ds)
  l      = ecliptic_longitude(m)
  dec    = declination(l, 0)
  j_noon = solar_transit_j(ds, m, l)
  { solar_noon: from_julian(j_noon),
    nadir: from_julian(j_noon - 0.5) }.tap do |result|
    @__sun_times__.each do |time|
      begin
        j_set = get_set_j(time[0] * ONE_RADIAN, lw, phi, dec, n, m, l)
        j_rise = j_noon - (j_set - j_noon)
        result[time[1].to_sym] = from_julian(j_rise)
        result[time[2].to_sym] = from_julian(j_set)
      rescue Math::DomainError
        result[time[1].to_sym] = nil
        result[time[2].to_sym] = nil
      end
    end
  end
end