Moon az/el calculations

moon.rb

# Source: https://github.com/gregseth/suncalc-php/blob/master/suncalc.php

require 'date'
RAD = Math::PI / 180.0
DEG = 180.0 / Math::PI
E = RAD * 23.4397 # obliquity of the Earth
J2000 = 2451545

def to_days(datetime)
  datetime.ajd.to_f - J2000
end

def declination(l, b)
  Math.asin(Math.sin(b) * Math.cos(E) + Math.cos(b) * Math.sin(E) * Math.sin(l))
end

def right_ascension(l, b)
  Math.atan2(Math.sin(l) * Math.cos(E) - Math.tan(b) * Math.sin(E), Math.cos(l))
end

def sidereal_time(d, lw)
  RAD * (280.16 + 360.9856235 * d) - lw
end

def altitude(h, phi, dec)
  Math.asin(Math.sin(phi) * Math.sin(dec) + Math.cos(phi) * Math.cos(dec) * Math.cos(h))
end

def azimuth(h, phi, dec)
  Math.atan2(Math.sin(h), Math.cos(h) * Math.sin(phi) - Math.tan(dec) * Math.cos(phi))
end

def moon_coords(d)
  el = (218.316 + 13.176396 * d) * RAD
  m = (134.963 + 13.064993 * d) * RAD
  f = (93.272 + 13.229350 * d) * RAD

  l = el + 6.289 * RAD * Math.sin(m) # latitude
  b = 5.128 * RAD * Math.sin(f) # longitude
  dt = 385001 - 20905 * Math.cos(m) # distance to the moon in km

  {
    declination: declination(l, b),
    right_ascension: right_ascension(l, b),
    distance: dt
  }
end

def moon_position(lat_deg, lon_deg, datetime)
  lw = -lon_deg * RAD
  phi = lat_deg * RAD
  d = to_days(datetime)

  c = moon_coords(d)
  big_h = sidereal_time(d, lw) - c[:right_ascension]
  h = altitude(big_h, phi, c[:declination])

  # Altitude correction for refraction
  h = h + RAD * 0.017 / Math.tan(h + RAD * 10.26 / (h + RAD * 5.10))

  {
    azimuth: azimuth(big_h, phi, c[:declination]),
    elevation: h,
    distance: c[:distance]
  }
end

my_lat = ARGV[0].to_f # 41.8640081
my_lon = ARGV[1].to_f # -72.1238411

pos = moon_position(my_lat, my_lon, DateTime.now)
puts "Azimuth:   #{((pos[:azimuth] * DEG + 180) % 360).round(1)}°"
puts "Elevation: #{(pos[:elevation] * DEG).round(1)}°"