mk-mode BLOG

このブログは自作の自宅サーバに構築した Debian GNU/Linux で運用しています。
PC・サーバ構築等の話題を中心に公開しております。(クローンサイト: GitHub Pages

ブログ開設日2009-01-05
サーバ連続稼働時間
Reading...
Page View 合計
Reading...
今日
Reading...
昨日
Reading...

Ruby - JPL DE430 データから太陽・月の視位置を計算!

[ プログラミング, 暦・カレンダー ] [ Ruby, カレンダー ]

こんばんは。

NASA の機関である JPL(Jet Propulsion Laboratory) が惑星探査用に編纂・発行している太陽・月・惑星の暦の最新版 DE430 からデータを取得し、太陽と月の視位置を高精度で計算してみました。

0. 前提条件

  • Ruby 2.3.1-p112 での作業を想定。
  • 以下の自作 RubyGems ライブラリをインストールしておく。
    • ICRS 座標取得には eph_jpl を使用する。
    • バイアス・歳差・章動適用には eph_bpn を使用する。
    • 時刻系変換には mk_time を使用する。
    • 座標系変換には mk_coord を使用する。

1. 事前準備

JPL DE430 のデータを使用するので、バイナリデータ “linux_p1550p2650.430” をこちらから取得し “JPLEPH” とリネームして適当なディレクトリに配置しておく。

ただし、このバイナリファイルはサイズが大きいので、必要な年代のテキストデータのみをマージ&バイナリ化してもよい。(参考:JPL 天文暦データのバイナリ化!

2. Ruby スクリプトの作成

  • FILE_BIN の値は、 “JPLEPH” の配置場所に合わせて編集すること。
  • 重力の大きい太陽・木星・土星・天王星・海王星の重力場による光の曲がりは考慮していない。(太陽や月の視位置計算には影響がないと判断したので)
apparent_sun_moon_jpl.rb
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
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= 太陽・月の視位置計算
#
#  * JPLEPH(JPL の DE430 バイナリデータ)読み込み、視位置を計算する
#    - ICRS 座標取得には自作 RubyGems ライブラリ eph_jpl を使用する
#    - バイアス・歳差・章動適用には自作 RubyGems ライブラリ eph_bpn を使用する
#    - 時刻系変換には自作 RubyGems ライブラリ mk_time を使用する
#    - 座標系変換には自作 RubyGems ライブラリ mk_coord を使用する
#
# date          name            version
# 2016.08.23    mk-mode.com     1.00 新規作成
#
# Copyright(C) 2016 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# 引数 : 日時(JST)
#          書式:YYYYMMDD or YYYYMMDDHHMMSS
#          無指定なら現在(システム日時)と判断。
#---------------------------------------------------------------------------------
#++
require 'date'
require 'eph_bpn'   # 自作ライブラリ
require 'eph_jpl'   # 自作ライブラリ
require 'mk_coord'  # 自作ライブラリ
require 'mk_time'   # 自作ライブラリ

class ApparentSunMoon
  USAGE      = "[USAGE] ./apparent_sun_moon.rb [YYYYMMDD|YYYYMMDDHHMMSS]"
  MSG_ERR_1  = "[ERROR] Newton method error!"
  JST_OFFSET = 9             # JST - UTC
  FILE_BIN   = "JPLEPH"      # JPL バイナリデータ
  AU         = 149597870700  # 1天文単位 (m)
  C          = 299792458     # 光速 (m/s)
  DAYSEC     = 86400         # 1日の秒数
  BODIES     = {earth: 3, moon: 10, sun: 11}  # 天体名と JPL での天体番号
  DEBUG      = false

  #=========================================================================
  # 初期処理
  #
  # @param:  <none>
  # @return: <none>
  #=========================================================================
  def initialize
    # === t1(= TDB), t2(= TDB) における位置・速度(ICRS 座標)用Hash
    @icrs_1, @icrs_2 = Hash.new, Hash.new
    # === 視位置を求めたい時刻 t2 (= JST; 日本標準時) の決定
    @t_2_jst = get_arg
    # === 時刻 t2 の変換(JST(日本標準時) -> UTC(協定世界時))
    @t_2_utc = jst2utc(@t_2_jst)
    # === 時刻 t2 の変換(UTC(協定世界時) -> TDB(太陽系力学時))
    @t_2 = utc2tdb(@t_2_utc)
    # === 時刻 t2 のユリウス日
    @jd_2 = get_jd(@t_2)
    # === 時刻 t2(= TDB) におけるの位置・速度(ICRS 座標)の計算 (地球, 月, 太陽)
    BODIES.each { |k, v| @icrs_2[k] = get_icrs(v, @jd_2) }
    # === 時刻 t2(= TDB) における地球と太陽・月の距離
    @r_e = get_r_e
  end

  #=========================================================================
  # 主処理
  #
  # @param:  <none>
  # @return: <none>
  #=========================================================================
  def exec
    disp_init
    @eq_pol_s, @ec_pol_s = calc_sun
    @eq_pol_m, @ec_pol_m = calc_moon
    disp_result
  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: <none>
  #=========================================================================
  def calc_sun
    # === 太陽が光を発した時刻 t1(JD) の計算
    t_1_jd = calc_t1(:sun, @t_2)
    # === 時刻 t1(= TDB) におけるの位置・速度(ICRS 座標)の計算 (地球, 月, 太陽)
    BODIES.each { |k, v| @icrs_1[k] = get_icrs(v, t_1_jd) }
    # === 時刻 t2 における地球重心から時刻 t1 における太陽への方向ベクトルの計算
    v_12 = calc_unit_vector(@icrs_2[:earth][0,3], @icrs_1[:sun][0,3])
    # === GCRS 座標系: 光行差の補正(方向ベクトルの Lorentz 変換)
    dd = conv_lorentz(v_12)
    pos_sun = calc_pos(dd, @r_e[:sun])
    # === 瞬時の真座標系: GCRS への bias & precession(歳差) & nutation(章動)の適用
    bpn = EphBpn.new(@t_2.strftime("%Y%m%d%H%M%S"))
    pos_sun_bpn = bpn.apply_bpn(pos_sun)
if DEBUG
    puts "---"
    puts "* ICRS@t2 = #{@icrs_2}"
    puts "* 時刻 t2 における地球と太陽・月の距離 r_e\n  = #{@r_e}"
    puts "* 太陽が光を発した時刻 t1(JD)\n  = #{t_1_jd}"
    puts "* ICRS@t1\n  = #{@icrs_1}"
    puts "* 時刻 t2 における地球重心から時刻 t1 における太陽への方向ベクトル v_12\n  = #{v_12}"
    puts "* 光行差補正後(方向ベクトルの Lorentz 変換) dd\n  = #{dd}"
    puts "* 太陽位置(GCRS 座標系)"
    puts "  = #{pos_sun}"
    puts "* 太陽位置(瞬時の真座標系)"
    puts "  = #{pos_sun_bpn}"
end
    # === 座標変換
    eq_pol_s  = MkCoord.rect2pol(pos_sun_bpn)
    ec_rect_s = MkCoord.rect_eq2ec(pos_sun_bpn, bpn.eps)
    ec_pol_s  = MkCoord.rect2pol(ec_rect_s)
    return [eq_pol_s, ec_pol_s]
  rescue => e
    raise
  end

  #=========================================================================
  # 視位置計算(月)
  #
  # @param:  <none>
  # @return: <none>
  #=========================================================================
  def calc_moon
    # === 月が光を発した時刻 t1(jd) の計算
    t_1_jd = calc_t1(:moon, @t_2)
    # === 時刻 t1(= TDB) におけるの位置・速度(ICRS 座標)の計算 (地球, 月, 太陽)
    BODIES.each { |k, v| @icrs_1[k] = get_icrs(v, t_1_jd) }
    # === 時刻 t2 における地球重心から時刻 t1 における月への方向ベクトルの計算
    v_12 = calc_unit_vector(@icrs_2[:earth][0,3], @icrs_1[:moon][0,3])
    # === GCRS 座標系: 光行差の補正(方向ベクトルの Lorentz 変換)
    dd = conv_lorentz(v_12)
    pos_moon = calc_pos(dd, @r_e[:moon])
    # === 瞬時の真座標系: GCRS への bias & precession(歳差) & nutation(章動)の適用
    bpn = EphBpn.new(@t_2.strftime("%Y%m%d%H%M%S"))
    pos_moon_bpn = bpn.apply_bpn(pos_moon)
if DEBUG
    puts "---"
    puts "* ICRS@t2 = #{@icrs_2}"
    puts "* 時刻 t2 における地球と太陽・月との距離 r_e\n  = #{@r_e}"
    puts "* 月が光を発した時刻 t1(JD)\n  = #{t_1_jd}"
    puts "* ICRS@t1\n  = #{@icrs_1}"
    puts "* 時刻 t2 における地球重心から時刻 t1 における月への方向ベクトル v_12\n  = #{v_12}"
    puts "* 光行差補正後(方向ベクトルの Lorentz 変換) dd\n  = #{dd}"
    puts "* 月位置(GCRS 座標系)"
    puts "  = #{pos_moon}"
    puts "* 月位置(瞬時の真座標系)"
    puts "  = #{pos_moon_bpn}"
end
    # === 座標変換
    eq_pol_m  = MkCoord.rect2pol(pos_moon_bpn)
    ec_rect_m = MkCoord.rect_eq2ec(pos_moon_bpn, bpn.eps)
    ec_pol_m  = MkCoord.rect2pol(ec_rect_m)
    return [eq_pol_m, ec_pol_m]
  rescue => e
    raise
  end

  #=========================================================================
  # ICRS 座標取得
  #
  # * JPL DE430 データを自作 RubyGems ライブラリ eph_jpl を使用して取得。
  #
  # @param:  target  (対象天体(Symbol))
  # @param:  jd      (ユリウス日(Numeric))
  # @return: [
  #   pos_x, pos_y, pos_z,
  #   vel_x, vel_y, vel_z
  # ]                (位置・速度(Array), 単位: AU, AU/day)
  #=========================================================================
  def get_icrs(target, jd)
    return EphJpl.new(FILE_BIN, target, 12, jd).calc
  rescue => e
    raise
  end

  #=========================================================================
  # 引数取得
  #
  # @param:  <none>
  # @return: JST(Time Object)
  #=========================================================================
  def get_arg
    unless arg = ARGV.shift
      t = Time.now
      return Time.new(t.year, t.month, t.day, t.hour, t.min, t.sec)
    end
    (puts USAGE; exit 0) unless arg =~ /^\d{8}$|^\d{14}$/
    year, month, day = arg[ 0, 4].to_i, arg[ 4, 2].to_i, arg[ 6, 2].to_i
    hour, min,   sec = arg[ 8, 2].to_i, arg[10, 2].to_i, arg[12, 2].to_i
    (puts USAGE; exit 0) unless Date.valid_date?(year, month, day)
    (puts USAGE; exit 0) if hour > 23 || min > 59 || sec > 59
    return Time.new(year, month, day, hour, min, sec)
  rescue => e
    raise
  end

  #=========================================================================
  # JST -> UTC
  #
  # @param:  jst  (Time Object)
  # @return: utc  (Time Object)
  #=========================================================================
  def jst2utc(jst)
    return jst - JST_OFFSET * 60 * 60
  rescue => e
    raise
  end

  #=========================================================================
  # UTC -> TDB
  #
  # @param:  utc  (Time Object)
  # @return: tdb  (Time Object)
  #=========================================================================
  def utc2tdb(utc)
    return MkTime.new(utc.strftime("%Y%m%d%H%M%S")).tdb
  rescue => e
    raise
  end

  #=========================================================================
  # ユリウス日取得
  #
  # @param:  t   (Time Object)
  # @return: jd  (Julian Day(Numeric))
  #=========================================================================
  def get_jd(t)
    return MkTime.new(t.strftime("%Y%m%d%H%M%S")).jd
  rescue => e
    raise
  end

  #=========================================================================
  # t2(= TDB) における地球と太陽・月の距離
  #
  # @param:  <none>
  # @return: r_e  (地球と太陽・月の距離(Hash))
  #=========================================================================
  def get_r_e
    r_e = Hash.new

    begin
      @icrs_2.each do |k, v|
        next if k == :earth
        r_e[k] = calc_dist(@icrs_2[:earth][0, 3], v[0, 3])
      end
      return r_e
    rescue => e
      raise
    end
  end

  #=========================================================================
  # 天体Aを中心とした場合の天体Bの位置計算
  #
  # @param:  pos_a  (位置ベクトル(天体A))
  # @param:  pos_b  (位置ベクトル(天体B))
  # @return: pos    (位置ベクトル)
  #=========================================================================
  def calc_pos_ab(pos_a, pos_b)
    return (0..2).map { |i| pos_b[i] - pos_a[i] }
  rescue => e
    raise
  end

  #=========================================================================
  # 天体Aと天体Bの距離計算
  #
  # @param:   pos_a  (位置ベクトル)
  # @param:   pos_a  (位置ベクトル)
  # @return:  r      (距離)
  #=========================================================================
  def calc_dist(pos_a, pos_b)
    r =  (0..2).inject(0.0) { |sum, i| sum + (pos_b[i] -pos_a[i]) ** 2 }
    return Math.sqrt(r)
  rescue => e
    raise
  end

  #=========================================================================
  # 対象天体が光を発した時刻 t1 の計算(太陽・月専用)
  #
  # * 計算式: c * (t2 - t1) = r12  (但し、 c: 光の速度。 Newton 法で近似)
  # * 太陽・月専用なので、太陽・木星・土星・天王星・海王星の重力場による
  #   光の曲がりは非考慮。
  #
  # @param:  target  (対象天体(Symbol))
  # @param:  tdb     (Time Object(観測時刻;TDB))
  # @return: t_1     (Time Ojbect)
  #=========================================================================
  def calc_t1(target, tdb)
    t_1 = MkTime.new(tdb.strftime("%Y%m%d%H%M%S")).jd
    t_2 = t_1
    pv_1 = @icrs_2[target]
    df, m = 1.0, 0
    while df > 1.0e-10
      r_12 = (0..2).map { |i| pv_1[i] - @icrs_2[:earth][i] }
      r_12_d = calc_dist(pv_1, @icrs_2[:earth])
      df = (C * DAYSEC / AU) * (t_2 - t_1) - r_12_d
      df_wk = (0..2).inject(0.0) { |s, i| s + r_12[i] * @icrs_2[target][i + 3] }
      df /= ((C * DAYSEC / AU) + (df_wk) / r_12_d)
      t_1 += df
      m += 1
      raise MSG_ERR_1 if m > 10
      pv_1 = get_icrs(BODIES[target], t_1)
    end
    return t_1
  rescue => e
    raise
  end

  #=========================================================================
  # 天体Aから見た天体Bの方向ベクトル計算(太陽・月専用)
  #
  # * 太陽・月専用なので、太陽・木星・土星・天王星・海王星の重力場による
  #   光の曲がりは非考慮。
  #
  # @param:   pos_a  (位置ベクトル(天体A))
  # @param:   pos_b  (位置ベクトル(天体B))
  # @return:  vec    (方向(単位)ベクトル)
  #=========================================================================
  def calc_unit_vector(pos_a, pos_b)
    vec = [0.0, 0.0, 0.0]

    begin
      w = calc_dist(pos_a, pos_b)
      vec = (0..2).map { |i| pos_b[i] - pos_a[i] }
      return vec.map { |v| v / w } unless w == 0.0
    rescue => e
      raise
    end
  end

  #=========================================================================
  # 光行差の補正(方向ベクトルの Lorentz 変換)
  #
  # * vec_dd = f * vec_d + (1 + g / (1 + f)) * vec_v
  #   但し、 f = vec_v * vec_d  (ベクトル内積)
  #          g = sqrt(1 - v^2)  (v: 速度)
  #
  # @param:  v_d   (方向(単位)ベクトル)
  # @return: v_dd  (補正後ベクトル)
  #=========================================================================
  def conv_lorentz(vec_d)
    vec_v = @icrs_2[:earth][3,3].map { |v| (v / DAYSEC) / (C / AU.to_f) }
    g = inner_prod(vec_v, vec_d)
    f = Math.sqrt(1.0 - calc_velocity(vec_v))
    vec_dd_1 = vec_d.map { |d| d * f }
    vec_dd_2 = vec_v.map { |v| (1.0 + g / (1.0 + f)) * v }
    vec_dd = (0..2).map { |i| vec_dd_1[i] + vec_dd_2[i] }.map { |a| a / (1.0 + g) }
    return vec_dd
  rescue => e
    raise
  end

  #=========================================================================
  # 単位(方向)ベクトルと距離から位置ベクトルを計算
  #
  # @param:  d    (単位(方向)ベクトル)
  # @param:  r    (距離)
  # @return  pos  (位置ベクトル)
  #=========================================================================
  def calc_pos(d, r)
    return d.map { |a| a * r }
  rescue => e
    raise
  end

  #=========================================================================
  # ベクトルの外積
  #
  # @param:  a  (ベクトル)
  # @param:  b  (ベクトル)
  # @return: w  (ベクトル)
  #=========================================================================
  def cross_prod(a, b)
    return [
      a[1] * b[2] - a[2] * b[1],
      a[2] * b[0] - a[0] * b[2],
      a[0] * b[1] - a[1] * b[0]
    ]
  rescue => e
    raise
  end

  #=========================================================================
  # ベクトルの内積
  #
  # @param:  a  (ベクトル)
  # @param:  b  (ベクトル)
  # @return: w  (スカラー(内積の値))
  #=========================================================================
  def inner_prod(a, b)
    return (0..2).inject(0.0) { |s, i| s + a[i] * b[i] }
  rescue => e
    raise
  end

  #=========================================================================
  # 天体の速度ベクトルから実際の速度を計算
  #
  # @param:   vec  (ベクトル)
  # @return:  v    (速度)
  #=========================================================================
  def calc_velocity(vec)
    v = vec.inject(0.0) { |s, i| s + i * i }
    return Math.sqrt(v)
  rescue => e
    raise
  end

  #=========================================================================
  # 結果出力(時刻パラメータ)
  #
  # @param:  <none>
  # @return: <none>
  #=========================================================================
  def disp_init
    print "* 計算対象時刻 t2(JST) = "
    puts @t_2_jst.instance_eval { '%s.%03d' % [strftime('%Y-%m-%d %H:%M:%S'), (usec / 1000.0).round] }
    print "                 (UTC) = "
    puts @t_2_utc.instance_eval { '%s.%03d' % [strftime('%Y-%m-%d %H:%M:%S'), (usec / 1000.0).round] }
    print "                 (TDB) = "
    puts @t_2.instance_eval { '%s.%03d' % [strftime('%Y-%m-%d %H:%M:%S'), (usec / 1000.0).round] }
    puts "                  (JD) = #{@jd_2}"
  rescue => e
    raise
  end

  #=========================================================================
  # 結果出力(計算結果)
  #
  # param:  <none>
  # return: <none>
  #=========================================================================
  def disp_result
    puts "* 視位置(太陽)"
    printf("  = [赤経: %14.10f rad, 赤緯: %14.10f rad]\n", @eq_pol_s[0], @eq_pol_s[1])
    printf("  = [赤経: %14.10f deg, 赤緯: %14.10f deg]\n", @eq_pol_s[0] * 180 / PI, @eq_pol_s[1] * 180 / PI)
    printf("  = [黄経: %14.10f rad, 黄緯: %14.10f rad]\n", @ec_pol_s[0], @ec_pol_s[1])
    printf("  = [黄経: %14.10f deg, 黄緯: %14.10f deg]\n", @ec_pol_s[0] * 180 / PI, @ec_pol_s[1] * 180 / PI)
    puts "* 視位置(月)"
    printf("  = [赤経: %14.10f rad, 赤緯: %14.10f rad]\n", @eq_pol_m[0], @eq_pol_m[1])
    printf("  = [赤経: %14.10f deg, 赤緯: %14.10f deg]\n", @eq_pol_m[0] * 180 / PI, @eq_pol_m[1] * 180 / PI)
    printf("  = [黄経: %14.10f rad, 黄緯: %14.10f rad]\n", @ec_pol_m[0], @ec_pol_m[1])
    printf("  = [黄経: %14.10f deg, 黄緯: %14.10f deg]\n", @ec_pol_m[0] * 180 / PI, @ec_pol_m[1] * 180 / PI)
    puts "* 視黄経差(太陽 - 月)"
    printf("  = %15.10f rad\n", @ec_pol_s[0] - @ec_pol_m[0])
    printf("  = %15.10f deg\n", @ec_pol_s[0] * 180 / PI - @ec_pol_m[0] * 180 / PI)
    puts "* 距離(太陽)"
    printf("  = %.10f AU = %.2f km\n", @ec_pol_s[2], @ec_pol_s[2] * AU / 1000.0)
    puts "* 距離(月)"
    printf("  = %.10f AU = %.2f km\n", @ec_pol_m[2], @ec_pol_m[2] * AU / 1000.0)
  rescue => e
    raise
  end
end

exit 0 unless __FILE__ == $0
ApparentSunMoon.new.exec

(いつもなら GitHub Gist にソーススクリプトを掲載するが、今回のスクリプトは近いうちに RubyGems ライブラリ化する予定なので Gist には掲載しない)

3. Ruby スクリプトの実行

コマンドライン引数に日時(JST)を YYYYMMDD or YYYYMMDDHHMMSS の書式で指定指定実行する。(指定しない場合はシステム時刻を JST とみなす)
以下は 2017年02月26日23時58分22秒(視黄経差が 0.0 に近づく日時) の例。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
$ ./apparent_sun_moon_jpl.rb 20170226235822
* 計算対象時刻 t2(JST) = 2017-02-26 23:58:22.000
                 (UTC) = 2017-02-26 14:58:22.000
                 (TDB) = 2017-02-26 14:59:31.184
                  (JD) = 2457811.124664352
* 視位置(太陽)
  = [赤経:   5.9314958658 rad, 赤緯:  -0.1482278532 rad]
  = [赤経: 339.8496793121 deg, 赤緯:  -8.4928303972 deg]
  = [黄経:   5.9027352417 rad, 黄緯:   0.0000119199 rad]
  = [黄経: 338.2018169337 deg, 黄緯:   0.0006829627 deg]
* 視位置(月)
  = [赤経:   5.9344195445 rad, 赤緯:  -0.1554065421 rad]
  = [赤経: 340.0171937614 deg, 赤緯:  -8.9041389710 deg]
  = [黄経:   5.9027340610 rad, 黄緯:  -0.0077266797 rad]
  = [黄経: 338.2017492834 deg, 黄緯:  -0.4427061343 deg]
* 視黄経差(太陽 - 月)
  =    0.0000011807 rad
  =    0.0000676503 deg
* 距離(太陽)
  = 0.9902304257 AU = 148136363.19 km
* 距離(月)
  = 0.0025279553 AU = 378176.74 km

4. 計算結果の検証

太陽の視位置については「太陽の地心座標 - 国立天文台暦計算室」、月の視位置については「月の地心座標 - 国立天文台暦計算室」で計算した結果と比較してみる。

今回の例の場合、太陽の視視黄緯以外は 1/1000〜1/100000 の精度で、太陽の視黄緯は 1/100 の精度で一致した。(単位:「度」)

そして、何よりも、視黄経差(太陽 - 月)が「2017年02月26日23時58分」に 0.0 に最も近づいたことで、精度が高いと判断できる。

なぜなら、視黄経差が 0.0 になる瞬間は「旧暦」を生成する際に非常に重要になるからで、「2017年の2月に訪れる旧暦の1日は、26日か27日か?」という問題も出てくるほど。
実際、国立天文台の計算では、その頃の新月(朔)は 26日23時58分とされており、旧暦の1日は「26日」である。
(精度の低いアルゴリズムで計算すると、新月の瞬間が「27日」になってしまう)

5. 参考サイト


これで、精度の高くないアルゴリズムで計算していた際に発生していた旧暦がズレて都度調整していた問題が解消されます。

近いうちに、今回のスクリプト自体も RubyGems ライブラリ化し、カレンダー計算に利用したいと考えております。

※当記事「執筆」時点では RubyGems ラリブラリ化は実現していませんでしたが、当記事「公開」時点では実現しています。次回、公開する予定。

以上。

Comments