Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Ruby script to calculate planetary distances.
#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= JPLEPH(JPL の DE430 バイナリデータ)読込後に座標(位置・速度)を計算し、さらに
# 地球から太陽・惑星・月までの距離を計算する。
# : 自作の RubyGems ライブラリ eph_jpl を使用する
# (JPL DE430 のバイナリデータが準備済みであること)
#
# date name version
# 2016.07.13 mk-mode.com 1.00 新規作成
#
# Copyright(C) 2016 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# 【引数】 ユリウス日
# (省略可。省略時は現在日時のユリウス日を地球時TT(≒太陽系力学時TDB)とみなす)
#---------------------------------------------------------------------------------
#++
require 'eph_jpl'
class JplDistanceEarth
FILE_BIN = "/path/to/JPLEPH"
USAGE = "【使用方法】 ./jpl_distance_earth.rb [ユリウス日]"
CENTER = 3
KM = false # true: km, false: AU
BODIES = {
1 => " 水星 (Mercury)", 2 => " 金星 (Venus) ", 3 => " 地球 (Earth) ",
4 => " 火星 (Mars) ", 5 => " 木星 (Jupiter)", 6 => " 土星 (Saturn) ",
7 => "天王星 (Uranus) ", 8 => "海王星 (Neptune)", 9 => "冥王星 (Pluto) ",
10 => " 月 (Moon) ", 11 => " 太陽 (Sun) "
}
def initialize
@jd = get_args
@unit = KM ? "KM" : "AU"
end
def exec
puts "JD(TT) = #{@jd}\n#{BODIES[3]}"
BODIES.each do |target, target_name|
next if target == 3
eph = EphJpl.new(FILE_BIN, target, CENTER, @jd, KM)
ps = eph.calc[0,3]
d = calc_distance(ps)
d = KM ? sprintf("%19.8f", d) : sprintf("%11.8f", d)
puts "- #{target_name} : #{d} #{@unit}"
end
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: jd
#=========================================================================
def get_args
jd = ARGV.shift
return jd ? jd.to_f : gc2jd({
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
puts USAGE
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` とする。
#
# @param: hash = {
# year(int), month(int), day(int), hour(int), min(int), sec(int)
# }
# @return: jd (ユリウス日)
#=========================================================================
def gc2jd(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 \
+ (year / 400.0).floor \
- (year / 100.0).floor \
+ (30.59 * (month - 2)).floor \
+ day \
+ 1721088
# 時間(小数)部分計算
t = (sec / 3600.0 + min / 60.0 + hour) / 24.0
return jd + t
rescue => e
raise
end
end
#=========================================================================
# 距離計算
#
# @param: ps (x-座標, y-座標, z-座標)
# @return: distance (距離)
#=========================================================================
def calc_distance(ps)
return Math.sqrt(ps.inject(0) { |s, a| s + a * a })
rescue => e
raise
end
end
JplDistanceEarth.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.