太陽の位置を計算するプログラムを書きました。

JPL(NASA Jet Propulsion Laboratory、NASAジェット推進研究所)が太陽や月、惑星の位置を計算するためのデータを公開しています。そのデータを読み込んで計算するRuby Gemがありますので、それを呼び出すだけです。

JPLのデータは、座標系は国際天文準拠系で、時刻系はTDB（太陽系力学時）です。

時刻系については前回の記事参照。

国際天文準拠系は、太陽系重心を原点にし、銀河系外天体を観測して得られる固定した方向を座標軸とした、回転しない慣性系です。J2000.0とほとんど一致するように決められているようです。

理科年表オフィシャルサイト／バックナンバー／平成13年／トピックス／天文部　現在の天文基準座標系

JPLのデータは西暦1550年から2549年までの天体の位置を計算できる係数が含まれています。1000年の期間を最大32日の期間で区切って、期間ごとにチェビシェフ級数展開の係数が与えられています。

TDB（太陽系力学時）での時刻から該当期間の係数を取り出し、チェビシェフ級数展開を計算すれば、国際天文準拠系での位置を求めることができます。

太陽系重心が原点になっていますので、地心にするには地球の位置も計算して引き算すればよさそうです。

eph_jpl というRuby Gemがそれをしているようです。以下のブログ記事の方がRuby Gemを作成されてました。

サンプルコード

JPLのデータは以下のコマンドでダウンロードできます。

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

作成したファイルは Gemfile と sample.rb の2つ。

Gemfile

source "https://rubygems.org"

gem 'mk_time'
gem 'eph_jpl'

sample.rb

require 'mk_time'
require 'eph_jpl'

$data_path = 'ssd.jpl.nasa.gov/pub/eph/planets/Linux/de430/linux_p1550p2650.430'

def outputEph(obj)
  puts("#{obj.center} -> #{obj.target}")
  puts("#{obj.center_name} -> #{obj.target_name}")
  puts("Julian Day: #{obj.jd}")
  puts("km flag: #{obj.km}")
  puts("unit: #{obj.unit}")
  value = obj.calc
  p(value)
  if obj.unit == "AU, AU/day"
    x = value[0]
    y = value[1]
    z = value[2]
    r = Math.sqrt(x * x + y * y + z * z)
    puts("r: #{r} AU")
    puts("r: #{r * 149_597_870.700} km")
  end
  puts()
end

def output(time_str)
  t = MkTime.new(time_str) # UTC
  tt = MkTime::Calc.new(t.tt) # TT
  puts(t.tt)
  puts()

  # 太陽
  obj = EphJpl.new($data_path, 11, 3, tt.jd)
  outputEph(obj)

  # 月
  obj = EphJpl.new($data_path, 10, 3, tt.jd)
  outputEph(obj)
end

output("20210101000000")   # UTCで2021年1月1日0時ちょうど
output("20201231235850816") # TTで2021年1月1日0時ちょうど

JPLのデータからの計算は太陽系力学時(TDB)を使うべきですが、TDBは地球時(TT)とほとんど同じなので、UTCからTTに変換したものをTDBの代わりに使いました。

地球時(TT)と太陽系力学時(TDB)の関係式

実行

$ bundle exec ruby sample.rb

以下は実行結果にコメントを入れたものです。

# UTC 2021年1月1日0時ちょうど
2021-01-01 00:01:09 +0000 # この表記はTT

# 地球から見た太陽の位置
3 -> 11
Earth -> Sun
Julian Day: 2459215.500800741
km flag: false
unit: AU, AU/day
[0.1791297879205522, -0.8870498072939582, -0.3845320934488644, 0.01719101705001951, 0.0029314742193810215, 0.0012712387581410812]
r: 0.9832648536548414 AU
r: 147094328.44091138 km

# 地球から見た月の位置
3 -> 10
Earth -> Moon
Julian Day: 2459215.500800741
km flag: false
unit: AU, AU/day
[-0.0013833376208353678, 0.0019323528675363104, 0.0010131350936447366, -0.0004831632689964738, -0.00032365358794430693, -9.882918359702189e-05]
r: 0.0025834189160623537 AU
r: 386473.96896903013 km

# TT 2021年1月1日0時ちょうど
# UTC表記にすると0時より69秒前
2021-01-01 00:00:00 +0000

# 地球から見た太陽の位置
3 -> 11
Earth -> Sun
Julian Day: 2459215.5
km flag: false
unit: AU, AU/day
[0.1791160223524058, -0.8870521545571489, -0.38453311134357027, 0.017191060993880306, 0.002931254100533473, 0.0012711434267272377]
r: 0.9832648616153199 AU
r: 147094329.631782 km

# 地球から見た月の位置
3 -> 10
Earth -> Moon
Julian Day: 2459215.5
km flag: false
unit: AU, AU/day
[-0.0013829507090954102, 0.0019326119983273523, 0.0010132142134440024, -0.00048322107911419196, -0.00032357397039131907, -9.87872953448828e-05]
r: 0.0025834366340576824 AU
r: 386476.6195434044 km

地球時(TT)2021年1月1日0時ちょうどの値は、天文年鑑に載っている有効桁数7桁の範囲で一致しました。

suzuki-navi’s blog

太陽の地心直交座標系での位置計算

サンプルコード

実行