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: 雑記 | No Comments »
ハクビシン(?) 目撃
By kmgs | 7 月 10, 2009
7/10 23:50ごろ、マンションの庭から"ケケケケケ、ケケケケケ"と奇妙な声がしました。最初は鳥かと思ったのですが、カミさんと覗いてみると、体長40cmほど、尻尾が20cm程度のタヌキのようなズングリした動物がゆっくり歩いていました。
動物の声は「ハクビシンの不安なときの声」(泣き声大全集)と似ていました。文京区に引っ越してきて3年半になりますが、こんな動物に遭遇するとは予想していませんでした。写真かビデオにでも残しておけば良かった...。
調べてみたら、23区内でも目撃例(区役所への聞き取り調査、東京タヌキ探検隊など)は少なくないようですね。
そういえば、今年の春先ごろから雨の日になるとマンションの駐車場の植え込みに体長10cmほどのヒキガエルが出現したり、先々週あたりは神田明神近辺でシマヘビに遭遇したりと、意外な場所で意外な動物に遭遇することが多いです。
Topics: 雑記 | No Comments »
PICで正弦波 (6)
By kmgs | 6 月 19, 2009
このへん(その1、その2、その3、その4、その5)の続き。
デルタシグマ変調を2次にしてみたところ、THD 0.005%、THD+N 0.23%となりました。もう少しTHD+Nを追い込めないものか、検討中です。ちなみにSPI出力は、現状で1.667Mbps (インストラクションサイクル40MHz の1/24)です。
Topics: PIC24HJ12GP202 | No Comments »
PICで正弦波 (5)
By kmgs | 6 月 17, 2009
このへん(その1、その2、その3、その4)の続き。
SPI出力を直接LPFで取り出した場合、2次の高調波が出現するため、原因が分からずにしばらく悩んでいました。PICのソースをいじってDDS化を試みたりしていたのですが、2次の高調波はデルタシグマ変調に関わる問題ではなかったようです。
オシロでSPI出力を追っていたら、High側の出力電圧が荒れていることに気づき、ためしにクランプダイオードで適当に落としてみたところ ダイオードでリミッタを掛けたところ、それなりの改善が見られました。
↑ リミット前: 2kHzが-64dBくらい
↑ リミット後: 2kHzが-84dBくらい
リミッタの条件は検討中ですが、現時点ではSPI出力に100Ωを噛ませて、小信号用シリコンダイオード->Vcc/2 (1kΩ x 2で作った中点)でクランプした出力を12kΩ+0.01uFのLPFに通し、ボルテージフォロアに入力しています。小信号用シリコンダイオードは1S1588と1N4148を試してみましたが、大差ないようです(手元の整流用ショットキーダイオードでは、どうやら波形がなまって逆効果)。
2009/06/19訂正: クランプ回路も試したのですが、最終的にはリミッタ回路を採用しました。図はリミッタ回路のものです。
[ その6 ]
Topics: PIC24HJ12GP202 | No Comments »
PICで正弦波 (4)
By kmgs | 6 月 5, 2009
このへん(その1、その2、その3)の続き。
だんだんマトモになってきました。現状でTHD 0.04%、THD+N 0.26%くらいです。
[ その5へ ]
Topics: PIC24HJ12GP202 | No Comments »
PICで三角波と鋸波
By kmgs | 6 月 2, 2009
このへん(その1、その2、その3)の続き。
2値量子化2次デルタシグマ変調(1.024MHz)をPIC上でリアルタイムで行い、SPIを使って1kHzの三角波と鋸波を出力してみました。三角波はそれなりの形をしていますが、少し鈍り気味です。鋸波はさらに鈍っており、ゼロ点付近で崩れています。コードがバグっているのか、LPFの調整が適当すぎるのか。
もう少し追い込みが必要そうです。
[ 続き ]
Topics: PIC24HJ12GP202 | No Comments »