満月・新月の日時をRubyで計算

太陽と月の地心直交座標系での位置を計算
黄道座標系に変換して太陽と月の黄経を計算
黄経差を計算
黄経差が0°や180°になるタイミングを計算

まずJPLからデータをダウンロード。

$ wget -r 'ftp://ssd.jpl.nasa.gov/pub/eph/planets/Linux/de430/linux_p1550p2650.430'

Gemfile は以下。eph_jplはJPLのデータを扱うためです。mk_timeはUTCからTTに時刻系を変換するために使っています。

source "https://rubygems.org"

gem 'mk_time'
gem 'eph_jpl'

Rubyのソースコードは以下です。

require 'time'
require 'mk_time'
require 'eph_jpl'

# JPLのデータ
$data_path = 'ssd.jpl.nasa.gov/pub/eph/planets/Linux/de430/linux_p1550p2650.430'

$pi57 = 180.0 / Math::PI
$pi2 = 2.0 * Math::PI

def ephXyz(target, time)
  # target: sun / moon
  # time: UTC DateTime型
  ttjd = MkTime::Calc.new(MkTime::Calc.new(time).tt).jd # UTCをTTに変換したユリウス日
  if target == :sun
    target_id = 11
  elsif target == :moon
    target_id = 10
  else
    raise
  end
  ephJpl = EphJpl.new($data_path, target_id, 3, ttjd)
  value = ephJpl.calc
  x = value[0]
  y = value[1]
  z = value[2]
  return [x, y, z]
end

def xyzEquatorialToEcliptic(xyz)
  # 地球の黄道傾斜角で座標を回転
  # バイアス・歳差・章動は考慮していない
  return rotationX(xyz, -0.4090926)
end

def rotationX(xyz, th)
  # X軸を中心に位置を反時計回りに回転
  x = xyz[0]
  y = xyz[1]
  z = xyz[2]
  cos = Math.cos(th)
  sin = Math.sin(th)
  return [x, y * cos - z * sin, y * sin + z * cos]
end

# 計算する期間
startTime = Time.parse("2021-01-01T00:00:00.000Z")
endTime = Time.parse("2021-02-01T00:00:00.000Z")

time = startTime
moonSunLng = nil

# 1時間ごとに月と太陽の黄経差を計算
while time < endTime
  sun = xyzEquatorialToEcliptic(ephXyz(:sun, time))
  moon = xyzEquatorialToEcliptic(ephXyz(:moon, time))

  # 黄経を計算
  sunLng = Math.atan2(sun[1], sun[0]);
  moonLng = Math.atan2(moon[1], moon[0]);

  # 黄経差を計算
  moonSunLng2 = moonLng - sunLng
  moonSunLng2 += $pi2 if (moonSunLng2 < 0)
  moonSunLng2 *= $pi57

  if moonSunLng != nil
    # 上弦、満月、下弦、新月の時刻を線形補間で計算
    if moonSunLng <= 90 && moonSunLng2 > 90
      t = time - (moonSunLng2 - 90) / (moonSunLng2 - moonSunLng) * 3600
      puts("#{t} 上弦")
    elsif moonSunLng <= 180 && moonSunLng2 > 180
      t = time - (moonSunLng2 - 180) / (moonSunLng2 - moonSunLng) * 3600
      puts("#{t} 満月")
    elsif moonSunLng <= 270 && moonSunLng2 > 270
      t = time - (moonSunLng2 - 270) / (moonSunLng2 - moonSunLng) * 3600
      puts("#{t} 下弦")
    elsif moonSunLng > 270 && moonSunLng2 < 90
      t = time - moonSunLng2 / (moonSunLng2 - moonSunLng + 360) * 3600
      puts("#{t} 新月")
    end
  end

  moonSunLng = moonSunLng2
  #puts("#{time} #{moonSunLng}")

  time = time + 3600
end

実行。

$ bundle install

$ bundle exec ruby sample.rb

実行結果。

2021-01-06 09:37:51 UTC 下弦
2021-01-13 05:00:48 UTC 新月
2021-01-20 21:02:22 UTC 上弦
2021-01-28 19:16:53 UTC 満月

時刻はUTCです。国立天文台天文情報センター暦計算室のページで満月などの時刻を調べると、分単位で一致しています。