Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Ruby script to calculate apparent position of Sun and Moon.
#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= 太陽・月の視位置計算
#
# * JPLEPH(JPL の DE430 バイナリデータ)読み込み、視位置を計算する
# (自作 RubyGems ライブラリ mk_apos を使用)
#
# date name version
# 2016.09.14 mk-mode.com 1.00 新規作成
# 2018.06.08 mk-mode.com 1.01 視半径/(地平)視差の取得を追加
#
# Copyright(C) 2016 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# 引数 : 日時(JST)
# 書式:YYYYMMDD or YYYYMMDDHHMMSS
# 無指定なら現在(システム日時)と判断。
#---------------------------------------------------------------------------------
#++
require 'date'
require 'mk_apos'
class ApparentSunMoon
USAGE = "[USAGE] ./apparent_sun_moon.rb [YYYYMMDD|YYYYMMDDHHMMSS|YYYYMMDDHHMMSSU...]"
FILE_BIN = "/home/masaru/src/calendar/ephemeris_jpl/JPLEPH" # JPL バイナリデータ
AU = 149597870.7 # 1天文単位 (km)
#=========================================================================
# 初期処理
#
# @param: <none>
# @return: <none>
#=========================================================================
def initialize
@jst = get_arg
@utc = jst2utc(@jst)
end
#=========================================================================
# 主処理
#
# @param: <none>
# @return: <none>
#=========================================================================
def exec
# 計算
@apos = MkApos.new(FILE_BIN, @utc.strftime("%Y%m%d%H%M%S%6N"))
sun = @apos.sun
moon = @apos.moon
# 結果出力
puts "* 計算対象時刻 (JST) = #{@jst.strftime("%Y-%m-%d %H:%M:%S.%3N")}"
puts " (UTC) = #{@apos.utc.strftime("%Y-%m-%d %H:%M:%S.%3N")}"
puts " (TDB) = #{@apos.tdb.strftime("%Y-%m-%d %H:%M:%S.%3N")}"
puts " (JD) = #{@apos.jd_tdb}"
puts "* 視位置:太陽"
printf(" = [赤経: %14.10f rad, 赤緯: %14.10f rad]\n", sun[0][0], sun[0][1])
printf(" = [赤経: %14.10f deg, 赤緯: %14.10f deg]\n", sun[0][0] * 180 / Math::PI, sun[0][1] * 180 / Math::PI)
printf(" = [黄経: %14.10f rad, 黄緯: %14.10f rad]\n", sun[1][0], sun[1][1])
printf(" = [黄経: %14.10f deg, 黄緯: %14.10f deg]\n", sun[1][0] * 180 / Math::PI, sun[1][1] * 180 / Math::PI)
puts "* 視位置:月"
printf(" = [赤経: %14.10f rad, 赤緯: %14.10f rad]\n", moon[0][0], moon[0][1])
printf(" = [赤経: %14.10f deg, 赤緯: %14.10f deg]\n", moon[0][0] * 180 / Math::PI, moon[0][1] * 180 / Math::PI)
printf(" = [黄経: %14.10f rad, 黄緯: %14.10f rad]\n", moon[1][0], moon[1][1])
printf(" = [黄経: %14.10f deg, 黄緯: %14.10f deg]\n", moon[1][0] * 180 / Math::PI, moon[1][1] * 180 / Math::PI)
puts "* 視黄経差:太陽 - 月"
printf(" = %.10f rad = %.10f deg\n", sun[1][0] - moon[1][0], (sun[1][0] - moon[1][0]) * 180 / Math::PI)
puts "* 距離:太陽"
printf(" = %.10f AU = %.2f km\n", sun[0][2], sun[0][2] * AU)
puts "* 距離:月"
printf(" = %.10f AU = %.2f km\n", moon[0][2], moon[0][2] * AU)
puts "* 視半径:太陽"
printf(" = %.2f ″\n", sun[2][0])
puts "* 視半径:月"
printf(" = %.2f ″\n", moon[2][0])
puts "* (地平)視差:太陽"
printf(" = %.2f ″\n", sun[2][1])
puts "* (地平)視差:月"
printf(" = %.2f ″\n", moon[2][1])
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: UTC (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
usec = arg.split(//).size > 14 ? arg[14..-1] : "0"
(puts USAGE; exit 0) unless Date.valid_date?(year, month, day)
(puts USAGE; exit 0) if hour > 23 || min > 59 || sec > 59
d = usec.split(//).size
return Time.new(year, month, day, hour, min, sec + Rational(usec.to_i, 10 ** d), "+09:00")
rescue => e
raise
end
#=========================================================================
# JST -> UTC
#
# @param: JST (Time Object)
# @return: UTC (Time Object)
#=========================================================================
def jst2utc(jst)
return jst - 9 * 60 * 60
rescue => e
raise
end
end
ApparentSunMoon.new.exec if __FILE__ == $0
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.