mk-mode BLOG

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

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

Ruby - 章動の計算(IAU2000A 理論)!

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

こんばんは。

天体位置の回転に使用する章動を Ruby で実装して計算してみました。(使用するのは IAU2000A 理論)

歳差・章動の変換については、過去記事をご参照ください。

0. 前提条件、事前知識

  • 章動の計算には、 IAU SOFA(International Astronomical Union, Standards of Fundamental Astronomy) の提供する C ソースコードに実装されているアルゴリズム “nut00a.c” を使用する。
  • IAU SOFA のソースコードには、 MHB2000(Mathews-Herring-Buffett, 2000) の理論や IERS2003(International Earth Rotation & Reference Systems service, 2003) の理論の使用が混在していることに留意。
  • ここでは「章動(しょうどう)」そのものが何かについては詳細には説明しないが、簡単に説明すると、章動には
    • 黄経における章動()
    • 黄道傾斜における章動()

    があり、それぞれが

    • 日月章動(luni-solar nutation)
    • 惑星章動(planetary nutation)

    で構成されている。

  • また、算出アルゴリズムについてもここでは詳細には説明しない。(と言うより、煩雑で自分には説明できない)
    参考サイトやソーススクリプトを参照のこと。
  • 今回は IAU2000A 理論に基づいて計算しているが、簡略版の IAU2000B 理論 “nut00b.c” も存在する。
  • IAU2000A 理論に若干の補正を加えた IAU2006 理論も存在するが、今回は非考慮。

1. 係数データの準備

計算に使用する係数データを用意する。

IERS Conventions Center のページから、日月章動用の係数データ Lunisolar components と惑星章動用の係数データ Planetary components を取得して、扱いやすいように整形する。

整形したものは以下にアップしている。

2. Ruby スクリプトの作成

IAU SOFA の提供する C ソースコードに実装されているアルゴリズム “nut00a.c” に則って作成している。

nutation_model.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
#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= 章動の計算
#  : IAU2000A 章動理論(MHB2000, IERS2003)による
#    黄経における章動(Δψ), 黄道傾斜における章動(Δε) の計算
#
#  * IAU SOFA(International Astronomical Union, Standards of Fundamental Astronomy)
#    の提供する C ソースコード "nut00a.c" で実装されているアルゴリズムを使用する。
#  * 係数データファイルの項目について
#    - 日月章動(luni-solar nutation, "dat_ls.txt")
#      (左から) L L' F D Om PS PST PC EC ECT ES
#    - 惑星章動(planetary nutation, "dat_pl.txt)
#      (左から) L L' F D Om Lm Lv Le LM Lj Ls Lu Ln Pa PS PC ES EC
#  * 参考サイト
#    - [SOFA Library Issue 2012-03-01 for ANSI C: Complete List](http://www.iausofa.org/2012_0301_C/sofa/)
#    - [USNO Circular 179](http://aa.usno.navy.mil/publications/docs/Circular_179.php)
#    - [IERS Conventions Center](http://62.161.69.131/iers/conv2003/conv2003_c5.html)
#
# Date          Author          Version
# 2016.05.26    mk-mode.com     1.00 新規作成
#
# Copyright(C) 2016 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# 引数 : 日時(TT(地球時))
#          書式:YYYYMMDD or YYYYMMDDHHMMSS
#          無指定なら現在(システム日時)を地球時とみなす。
#---------------------------------------------------------------------------------
#++
require 'date'

class NutationModel
  DAT_LS = "nut_ls.txt"
  DAT_PL = "nut_pl.txt"
  PI     = 3.141592653589793238462643     # PI
  PI2    = 6.283185307179586476925287     # 2 * PI
  AS2R   = 4.848136811095359935899141e-6  # Arcseconds to radians
  TURNAS = 1296000.0                      # Arcseconds in a full circle
  U2R    = AS2R / 1e7                     # Units of 0.1 microarcsecond to radians
  R2D    = 57.29577951308232087679815     # Radians to degrees
  D2S    = 3600.0                         # Degrees to seconds

  def initialize
    get_now
    get_arg
    @dat_ls = get_data(DAT_LS)
    @dat_pl = get_data(DAT_PL)
  end

  def calc
    begin
      jd  = calc_jd(@tt)
      t   = calc_t(jd)
      dpsi_ls, deps_ls = calc_lunisolar(t)
      dpsi_pl, deps_pl = calc_planetary(t)
      dpsi, deps = dpsi_ls + dpsi_pl, deps_ls + deps_pl
      dpsi_d, deps_d = dpsi * R2D, deps * R2D
      dpsi_s, deps_s = dpsi_d * D2S, deps_d * D2S
      # Results
      puts @tt.strftime("  [%Y-%m-%d %H:%M:%S TT]")
      puts "  DeltaPsi = #{dpsi} rad"
      puts "           = #{dpsi_d} °"
      puts "           = #{dpsi_s} ″"
      puts "  DeltaEps = #{deps} rad"
      puts "           = #{deps_d} °"
      puts "           = #{deps_s} ″"
    rescue => e
      msg = "[#{e.class}] #{e.message}\n"
      msg << e.backtrace.map { |tr| "\t#{tr}"}.join("\n")
      $stderr.puts msg
      exit 1
    end
  end

  private

  # Getting a current time
  def get_now
    t = Time.now
    @tt = Time.new(t.year, t.month, t.day, t.hour, t.min, t.sec)
  rescue => e
    raise
  end

  # コマンドライン引数の取得
  def get_arg
    return unless arg = ARGV.shift
    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
    exit 0 unless Date.valid_date?(year, month, day)
    exit 0 if hour > 23 || min > 59 || sec > 59
    @tt = Time.new(year, month, day, hour, min, sec)
  rescue => e
    raise
  end

  # ユリウス日の計算
  def calc_jd(tt)
    year, month, day = tt.year, tt.month, tt.day
    hour, min, sec   = tt.hour, tt.min, tt.sec

    begin
      # 1月,2月は前年の13月,14月とする
      if month < 3
        year  -= 1
        month += 12
      end
      # 日付(整数)部分計算
      jd  = (365.25 * year).floor       \
          + (year / 400.0).floor        \
          - (year / 100.0).floor        \
          + (30.59 * (month - 2)).floor \
          + day                         \
          + 1721088.5
      # 時間(小数)部分計算
      t  = (sec / 3600.0 + min / 60.0 + hour) / 24.0
      return jd + t
    rescue => e
      raise
    end
  end

  # ユリウス世紀数の計算
  def calc_t(jd)
    return (jd - 2451545) / 36525.0
  rescue => e
    raise
  end

  # テキストファイルからデータ取得
  # * luni-solar の最初の5列、planetary の最初の14列は整数に、
  #   残りの列は浮動小数点*10000にする
  def get_data(file)
    c = file == DAT_LS ? 5 : 14
    return File.open(file, "r").read.split("\n")[1..-1].map do |l|
      l = l.sub(/^\s+/, "").split(/\s+/)
      l[0, c].map(&:to_i) + l[c..-1].map { |a| a.sub(/\./, "").to_i }
    end
  rescue => e
    raise
  end

  # 日月章動(luni-solar nutation)の計算
  def calc_lunisolar(t)
    dp, de = 0.0, 0.0

    begin
      l  = calc_l_iers2003(t)
      lp = calc_lp_mhb2000(t)
      f  = calc_f_iers2003(t)
      d  = calc_d_mhb2000(t)
      om = calc_om_iers2003(t)
      @dat_ls.reverse.each do |x|
        arg = (x[0] * l + x[1] * lp + x[2] * f +
               x[3] * d + x[4] * om) % PI2
        sarg, carg = Math.sin(arg), Math.cos(arg)
        dp += (x[5] + x[6] * t) * sarg + x[ 7] * carg
        de += (x[8] + x[9] * t) * carg + x[10] * sarg
      end
      return [dp * U2R, de * U2R]
    rescue => e
      raise
    end
  end

  # 惑星章動(planetary nutation)
  def calc_planetary(t)
    dp, de = 0.0, 0.0

    begin
      l   = calc_l_mhb2000(t)
      f   = calc_f_mhb2000(t)
      d   = calc_d_mhb2000_2(t)
      om  = calc_om_mhb2000(t)
      pa  = calc_pa_iers2003(t)
      lme = calc_lme_iers2003(t)
      lve = calc_lve_iers2003(t)
      lea = calc_lea_iers2003(t)
      lma = calc_lma_iers2003(t)
      lju = calc_lju_iers2003(t)
      lsa = calc_lsa_iers2003(t)
      lur = calc_lur_iers2003(t)
      lne = calc_lne_mhb2000(t)
      @dat_pl.reverse.each do |x|
        arg = (x[ 0] * l   + x[ 2] * f   + x[ 3] * d   + x[ 4] * om  +
               x[ 5] * lme + x[ 6] * lve + x[ 7] * lea + x[ 8] * lma +
               x[ 9] * lju + x[10] * lsa + x[11] * lur + x[12] * lne +
               x[13] * pa) % PI2
        sarg, carg = Math.sin(arg), Math.cos(arg)
        dp += x[14] * sarg + x[15] * carg
        de += x[16] * sarg + x[17] * carg
      end
      return [dp * U2R, de * U2R]
    rescue => e
      raise
    end
  end

  # Mean anomaly of the Moon (IERS 2003)
  def calc_l_iers2003(t)
    return ((    485868.249036  + \
            (1717915923.2178    + \
            (        31.8792    + \
            (         0.051635  + \
            (        -0.00024470) \
            * t) * t) * t) * t) % TURNAS) * AS2R
  rescue => e
    raise
  end

  # Mean anomaly of the Sun (MHB2000)
  def calc_lp_mhb2000(t)
    return ((  1287104.79305   + \
            (129596581.0481    + \
            (       -0.5532    + \
            (        0.000136  + \
            (       -0.00001149) \
            * t) * t) * t) * t) % TURNAS) * AS2R
  rescue => e
    raise
  end

  # Mean longitude of the Moon minus that of the ascending node (IERS 2003)
  def calc_f_iers2003(t)
    return ((     335779.526232 + \
            (1739527262.8478    + \
            (       -12.7512    + \
            (        -0.001037  + \
            (         0.00000417) \
            * t) * t) * t) * t) % TURNAS) * AS2R
  rescue => e
    raise
  end

  # Mean elongation of the Moon from the Sun (MHB2000)
  def calc_d_mhb2000(t)
    return ((   1072260.70369   + \
            (1602961601.2090    + \
            (        -6.3706    + \
            (         0.006593  + \
            (        -0.00003169) \
            * t) * t) * t) * t) % TURNAS) * AS2R
  rescue => e
    raise
  end

  # Mean longitude of the ascending node of the Moon (IERS 2003)
  def calc_om_iers2003(t)
   return ((    450160.398036  + \
           (  -6962890.5431    + \
           (         7.4722    + \
           (         0.007702  + \
           (        -0.00005939) \
           * t) * t) * t) * t) % TURNAS) * AS2R
  rescue => e
    raise
  end

  # Mean anomaly of the Moon (MHB2000)
  def calc_l_mhb2000(t)
    return (2.35555598 + 8328.6914269554 * t) % PI2
  rescue => e
    raise
  end

  # Mean longitude of the Moon minus that of the ascending node (MHB2000)
  def calc_f_mhb2000(t)
    return (1.627905234 + 8433.466158131 * t) % PI2
  rescue => e
    raise
  end

  # Mean elongation of the Moon from the Sun (MHB2000)
  def calc_d_mhb2000_2(t)
    return (5.198466741 + 7771.3771468121 * t) % PI2
  rescue => e
    raise
  end

  # Mean longitude of the ascending node of the Moon (MHB2000)
  def calc_om_mhb2000(t)
    return (2.18243920 - 33.757045 * t) % PI2
  rescue => e
    raise
  end

  # General accumulated precession in longitude (IERS 2003)
  def calc_pa_iers2003(t)
    return (0.024381750 + 0.00000538691 * t) * t
  rescue => e
    raise
  end

  # Mercury longitudes (IERS 2003)
  def calc_lme_iers2003(t)
    return (4.402608842 + 2608.7903141574 * t) % PI2
  rescue => e
    raise
  end

  # Venus longitudes (IERS 2003)
  def calc_lve_iers2003(t)
    return (3.176146697 + 1021.3285546211 * t) % PI2
  rescue => e
    raise
  end

  # Earth longitudes (IERS 2003)
  def calc_lea_iers2003(t)
    return (1.753470314 + 628.3075849991 * t) % PI2
  rescue => e
    raise
  end

  # Mars longitudes (IERS 2003)
  def calc_lma_iers2003(t)
    return (6.203480913 + 334.0612426700 * t) % PI2
  rescue => e
    raise
  end

  # Jupiter longitudes (IERS 2003)
  def calc_lju_iers2003(t)
    return (0.599546497 + 52.9690962641 * t) % PI2
  rescue => e
    raise
  end

  # Saturn longitudes (IERS 2003)
  def calc_lsa_iers2003(t)
    return (0.874016757 + 21.3299104960 * t) % PI2
  rescue => e
    raise
  end

  # Uranus longitudes (IERS 2003)
  def calc_lur_iers2003(t)
    return (5.481293872 + 7.4781598567 * t) % PI2
  rescue => e
    raise
  end

  # Neptune longitude (MHB2000)
  def calc_lne_mhb2000(t)
    return (5.321159000 + 3.8127774000 * t) % PI2
  rescue => e
    raise
  end
end

nm = NutationModel.new
nm.calc

参考までに、スクリプト中に出現する用語の意味は以下のとおり。

1
2
3
4
5
6
7
8
9
10
11
* Mean anomaly of the Moon                                    (平均近点角:月)
* Mean longitude of the Moon minus that of the ascending node (平均黄経:月(昇交点を引いたもの))
* Mean longitude of the Moon's ascending node                 (平均黄経:月の昇交点)
* Mean longitude of Mercury                                   (平均黄経:水星)
* Mean longitude of Venus                                     (平均黄経:金星)
* Mean longitude of Earth                                     (平均黄経:地球)
* Mean longitude of Mars                                      (平均黄経:火星)
* Mean longitude of Jupiter                                   (平均黄経:木星)
* Mean longitude of Saturn                                    (平均黄経:土星)
* Mean longitude of Uranus                                    (平均黄経:天王星)
* General accumulated precession in longitude                 (一般歳差:J2000.0 からの蓄積)

3. Ruby スクリプトの実行

YYYYMMDD or YYYYMMDDHHMMSS で地球時(TT) を指定して実行する。(引数なしでシステム時刻を地球時とみなす)

1
2
3
4
5
6
7
8
./nutation_model.rb 20160528000000
  [2016-05-28 00:00:00 TT]
  DeltaPsi = -2.1633261603237626e-05 rad
           = -0.0012394945869679329 °
           = -4.462180513084558 ″
  DeltaEps = -4.591610134570194e-05 rad
           = -0.002630798818803681 °
           = -9.47087574769325 ″

4. 計算結果の検証

国立天文台・暦象年表のツールで計算した値と比較してみたが、一致することが確認できた。

5. 参考サイト

(USNO(米・海軍天文台)のサイトは、環境(プロキシ設定)によっては閲覧不可)


今回の章動の計算も、天体位置やこよみを正確に計算する際に必要になってきます。

当方、今後挑戦してみようと考えております。(可能ならば)

以上。

Comments