mk-mode BLOG

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

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

Ruby - 平均黄道傾斜角の計算!

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

こんばんは。

当ブログの以前の記事「黄道傾斜角について!」を元に、平均黄道傾斜角の計算を Ruby で実装してみました。(ただそれだけ)

0. 前提条件

1. Ruby スクリプトの作成

mean_obliquity_ecliptic.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
#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= 平均黄道傾斜角の計算
#
# date          name            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 MeanObliquityEcliptic
  def initialize
    get_now
    get_arg
  end

  def calc
    jd  = calc_jd(@tt)
    t   = calc_t(jd)
    moe = calc_eps_a(t)
    puts "           地球時(TT): #{@tt.strftime("%Y-%m-%d %H:%M:%S")}"
    puts "       ユリウス日(JD): #{jd}"
    puts "   ユリウス世紀数(JD): #{t}"
    puts "平均黄道傾斜角(eps_a): #{moe} °"
  rescue => e
    $stderr.puts "[#{e.class}] #{e.message}"
    e.backtrace.each { |tr| $stderr.puts "\t#{tr}"}
    exit 1
  end

  private

  # 現在日時の取得
  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

  # 黄道傾斜角(εA)の計算
  def calc_eps_a(t)
    return (84381.406           + \
           (  -46.836769        + \
           (   -0.0001831       + \
           (    0.00200340      + \
           (   -5.76 * 10 ** -7 + \
           (   -4.34 * 10 ** -8)  \
           * t) * t) * t) * t) * t) / 3600.0
  rescue => e
    raise
  end
end

exit unless __FILE__ == $0
MeanObliquityEcliptic.new.calc

calc_eps_a メソッド内の t のべき乗計算部分を特殊な記述にしているのは、演算コストのかかる乗算回数を減らして高速化を図るため。(例えば、1 + 2 * t + 3 * t**2 + 4 * t**3 だと乗算回数が「6回」だが、1 + (2 + (3 + 4 * t) * t) * t と書き換えることで乗算回数が「3回」になる)

2. Ruby スクリプトの実行

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

1
2
3
4
5
$ ./mean_obliquity_ecliptic.rb 20160526
           地球時(TT): 2016-05-26 00:00:00
       ユリウス日(JD): 2457534.5
   ユリウス世紀数(JD): 0.1639835728952772
平均黄道傾斜角(eps_a): 23.437145984218514 °

3. 計算結果の検証

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

4. 参考サイト


高精度な結果が得られるので、天体位置や暦の正確な計算に利用できそうです。

以上。

Comments