Skip to content

Instantly share code, notes, and snippets.

@komasaru komasaru/jpl_calc_de430.rb
Last active Apr 19, 2018

Embed
What would you like to do?
Ruby script to calc ICRS coordinates with JPL binary data DE430.
#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= JPLEPH(JPL の DE430 バイナリデータ)読み込み、座標(位置・速度)を計算する
#
# date name version
# 2016.03.12 mk-mode.com 1.00 新規作成
#
# Copyright(C) 2016 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# 【引数】
# [第1] 対象天体番号(必須)
# 1: 水星 (Mercury)
# 2: 金星 (Venus)
# 3: 地球 (Earth)
# 4: 火星 (Mars)
# 5: 木星 (Jupiter)
# 6: 土星 (Saturn)
# 7: 天王星 (Uranus)
# 8: 海王星 (Neptune)
# 9: 冥王星 (Pluto)
# 10: 月 (Moon)
# 11: 太陽 (Sun)
# 12: 太陽系重心 (Solar system Barycenter)
# 13: 地球 - 月の重心 (Earth-Moon Barycenter)
# 14: 地球の章動 (Earth Nutations)
# 15: 月の秤動 (Lunar mantle Librations)
# [第2] 基準天体番号(必須。 0, 1 - 13)
# ( 0 は、対象天体番号が 14, 15 のときのみ)
# [第3] ユリウス日(省略可。省略時は現在日時のユリウス日)
#---------------------------------------------------------------------------------
# 【注意事項】
# * 求める座標は「赤道直角座標(ICRS)」
# * 天体番号は、係数データの番号(並び順)と若干異なるので注意する。
# (特に、天体番号 3, 10, と 12 以降)
# 係数データの並び順:
# 水星(1)、金星(2)、地球 - 月の重心(3)、火星(4)、木星(5)、土星(6)、
# 天王星(7)、海王星(8)、冥王星(9)、月(地心)(10)、太陽(11)、
# 地球の章動(12)、月の秤動(13)
# * 時刻系は「太陽系力学時(TDB)」である。(≒地球時(TT))
# * 天体番号が 1 〜 13 の場合は、 x, y, z の位置・速度(6要素)、
# 天体番号が 14 の場合は、黄経における章動 Δψ, 黄道傾斜における章動 Δε の
# 角位置・角速度(4要素)、
# 天体番号が 15 の場合は、 φ, θ, ψ の角位置・角速度(6要素)。
# * 対象天体番号 = 基準天体番号 は、無意味なので処理しない。
# * 天体番号が 12 の場合は、 x, y, z の位置・速度の値は全て 0.0 とする。
# * その他、JPL 提供の FORTRAN プログラム "testeph.f" を参考にした。
#---------------------------------------------------------------------------------
#++
class JplCalcDe430
USAGE = <<-"EOS"
【使用方法】 ./jpl_calc_de430.rb 対象天体番号 基準天体番号 [ユリウス日]
【天体番号】(対象天体番号: 1 - 15, 基準天体番号: 0 - 13)
1: 水星, 2: 金星, 3: 地球, 4: 火星, 5: 木星,
6: 土星, 7: 天王星, 8: 海王星, 9: 冥王星, 10: 月,
11: 太陽, 12: 太陽系重心, 13: 地球 - 月の重心,
14: 地球の章動, 15: 月の秤動,
0: 対象天体番号 14, 15 の場合に指定
※対象天体番号≠基準天体番号
EOS
FILE_BIN = "JPLEPH"
KSIZE = 2036
RECL = 4
ASTRS = [
"Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus",
"Neptune", "Pluto", "Moon", "Sun", "Solar system Barycenter",
"Earth-Moon barycenter", "Earth Nutations", "Lunar mantle Librations"
]
KIND = 2 # 計算区分(0: 計算しない、1: 位置のみ計算、2: 位置・速度を計算)
# (但し、現時点では 1 の場合の処理は実装していない)
BARY = true # 基準フラグ(true: 太陽系重心が基準 false: 太陽が基準)
KM = false # 単位フラグ(true: km, km/sec, false: AU, AU/day)
def initialize
get_args # コマンドライン引数取得
@pos = 0 # レコード位置
@list = Array.new(12, 0) # 計算対象フラグ一覧配列
@ipts = Array.new # 定数一覧配列
@coeffs = Array.new # 係数一覧配列
@idx_subs = Array.new(12) # サブ区間のインデックス
@tcs = Array.new(12) # サブ区間内の時間(チェビシェフ多項式用に正規化)
@pvs = Array.new(11).map { |a| Array.new(6, 0.0) } # 位置・角度データ配列
@pvs_2 = Array.new(13) # 位置・角度データ配列(対象天体と基準天体の差算出用)
@rrds = Array.new(6, 0.0) # 算出データ(対象天体と基準天体の差)配列
end
def exec
begin
# === ヘッダ情報取得
get_ttl # TTL(タイトル)
get_cnam # CNAM(定数名)
get_ss # SS(ユリウス日(開始、終了)、分割日数)
check_jd # 引数のユリウス日をチェック
get_ncon # NCON(定数の数)
get_au # AU(天文単位)
get_emrat # EMRAT(地球と月の質量比)
get_ipt # IPT(オフセット、係数の数、サブ区間数)(水星〜地球の章動)
get_numde # NUMDE(DEバージョン番号)取得
get_ipt_13 # IPT(オフセット、係数の数、サブ区間数)(月の秤動)
get_cval # CVAL(定数値)
get_jdepoc # JDEPOC(元期)
# === 各種計算
get_list # 計算対象フラグ一覧(係数データの並びに対応)取得
get_coeff # 係数取得(対象、基準天体のインデックス分を取得)
#p @ttl
#p @cnams
#p @sss
#p @ncon
#p @au
#p @emrat
#p @ver_de
#p @ipts
#p @cvals
#p @jdepoc
#p @list
#p @coeffs
#exit
# 補間(11:太陽)
@pv_sun = interpolate(11)
# 補間(1:水星〜10:月)
0.upto(9) do |i|
next if @list[i] == 0
@pvs[i] = interpolate(i + 1)
next if i > 8
next if BARY
@pvs[i] = @pvs[i].map.with_index do |pv, j|
pv - @pv_sun[j]
end
end
# 補間(14:地球の章動)
if @list[10] > 0 && @ipts[11][1] > 0
@p_nut = interpolate(14)
end
# 補間(15:月の秤動)
if @list[11] > 0 && @ipts[12][1] > 0
@pvs[10] = interpolate(15)
end
# 対象天体と基準天体の差
case
when @astrs[0] == 14
@rrds = @p_nut + [0.0, 0.0] if @ipts[11][1] > 0
when @astrs[0] == 15
@rrds = @pvs[10] if @ipts[12][1] > 0
else
0.upto(9) { |i| @pvs_2[i] = @pvs[i] }
@pvs_2[10] = @pv_sun if @astrs.include?(11)
@pvs_2[11] = Array.new(6, 0.0) if @astrs.include?(12)
@pvs_2[12] = @pvs[2] if @astrs.include?(13)
if @astrs[0] * @astrs[1] == 30 || @astrs[0] + @astrs[1] == 13
@pvs_2[2] = Array.new(6, 0.0)
else
@pvs_2[2] = @pvs[2].map.with_index do |pv, i|
pv - @pvs[9][i] / (1.0 + @emrat)
end unless @list[2] == 0
@pvs_2[9] = @pvs_2[2].map.with_index do |pv, i|
pv + @pvs[9][i]
end unless @list[9] == 0
end
0.upto(5) do |i|
@rrds[i] = @pvs_2[@astrs[0] - 1][i] - @pvs_2[@astrs[1] - 1][i]
end
end
# === 結果出力
display
rescue => e
$stderr.puts "[#{e.class}] #{e.message}"
e.backtrace.each { |tr| $stderr.puts "\t#{tr}"}
exit 1
end
end
private
#=========================================================================
# コマンドライン引数取得
#
# * 第1引数の値をインスタンス変数 @astrs[0] に設定する
# * 第2引数の値をインスタンス変数 @astrs[1] に設定する
# * 第3引数の値をインスタンス変数 @jd に設定する
# (但し、第3引数が存在しない場合は現在日時を設定する)
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_args
# 対象、基準天体番号
unless ARGV[0] && ARGV[1]
puts USAGE
exit 0
end
target = ARGV[0].to_i
center = ARGV[1].to_i
if (target < 1 || target > 15) ||
(center < 0 || center > 13) ||
target == center ||
(target == 14 && center != 0) ||
(target == 15 && center != 0) ||
(target != 14 && target != 15 && center == 0)
puts USAGE
exit 0
end
@astrs = [target, center]
# ユリウス日
@jd = ARGV[2]
@jd = @jd ? @jd.to_f : gc_to_jd({
year: Time.now.strftime("%Y").to_i,
month: Time.now.strftime("%m").to_i,
day: Time.now.strftime("%d").to_i,
hour: Time.now.strftime("%H").to_i,
min: Time.now.strftime("%M").to_i,
sec: Time.now.strftime("%S").to_i
})
rescue => e
raise
end
#=========================================================================
# 年月日(グレゴリオ暦)からユリウス日(JD)を計算する
#
# * フリーゲルの公式を使用する
# JD = int(365.25 * year)
# + int(year / 400)
# - int(year / 100)
# + int(30.59 (month - 2))
# + day
# + 1721088
# * 上記の int(x) は厳密には、 x を超えない最大の整数
# * 「ユリウス日」でなく「準ユリウス日」を求めるなら、
# `+ 1721088` を `- 678912` とする。
#
# @params: hash = {
# year(int), month(int), day(int), hour(int), min(int), sec(int)
# }
# @return: jd (ユリウス日)
#=========================================================================
def gc_to_jd(hash)
year = hash[:year]
month = hash[:month]
day = hash[:day]
hour = hash[:hour]
min = hash[:min]
sec = hash[:sec]
begin
# 1月,2月は前年の13月,14月とする
if month < 3
year -= 1
month += 12
end
# 日付(整数)部分計算
jd = (365.25 * year).floor
jd += (year / 400.0).floor
jd -= (year / 100.0).floor
jd += (30.59 * (month - 2)).floor
jd += day
jd += 1721088
# 時間(小数)部分計算
t = sec / 3600.0
t += min / 60.0
t += hour
t /= 24.0
return jd + t
rescue => e
raise
end
end
#=========================================================================
# 引数のユリウス日をチェック
#
# * @jd の値をヘッダ情報のユリウス日(開始、終了)と比較・チェックする
#
# @params: <none>
# @return: <none>
#=========================================================================
def check_jd
if @jd < @sss[0] || @jd >= @sss[1]
puts "Please input JD s.t. #{@sss[0]} <= JD < #{@sss[1]}."
exit 0
end
rescue => e
raise
end
#=========================================================================
# TTL 取得
#
# * タイトルをインスタンス変数 @ttl に設定する
# * 84 byte * 3
# * ASCII文字列(スペースを詰める/後続するnull文字やスペースを削除)
# * 取得したバイト数だけレコード位置 @pos を進める
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_ttl
recl = 84
begin
@ttl = (0..2).map do |i|
File.binread(FILE_BIN, recl, @pos + recl * i).unpack("A*")[0]
end.join("\n")
@pos += recl * 3
rescue => e
raise
end
end
#=========================================================================
# CNAM 取得
#
# * 定数名をインスタンス変数(配列) @cnams に設定する
# * 6 byte * 400
# * ASCII文字列(スペースを詰める/後続するnull文字やスペースを削除)
# * 取得したバイト数だけレコード位置 @pos を進める
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_cnam
recl = 6
begin
@cnams = (0..399).map do |i|
File.binread(FILE_BIN, recl, @pos + recl * i).unpack("A*")[0]
end
@pos += recl * 400
rescue => e
raise
end
end
#=========================================================================
# SS 取得
#
# * ユリウス日(開始、終了)、分割日数をインスタンス変数(配列) @sss に
# 設定する
# * 8 byte * 3
# * 倍精度浮動小数点数(機種依存)
# * 取得したバイト数だけレコード位置 @pos を進める
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_ss
recl = 8
begin
@sss = (0..2).map do |i|
File.binread(FILE_BIN, recl, @pos + recl * i).unpack("d*")[0]
end
@pos += recl * 3
rescue => e
raise
end
end
#=========================================================================
# NCON 取得
#
# * 定数の数をインスタンス変数 @ncon に設定する
# * 4 byte * 1
# * unsigned int (符号なし整数, エンディアンと int のサイズに依存)
# * 取得したバイト数だけレコード位置 @pos を進める
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_ncon
recl = 4
begin
@ncon = File.binread(FILE_BIN, recl, @pos).unpack("I*")[0]
@pos += recl
rescue => e
raise
end
end
#=========================================================================
# AU 取得
#
# * AU(天文単位)の値をインスタンス変数 @au に設定する
# * 8 byte * 1
# * 倍精度浮動小数点数(機種依存)
# * 取得したバイト数だけレコード位置 @pos を進める
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_au
recl = 8
begin
@au = File.binread(FILE_BIN, recl, @pos).unpack("d*")[0]
@pos += recl
rescue => e
raise
end
end
#=========================================================================
# EMRAT 取得
#
# * EMRAT(地球と月の質量比)の値をインスタンス変数 @emrat に設定する
# * 8 byte * 1
# * 倍精度浮動小数点数(機種依存)
# * 取得したバイト数だけレコード位置 @pos を進める
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_emrat
recl = 8
begin
@emrat = File.binread(FILE_BIN, recl, @pos).unpack("d*")[0]
@pos += recl
rescue => e
raise
end
end
#=========================================================================
# NUMDE 取得
#
# * DEバージョン番号(NUMDE)をインスタンス変数 @ver_de に設定する
# * 4 byte * 1
# * unsigned int (符号なし整数, エンディアンと int のサイズに依存)
# * 取得したバイト数だけレコード位置 @pos を進める
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_numde
recl = 4
begin
@ver_de = File.binread(FILE_BIN, recl, @pos).unpack("I*")[0]
@pos += recl
rescue => e
raise
end
end
#=========================================================================
# IPT 取得
#
# * IPT(オフセット、係数の数、サブ区間数)(水星〜地球の章動)の値をイ
# ンスタンス変数(配列) @ipts に設定する
# * 4 byte * 12 * 3
# * unsigned int (符号なし整数, エンディアンと int のサイズに依存)
# * 取得したバイト数だけレコード位置 @pos を進める
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_ipt
recl = 4
begin
@ipts = (0..11).map do |i|
ary = (0..2).map do |j|
File.binread(FILE_BIN, recl, @pos + recl * j).unpack("I*")[0]
end
@pos += recl * 3
ary
end
rescue => e
raise
end
end
#=========================================================================
# IPT_13(月の秤動)
#
# * IPT(オフセット、係数の数、サブ区間数)(月の秤動)の値をインスタン
# ス変数(配列) @ipts に設定する
# * 4 byte * 1 * 3
# * unsigned int (符号なし整数, エンディアンと int のサイズに依存)
# * 取得したバイト数だけレコード位置 @pos を進める
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_ipt_13
recl = 4
begin
@ipts << (0..2).map do |i|
File.binread(FILE_BIN, recl, @pos + recl * i).unpack("I*")[0]
end
@pos += recl * 3
rescue => e
raise
end
end
#=========================================================================
# CVAL 取得
#
# * レコード位置 @pos を2レコード目の先頭へ進める
# * 定数の値をインスタンス変数(配列) @cvals に設定する
# * 8 byte * @ncon
# * 倍精度浮動小数点数(機種依存)
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_cval
@pos = KSIZE * RECL
recl = 8
begin
@cvals = (0..@ncon - 1).map do |i|
File.binread(FILE_BIN, recl, @pos + recl * i).unpack("d*")[0]
end
rescue => e
raise
end
end
#=========================================================================
# JDEPOC 取得
#
# * @cvals(定数値配列)の中から「元期」をインスタンス変数 @jdepoc に設
# 定する
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_jdepoc
@jdepoc = @cvals[4]
rescue => e
raise
end
#=========================================================================
# 計算対象フラグ一覧取得
#
# * チェビシェフ多項式による計算が必要な天体の一覧をインスタンス変数(配
# 列) @list に設定する
# * 配列の並び順(係数データの並び順から「太陽」を除外した13個)
# @list = [
# 水星, 金星, 地球 - 月の重心, 火星, 木星, 土星, 天王星, 海王星,
# 冥王星, 月(地心), 地球の章動, 月の秤動
# ]
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_list
if @astrs[0] == 14
@list[10] = KIND if @ipts[11][1] > 0
return
end
if @astrs[0] == 15
@list[11] = KIND if @ipts[12][1] > 0
return
end
@astrs.each do |k|
@list[k - 1] = KIND if k <= 10
@list[2] = KIND if k == 10
@list[9] = KIND if k == 3
@list[2] = KIND if k == 13
end
rescue => e
raise
end
#=========================================================================
# COEFF 取得
#
# * レコード位置 @pos を3レコード目の先頭へ進める
# * 対象区間のユリウス日(開始、終了)をインスタンス変数(配列) @jds に
# 設定する
# * 係数の値をインスタンス変数(配列) @coeffs に設定する
# (x, y, z やサブ区間毎に分割して格納)
# * 8 byte * ?
# * 倍精度浮動小数点数(機種依存)
# * @coeffs に対象のインデックスの全ての係数を格納
# * 最初の2要素は当該データの開始・終了ユリウス日
#
# @params: <none>
# @return: <none>
#=========================================================================
def get_coeff
idx = ((@jd - @sss[0]) / @sss[2]).floor # レコードインデックス
@pos = KSIZE * RECL * (2 + idx)
recl = 8
begin
items = (0..(KSIZE / 2) - 1).map do |i|
File.binread(FILE_BIN, recl, @pos + recl * i).unpack("d*")[0]
end
@jds = [items.shift, items.shift]
@ipts.each_with_index do |ipt, i|
n = i == 11 ? 2 : 3 # 要素数
ary_1 = Array.new
ipt[2].times do |j|
ary_0 = Array.new
n.times do |k|
ary_0 << items.shift(ipt[1])
end
ary_1 << ary_0
end
@coeffs << ary_1
end
rescue => e
raise
end
end
#=========================================================================
# 補間
#
# * 使用するチェビシェフ多項式の係数は、
# * 天体番号が 1 〜 13 の場合は、 x, y, z の位置・速度(6要素)、
# 天体番号が 14 の場合は、 Δψ, Δε の角位置・角速度(4要素)、
# 天体番号が 15 の場合は、 φ, θ, ψ の角位置・角速度(6要素)。
# * 天体番号が 12 の場合は、 x, y, z の位置・速度の値は全て 0.0 とする。
#
# @params: astr(天体番号)
# @return: pvs = [
# x 位置, y 位置, z 位置,
# x 速度, y 速度, z 速度
# ]
# - 14(地球の章動)の場合は、
# pvs = [
# Δψ の角位置, Δε の角位置,
# Δψ の角速度, Δε の角速度
# ]
# - 15(月の秤動)の場合は、
# pvs = [
# φ の角位置, θ の角位置, ψ の角位置,
# φ の角速度, θ の角速度, ψ の角速度
# ]
#=========================================================================
def interpolate(astr)
pvs = Array.new
begin
tc, idx_sub = norm_time(astr)
n_item = astr == 14 ? 2 : 3 # 要素数
i_ipt = astr > 13 ? astr - 3 : astr - 1
i_coef = astr > 13 ? astr - 3 : astr - 1
# 位置
ary_p = [1, tc]
2.upto(@ipts[i_ipt][1] - 1) do |i|
ary_p << 2 * tc * ary_p[-1] - ary_p[-2]
end # 各項
0.upto(n_item - 1) do |i|
val = 0
0.upto(@ipts[i_ipt][1] - 1) do |j|
val += @coeffs[i_coef][idx_sub][i][j] * ary_p[j]
end
val /= @au if !KM && astr < 14
pvs << val
end #
# 速度
ary_v = [0, 1, 2 * 2 * tc]
3.upto(@ipts[i_ipt][1] - 1) do |i|
ary_v << 2 * tc * ary_v[-1] + 2 * ary_p[i - 1] - ary_v[-2]
end # 各項
0.upto(n_item - 1) do |i|
val = 0
0.upto(@ipts[i_ipt][1] - 1) do |j|
val += @coeffs[i_coef][idx_sub][i][j] * ary_v[j] * 2 * @ipts[i_ipt][2] / @sss[2].to_f
end
val /= KM ? 86400.0 : @au if astr < 14
pvs << val
end #
return pvs
rescue => e
raise
end
end
#=========================================================================
# チェビシェフ多項式用に時刻を正規化、サブ区間のインデックス算出
#
# @params: astr(天体番号)
# @return: [チェビシェフ時間, サブ区間のインデックス]
#=========================================================================
def norm_time(astr)
idx = astr > 13 ? astr - 2 : astr
jd_start = @jds[0]
tc = (@jd - jd_start) / @sss[2].to_f
temp = tc * @ipts[idx - 1][2]
idx = (temp - tc.floor).floor # サブ区間のインデックス
tc = (temp % 1.0 + tc.floor) * 2 - 1 # チェビシェフ時間
return [tc, idx]
rescue => e
raise
end
#=========================================================================
# 結果出力
#
# @params: <none>
# @return: <none>
#=========================================================================
def display
strs = " Target: %2d" % @astrs[0]
strs << " (#{ASTRS[@astrs[0] - 1]})\n"
strs << " Center: %2d" % @astrs[1]
strs << " (#{ASTRS[@astrs[1] - 1]})" unless @astrs[1] == 0
strs << "\n"
strs << " JD: %16.8f day\n" % @jd
strs << " 1AU: %11.1f km\n" % @au
puts strs
strs = ""
case @astrs[0]
when 14
strs << " Position(Δψ) = %32.20f rad\n" % @rrds[0]
strs << " Position(Δε) = %32.20f rad\n" % @rrds[1]
strs << " Velocity(Δψ) = %32.20f rad/day\n" % @rrds[2]
strs << " Velocity(Δε) = %32.20f rad/day\n" % @rrds[3]
when 15
strs << " Position(φ) = %32.20f rad\n" % @rrds[0]
strs << " Position(θ) = %32.20f rad\n" % @rrds[1]
strs << " Position(ψ) = %32.20f rad\n" % @rrds[2]
strs << " Velocity(φ) = %32.20f rad/day\n" % @rrds[3]
strs << " Velocity(θ) = %32.20f rad/day\n" % @rrds[4]
strs << " Velocity(ψ) = %32.20f rad/day\n" % @rrds[5]
else
unit_0, unit_1 = KM ? ["km", "sec"] : ["AU", "day"]
strs << " Position(x) = %32.20f #{unit_0}\n" % @rrds[0]
strs << " Position(y) = %32.20f #{unit_0}\n" % @rrds[1]
strs << " Position(z) = %32.20f #{unit_0}\n" % @rrds[2]
strs << " Velocity(x) = %32.20f #{unit_0}/#{unit_1}\n" % @rrds[3]
strs << " Velocity(y) = %32.20f #{unit_0}/#{unit_1}\n" % @rrds[4]
strs << " Velocity(z) = %32.20f #{unit_0}/#{unit_1}\n" % @rrds[5]
end
puts strs
rescue => e
raise
end
end
JplCalcDe430.new.exec if __FILE__ == $0
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.