« ハクビシン(?) 目撃 | Home
Rで部分日食の食の入り/食の終わりの時間を計算してみる
By kmgs | 7 月 21, 2009
ふと、部分日食の時間を計算で予測できないものかと考え、調べてみたところ山口大学 藤沢研究室に丁寧な解説がありました(PDF)。このドキュメントでは、日食の計算に必要なベッセル要素を数表から出しているのですが、NASAにベッセル要素の近似式パラメータがあったので、試しにそれを使ってみることにしました。作製した関数は以下の通り。
test <- function(phi0=35.658, lambda0=139.741, h=0) {
# NASAから取ってきたパラメータ (Besselian elements)
# http://eclipse.gsfc.nasa.gov/SEbeselm/SEbeselm2001/SE2009Jul22Tbeselm.html
tanf1 <- 0.0046014
tanf2 <- 0.0045784
deltaT <- 66.5
x <- function(t) {
0.239961 + 0.5563955 * t - 0.0000576 * t^2 - 0.0000094 * t^3
}
y <- function(t) {
-0.003276 - 0.1774579 * t - 0.0001344 * t^2 + 0.0000032 * t^3
}
d <- function(t) {
(20.26424 - 0.007874 * t - 0.000005 * t^2) / 180 * pi
}
l1 <- function(t) {
0.530427 + 0.0000063 * t - 0.0000128 * t^2
}
l2 <- function(t) {
-0.015633 + 0.0000063 * t - 0.0000127 * t^2
}
mu <- function(t) {
223.38821 + 15.001004 * t
}
# NASAパラメータ 終わり
# 日食の計算、全面的に以下URLを参考 (山口大学 藤沢研究室)
# http://www.astro.sci.yamaguchi-u.ac.jp/~kenta/eclipse/
Q <- function(t, lambda, s2, c2) {
mu_lambda <- (mu(t) + lambda) / 180 * pi
dt <- d(t)
cosd <- cos(dt)
sind <- sin(dt)
xi <- c2 * sin(mu_lambda)
eta <- s2 * cosd - c2 * sind * cos(mu_lambda)
zeta <- s2 * sind + c2 * cosd * cos(mu_lambda)
L1 <- l1(t) - zeta * tanf1
L2 <- l2(t) - zeta * tanf2
delta <- sqrt((x(t) - xi)^2 + (y(t) - eta)^2)
Q1 <- L1^2 - delta^2
Q2 <- L2^2 - delta^2
D <- (L1 - delta) / (L1 + L2)
t1 <- (t + 3 + 9) * 3600 - deltaT + strptime("2009-07-22", format="%Y-%m-%d")
op <- options(digits.secs=2)
z <- format(t1, "%H:%M:%OS")
options(op)
data.frame(t, z, xi, eta, zeta, L1, L2, delta, Q1, Q2, D)
}
# http://www.astro.sci.yamaguchi-u.ac.jp/~kenta/eclipse/Eclipse.pdf
# 11ページのS値、C値の計算
phi <- phi0 / 180 * pi
s0 <- 0.99497430 - 0.00167078 * cos(phi * 2) + 0.00000210 * cos(phi * 4)
c0 <- 1.00167993 - 0.00168204 * cos(phi * 2) + 0.00000212 * cos(phi * 4)
s1 <- s0 + h * 0.0000001568
c1 <- c0 + h * 0.0000001568
s2 <- s1 * sin(phi)
c2 <- c1 * cos(phi)
# 11ページの暦表経度の計算
lambda <- lambda0 - (15.0410685 * deltaT) / 3600
# 日食の計算を実行
res <- Q(seq(-3, 3, 1/14400), lambda, s2, c2)
# 食の始め、食の最大、食の終わり
wm <- which.max(res$Q1)
start <- which.min(abs(res$Q1[1:wm]))
end <- which.min(abs(res$Q1[wm:nrow(res)])) + wm - 1
res[c(start, wm, end),]
}こいつをRのコンソールにペーストするなり、ファイルに書き込んでsource()するなりして読み込み、
> test()
t z xi eta zeta L1 L2 delta Q1 Q2 D
13603 -2.0554167 09:55:34.00 -0.3816591 0.2947777 0.8747443 0.5263349 -0.01970453 0.5263334 1.650312e-06 -0.27663854 3.094449e-06
32176 -0.7656250 11:12:57.25 -0.1221117 0.2651806 0.9552503 0.5260192 -0.02001879 0.1472778 2.550054e-01 -0.02128999 7.485002e-01
50749 0.5241667 12:30:20.50 0.1512289 0.2670687 0.9505468 0.5260529 -0.01998517 0.5260528 1.800838e-07 -0.27633211 3.382256e-07のように、デフォルトでは東京(緯度35.658、経度139.741)に対して、食の入り(1行目)、最大食分(2行目)、食の終わり(3行目)を計算します。ここで、zが2009/7/22に対する日本標準時、Dが食分となります。ここでの計算結果は食の入りが09:55:34.00、最大職分0.749 @ 11:12:57.25、食の終わりが12:30:20.50となります(実際には、月面の凹凸や誤差の伝播を考慮すると、有効数字はここまでないと思われます)。
また、他の場所の予測を行うには、引数に緯度/経度/高度を与えます。
# 福岡: 最大食分の時間が国立天文台の予測時間(10:56:03)と4秒くらいずれている。
> test(33.583, 130.400, 0)
t z xi eta zeta L1 L2 delta
9292 -2.3547917 09:37:36.25 -0.55599852 0.3004395 0.7736686 0.5267812 -0.01926042 0.52677753
28104 -1.0484028 10:55:59.25 -0.31532190 0.2484387 0.9147741 0.5261971 -0.01984177 0.07156466
47743 0.3154167 12:17:49.00 -0.02558462 0.2273055 0.9724397 0.5259531 -0.02008449 0.52596000
Q1 Q2 D
9292 3.902627e-06 -0.277123598 7.298681e-06
28104 2.717619e-01 -0.004727805 8.978526e-01
47743 -7.227812e-06 -0.276230536 -1.358280e-05なお、国立天文台の予測と比較すると、食の入り・食の終わりは±1秒、最大食分時間には±5秒くらいのずれがあるようです。どこかでミスっているかもしれませんので、上記プログラムの結果はあくまで目安としてお考えください。
Topics: 雑記 |