1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
| #! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= JPLEPH(JPL の DE430 バイナリデータ)読み込み、座標（位置・速度）を計算する
# : 自作の RubyGems ライブラリ eph_jpl を使用する
#
# date name version
# 2016.06.19 mk-mode.com 1.00 新規作成
#
# Copyright(C) 2016 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# 【引数】
# [第１] 対象天体番号（必須）
# 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)
# [第２] 基準天体番号（必須。 0, 1 - 13）
# （ 0 は、対象天体番号が 14, 15 のときのみ）
# [第３] ユリウス日（省略可。省略時は現在日時のユリウス日を地球時とみなす）
#---------------------------------------------------------------------------------
# 【注意事項】
# * 求める座標は「赤道直角座標(ICRS)」
# * 時刻系は「太陽系力学時(TDB)」である。（≒地球時(TT)）
# * 天体番号が 1 〜 13 の場合は、 x, y, z の位置・速度（6要素）、
# 天体番号が 14 の場合は、黄経における章動 Δψ, 黄道傾斜における章動 Δε の
# 角位置・角速度（4要素）、
# 天体番号が 15 の場合は、 φ, θ, ψ の角位置・角速度（6要素）。
# * 対象天体番号 = 基準天体番号 は、無意味なので処理しない。
#---------------------------------------------------------------------------------
#++
require 'eph_jpl'
class EphemerisJpl
FILE_BIN = "/home/masaru/src/ephemeris_jpl/JPLEPH"
USAGE = <<-"EOS"
【使用方法】 ./ephemeris_jpl.rb 対象天体番号 基準天体番号 [ユリウス日]
【天体番号】（対象天体番号: 1 - 15, 基準天体番号: 0 - 13）
1: 水星, 2: 金星, 3: 地球, 4: 火星, 5: 木星,
6: 土星, 7: 天王星, 8: 海王星, 9: 冥王星, 10: 月,
11: 太陽, 12: 太陽系重心, 13: 地球 - 月の重心,
14: 地球の章動, 15: 月の秤動,
0: 対象天体番号 14, 15 の場合に指定
※対象天体番号≠基準天体番号
EOS
KM = false # true: km, false: AU
def initialize
target, center, jd = get_args
@eph = EphJpl.new(FILE_BIN, target, center, jd, KM)
end
def exec
display(@eph.calc)
rescue => e
msg = "[#{e.class}] #{e.message}\n"
msg << e.backtrace.each { |tr| "\t#{tr}"}.join("\n")
$stderr.puts msg
exit 1
end
private
#=========================================================================
# コマンドライン引数取得
#
# @param: <none>
# @return: [target, center, jd]
#=========================================================================
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
# ユリウス日
jd = ARGV[2]
if jd
jd = jd.to_f
else
now = Time.now
jd = gc2jd({
year: now.year, month: now.month, day: now.day,
hour: now.hour, min: now.min, sec: now.sec
})
end
return [target, center, jd]
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 gc2jd(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.5
# 時間(小数)部分計算
t = sec / 3600.0
t += min / 60.0
t += hour
t /= 24.0
return jd + t
rescue => e
raise
end
end
#=========================================================================
# 結果出力
#
# @param: pvs (= Position-Velocity array)
# @return: <none>
#=========================================================================
def display(pvs)
strs = sprintf(" Target: %2d", @eph.target)
strs << " (#{@eph.target_name})\n"
strs << sprintf(" Center: %2d", @eph.center)
strs << " (#{@eph.center_name})" unless @eph.center == 0
strs << "\n"
strs << sprintf(" JD: %16.8f day\n", @eph.jd)
strs << sprintf(" 1AU: %11.1f km\n", @eph.bin[:au])
case @eph.target
when 14
strs << sprintf(" Position(Δψ) = %30.20f rad\n", pvs[0])
strs << sprintf(" Position(Δε) = %30.20f rad\n", pvs[1])
strs << sprintf(" Velocity(Δψ) = %30.20f rad/day\n", pvs[2])
strs << sprintf(" Velocity(Δε) = %30.20f rad/day\n", pvs[3])
when 15
strs << sprintf(" Position(φ) = %30.20f rad\n", pvs[0])
strs << sprintf(" Position(θ) = %30.20f rad\n", pvs[1])
strs << sprintf(" Position(ψ) = %30.20f rad\n", pvs[2])
strs << sprintf(" Velocity(φ) = %30.20f rad/day\n", pvs[3])
strs << sprintf(" Velocity(θ) = %30.20f rad/day\n", pvs[4])
strs << sprintf(" Velocity(ψ) = %30.20f rad/day\n", pvs[5])
else
unit_p, unit_v = @eph.unit.split(", ")
strs << sprintf(" Position(x) = %30.20f #{unit_p}\n", pvs[0])
strs << sprintf(" Position(y) = %30.20f #{unit_p}\n", pvs[1])
strs << sprintf(" Position(z) = %30.20f #{unit_p}\n", pvs[2])
strs << sprintf(" Velocity(x) = %30.20f #{unit_v}\n", pvs[3])
strs << sprintf(" Velocity(y) = %30.20f #{unit_v}\n", pvs[4])
strs << sprintf(" Velocity(z) = %30.20f #{unit_v}\n", pvs[5])
end
puts strs
rescue => e
raise
end
end
EphemerisJpl.new.exec if __FILE__ == $0
|