Ruby - JPL 天文暦データから ICRS 座標を計算!

Updated:


前回、NASA の機関である JPL(Jet Propulsion Laboratory) が惑星探査用に編纂・発行している月・惑星の暦の最新版 DE430 のバイナリ形式のデータを Ruby で読み込みました。

今回は、読み込んだデータから ICRS 座標を計算してみました。(読み込んだデータとは、 ICRS 座標の計算に必要なチェビシェフ多項式の係数データ)

0. 前提条件、注意事項

  • Ruby 2.3.0-p0 での作業を想定。
  • 使用するバイナリ形式データは、テキスト形式データを「JPL 天文暦データのバイナリ化!」の方法でバイナリ化したもの。(今回は DE430 用の “ascp1550.430” 〜 “ascp2550.430” 全てをバイナリ化。ファイル名は “JPLEPH” に変更)
  • あらかじめ用意されているバイナリデータを使用する場合は KSIZE 値が組み込まれていないので、別途テキスト形式のヘッダファイルを読み込む必要がある。
  • テキスト形式データの仕様については「月・惑星の暦 JPL DE430 について!」を参照。
  • バイナリ形式データの仕様については「JPL 天文暦バイナリデータの仕様!」を参照。

1. Ruby スクリプトの作成

(バイナリデータ読み込みには IO.binread メソッドを使用)

File: jpl_calc_de430

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
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
#! /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

(ソーススクリプトがより実用に耐えられるものに洗練されれば、Gist でなくリポジトリとしても登録するつもり)

2. Ruby スクリプトの実行

バイナリ形式データ “JPLEPH” を Ruby スクリプトと同じディレクトリに配置して以下のように実行する。(以下は、「2016年3月18日0時0分0秒」の「地球」を基準とした「太陽」の ICRS 座標を計算する例)

$ ./jpl_calc_de430.rb 11 3 2457465.5
  Target: 11 (Sun)
  Center:  3 (Earth)
      JD: 2457465.50000000 day
     1AU: 149597870.7 km
  Position(x) =         0.99443659220700997281 AU
  Position(y) =        -0.03816291768957833647 AU
  Position(z) =        -0.01655177670960058384 AU
  Velocity(x) =         0.00099333756561238804 AU/day
  Velocity(y) =         0.01582779844821733881 AU/day
  Velocity(z) =         0.00686186627679565356 AU/day

3. データの検証

“testpo.430” のデータと突き合わせながら Ruby スクリプトを実行して検証してみる。

実際には、”testpo.430” に記述されている任意の行に注目し、その行に記載の2つの天体の ICRS 座標を計算して差を求め、その値が coordinates の値と等しいかを確認する。(但し、”testpo.430” と天体番号が異なるものがあるため、適宜置き換えて考えること)

例えば、 “testpo.430” の 2016.12.01 の行は、以下のようになっている。(ヘッダと当該行のみ抜粋)

File: tespo.430

1
2
3
de#  -- date -- -- jed -- t# c# x# -- coordinate ---

430  2016.12.01 2457723.5  8  5  2       -8.02225801549924900000

ユリウス日 2457723.5 の天体番号 5(木星)を基準とした天体番号 8(海王星)の ICRS 座標を計算すると次のようになる。

$ ./jpl_calc_de430.rb 8 5 2457723.5
  Target:  8 (Neptune)
  Center:  5 (Jupiter)
      JD: 2457723.50000000 day
     1AU: 149597870.7 km
  Position(x) =        33.70070119505614059108 AU
  Position(y) =        -8.02225801549925066070 AU
  Position(z) =        -4.10459769835983934172 AU
  Velocity(x) =        -0.00000061831278598047 AU/Day
  Velocity(y) =         0.00929863072586037051 AU/Day
  Velocity(z) =         0.00393151443263963212 AU/Day

“testpo.430” のインデックス番号(x# 列の値)が 2 であるので、coordinate 列の値と算出結果の 2 番目の値(Position(y) の値)を比較してみる。
-8.02225801549924900000-8.02225801549925066070 であるので、かなりの精度で一致している。

他についても同様に確認してみる。(特殊なパターンもあるので要注意。JPL から提供されている “testeph.f” や上記のソーススクリプトを熟知していれば分かるはず)

“testpo.430” の仕様(ヘッダファイルやデータファイルと検証用ファイルでの天体番号の相違等)については「月・惑星の暦 JPL DE430 について!」で説明しているので、そちらを参照のこと。

当方、 “testpo.430” の任意の行について検証してみた結果、かなりの精度差異がないことを確認済み。(小数点以下約13桁の精度。それ以下の誤差は(丸め)誤差の範囲内という認識)

また、地球を基準とした太陽の位置は DE430 データを使用している日本の国立天文台の計算結果とも一致しました。(「太陽の赤道直角座標と黄経・黄緯」で計算)


バイナリデータから目的の惑星・日付の「ICRS 座標系(国際天文基準座標系、原点が太陽系重心の座標系)」の座標が計算できるようになったので、次は「赤道座標(地球の自転軸によって定まる座標)」への変換や「黄道座標(地球の公転軌道によって定まる座標)」の計算にも挑戦してみたいです。

以上。





 

Sponsored Link

 

Comments