Ruby - 地球自転速度補正値 ΔT の計算！

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

地球自転速度補正値 delta T (ΔT) の計算式が NASA - Polynomial Expressions for Delta T に掲載されていますので、 Ruby で実装してみました。(-1900年〜3000年対応）

さらに、1972年から2018年までは、うるう秒実施日一覧: NICT - 日本標準時プロジェクトで公開されている「うるう秒」に 32.184(= TT（地球時） - TAI（国際原子時）) を加算した値とも比較できるようにしてみました。

各種時刻系については、当ブログ過去記事等をご参照ください。

0. 前提条件

  • Ruby 2.3.1-p112 での作業を想定。
  • NASA の計算式で計算できる西暦は -1900年から3000年。
    （AD元年の前年はBC元年なので、今回の計算上の 0年 は BC1年、 -1900年 は BC1901年）
  • うるう秒が実際されたのは1972年から2017年（当記事執筆時点）。（うるう秒実施日一覧: NICT - 日本標準時プロジェクト）
    次回、いつ「うるう秒」が挿入（or 削除）されるかは不明なので、以下のスクリプト内の NICT 版メソッドでの対応は2019年直前までとする。
  • 「ΔT = TT - UT1 = TAI - UTC - DUT1 + 32.184 = 32.184 - (UTC - TAI) - DUT1 = 32.184 - うるう秒総和 - DUT1」であることから、本来は DUT1(= UT1 - UTC) も考慮すべきであるが、今回は考慮しない。（誤差は大きくても 1 秒未満。参照「日本標準時プロジェクト　Announcement of DUT1」）

1. Ruby スクリプトの作成

calc_delta_t.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

#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= 地球自転速度の補正値 delta T(ΔT)の計算
#  : [NASA - Polynomial Expressions for Delta T](http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html)
#    の計算式を使用する。
#   1972年 - 2018年は、比較対象として「うるう年総和 + 32.184(TT - TAI)」の値も計算する。
#
# date          name            version
# 2016.06.15    mk-mode.com     1.00 新規作成
# 2016.07.20    mk-mode.com     1.01 第27回うるう秒調整に関する判定を追加
#
# Copyright(C) 2016 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# 引数 : YYYYMM
#        * YYYYMM は UTC
#        * 無指定なら現在年月を UTC とみなす。
#---------------------------------------------------------------------------------
#++
require 'date'

class CalcDeltaT
  USAGE     = "[USAGE] ./calc_delta_t.rb [[+-]YYYYMM]"
  MSG_ERR_1 = "[ERROR] Year must be between -1999 and 3000!"
  MSG_ERR_2 = "[ERROR] Month must be between 1 and 12!"
  TT_TAI    = 32.184

  def initialize
    ym = ARGV.shift
    ym ||= Time.now.strftime("%Y%m")
    unless ym =~ /^[+-]?\d{6}$/
      puts USAGE
      exit 0
    end
    @year, @month = ym.scan(/([+-]?\d{4})(\d{2})/)[0].map(&:to_i)
    if @year < -1999 || @year > 3000
      puts MSG_ERR_1
      exit 0
    end
    if @month < 1 || @month > 12
      puts MSG_ERR_2
      exit 0
    end
    @y = @year + (@month - 0.5) / 12
  end

  def calc
    # NASA Ver. (-1999 - 3000)
    case
    when                  @year <  -500; dt = calc_before_m500
    when -500 <= @year && @year <   500; dt = calc_before_500
    when  500 <= @year && @year <  1600; dt = calc_before_1600
    when 1600 <= @year && @year <  1700; dt = calc_before_1700
    when 1700 <= @year && @year <  1800; dt = calc_before_1800
    when 1800 <= @year && @year <  1860; dt = calc_before_1860
    when 1860 <= @year && @year <  1900; dt = calc_before_1900
    when 1900 <= @year && @year <  1920; dt = calc_before_1920
    when 1920 <= @year && @year <  1941; dt = calc_before_1941
    when 1941 <= @year && @year <  1961; dt = calc_before_1961
    when 1961 <= @year && @year <  1986; dt = calc_before_1986
    when 1986 <= @year && @year <  2005; dt = calc_before_2005
    when 2005 <= @year && @year <  2050; dt = calc_before_2050
    when 2050 <= @year && @year <= 2150; dt = calc_until_2150
    when 2150 <  @year                 ; dt = calc_after_2150
    end
    str = sprintf("[%04d-%02d]", @year, @month)
    str << " delta T = #{dt}"

    # NICT Ver.  (1972-01-01 - 2018-12-31)
    if 1972 <= @year && @year <= 2018
      str << " (NICT: #{calc_nict})"
    end
    puts str
  rescue => e
    $stderr.puts "[#{e.class}] #{e.message}"
    e.backtrace.each { |tr| $stderr.puts "\t#{tr}"}
    exit 1
  end

  private

  # year < -500
  def calc_before_m500
    t = (@y - 1820) / 100.0
    dt = -20 + 32 * t ** 2
    return dt
  rescue => e
    raise
  end

  # -500 <= year && year < 500
  def calc_before_500
    t = @y / 100.0
    dt  = 10583.6         + \
         (-1014.41        + \
         (   33.78311     + \
         (   -5.952053    + \
         (   -0.1798452   + \
         (    0.022174192 + \
         (    0.0090316521) \
         * t) * t) * t) * t) * t) * t
    return dt
  rescue => e
    raise
  end

  # 500 <= year && year < 1600
  def calc_before_1600
    t = (@y - 1000) / 100.0
    dt  = 1574.2         + \
         (-556.01        + \
         (  71.23472     + \
         (   0.319781    + \
         (  -0.8503463   + \
         (  -0.005050998 + \
         (   0.0083572073) \
         * t) * t) * t) * t) * t) * t
    return dt
  rescue => e
    raise
  end

  # 1600 <= year && year < 1700
  def calc_before_1700
    t = @y - 1600
    dt  = 120           + \
         ( -0.9808      + \
         ( -0.01532     + \
         (  1.0 / 7129.0) \
         * t) * t) * t
    return dt
  rescue => e
    raise
  end

  # 1700 <= year && year < 1800
  def calc_before_1800
    t = @y - 1700
    dt  =  8.83           + \
         ( 0.1603         + \
         (-0.0059285      + \
         ( 0.00013336     + \
         (-1.0 / 1174000.0) \
         * t) * t) * t) * t
    return dt
  rescue => e
    raise
  end

  # 1800 <= year && year < 1860
  def calc_before_1860
    t = @y - 1800
    dt  = 13.72          + \
         (-0.332447      + \
         ( 0.0068612     + \
         ( 0.0041116     + \
         (-0.00037436    + \
         ( 0.0000121272  + \
         (-0.0000001699  + \
         ( 0.000000000875) \
         * t) * t) * t) * t) * t) * t) * t
    return dt
  rescue => e
    raise
  end

  # 1860 <= year && year < 1900
  def calc_before_1900
    t = @y - 1860
    dt  =  7.62          + \
         ( 0.5737        + \
         (-0.251754      + \
         ( 0.01680668    + \
         (-0.0004473624  + \
         ( 1.0 / 233174.0) \
         * t) * t) * t) * t) * t
    return dt
  rescue => e
    raise
  end

  # 1900 <= year && year < 1920
  def calc_before_1920
    t = @y - 1900
    dt  = -2.79      + \
         ( 1.494119  + \
         (-0.0598939 + \
         ( 0.0061966 + \
         (-0.000197  ) \
         * t) * t) * t) * t
    return dt
  rescue => e
    raise
  end

  # 1920 <= year && year < 1941
  def calc_before_1941
    t = @y - 1920
    dt  = 21.20     + \
         ( 0.84493  + \
         (-0.076100 + \
         ( 0.0020936) \
         * t) * t) * t
    return dt
  rescue => e
    raise
  end

  # 1941 <= year && year < 1961
  def calc_before_1961
    t = @y - 1950
    dt  = 29.07      + \
         ( 0.407     + \
         (-1 / 233.0 + \
         ( 1 / 2547.0) \
         * t) * t) * t
    return dt
  rescue => e
    raise
  end

  # 1961 <= year && year < 1986
  def calc_before_1986
    t = @y - 1975
    dt = 45.45      + \
        ( 1.067     + \
        (-1 / 260.0 + \
        (-1 / 718.0)  \
        * t) * t) * t
    return dt
  rescue => e
    raise
  end

  # 1986 <= year && year < 2005
  def calc_before_2005
    t = @y - 2000
    dt  = 63.86         + \
         ( 0.3345       + \
         (-0.060374     + \
         ( 0.0017275    + \
         ( 0.000651814  + \
         ( 0.00002373599) \
         * t) * t) * t) * t) * t
    return dt
  rescue => e
    raise
  end

  # 2005 <= year && year < 2050
  def calc_before_2050
    t = @y - 2000
    dt  = 62.92    + \
         ( 0.32217 + \
         ( 0.005589) \
         * t) * t
    return dt
  rescue => e
    raise
  end

  # 2050 <= year && year <= 2150
  def calc_until_2150
    dt  = -20 \
        + 32 * ((@y - 1820) / 100.0) ** 2
        - 0.5628 * (2150 - @y)
    return dt
  rescue => e
    raise
  end

  # 2150 < year
  def calc_after_2150
    t = (@y - 1820) / 100.0
    dt  = -20 + 32 * t ** 2
    return dt
  rescue => e
    raise
  end

  # NICT Ver.  (1972-01-01 - 2017-12-31)
  def calc_nict
    ym = sprintf("%04d%02d", @year, @month)
    case
    when ym < "197207"; dt = TT_TAI + 10
    when ym < "197301"; dt = TT_TAI + 11
    when ym < "197401"; dt = TT_TAI + 12
    when ym < "197501"; dt = TT_TAI + 13
    when ym < "197601"; dt = TT_TAI + 14
    when ym < "197701"; dt = TT_TAI + 15
    when ym < "197801"; dt = TT_TAI + 16
    when ym < "197901"; dt = TT_TAI + 17
    when ym < "198001"; dt = TT_TAI + 18
    when ym < "198107"; dt = TT_TAI + 19
    when ym < "198207"; dt = TT_TAI + 20
    when ym < "198307"; dt = TT_TAI + 21
    when ym < "198507"; dt = TT_TAI + 22
    when ym < "198801"; dt = TT_TAI + 23
    when ym < "199001"; dt = TT_TAI + 24
    when ym < "199101"; dt = TT_TAI + 25
    when ym < "199207"; dt = TT_TAI + 26
    when ym < "199307"; dt = TT_TAI + 27
    when ym < "199407"; dt = TT_TAI + 28
    when ym < "199601"; dt = TT_TAI + 29
    when ym < "199707"; dt = TT_TAI + 30
    when ym < "199901"; dt = TT_TAI + 31
    when ym < "200601"; dt = TT_TAI + 32
    when ym < "200901"; dt = TT_TAI + 33
    when ym < "201207"; dt = TT_TAI + 34
    when ym < "201507"; dt = TT_TAI + 35
    when ym < "201701"; dt = TT_TAI + 36
    when ym < "201901"; dt = TT_TAI + 37  # <= 第28回うるう秒実施までの暫定措置
    end
    return dt
  rescue => e
    raise
  end
end

CalcDeltaT.new.calc if __FILE__ == $0

t のべき乗計算部分を特殊な記述に変更しているのは、乗算回数を減らして高速化するため。
（例えば、1 + 2 * t + 3 * t**2 + 4 * t**3 だと乗算回数が「6回」だが、1 + (2 + (3 + 4*t)*t)*t と書き換えることで乗算回数が「3回」になる）

2. Ruby スクリプトの実行

1
2
3
4
5
6
7
8
9
10
11
12
13
14

$ ./calc_delta_t.rb 000112
[0001-12] delta T = 10563.74738216227

$ ./calc_delta_t.rb 190007
[1900-07] delta T = -1.997290733102346

$ ./calc_delta_t.rb 200007
[2000-07] delta T = 64.02380536456008 (NICT: 64.184)

$ ./calc_delta_t.rb 201607
[2016-07] delta T = 69.77852857812505 (NICT: 68.184)

$ ./calc_delta_t.rb 300007
[3000-07] delta T = 4439.771605555555

NASA 版の計算式では細かい補正が行われてはいますが、当方は、うるう秒挿入（or 削除）が確実な場合は「うるう秒挿入（or 削除）総和 + 32.184(TT - TAI)」を ΔT とするようにしています。

以上。