EASolverで遊ぶ
EASolverの試しがてら、色々の星の極小データだけから周期を出して見ています。大部分の星の周期はもうほぼ確定しているものばかりでしたが、それがちゃんと出てくるかどうか試したわけです。
今回EASolverを作ろうと思い立った原因はOQ Gemだったのですが、私がたまたま観測した極小（1月20日）がなかったらどうだっただろうかと試してみました。
すでにこの4つの極小時刻が観測されていました。
2438406.6500 15
2439063.5300 15
2453056.4122 15
2453762.4595 15
これをEASolverに入れて周期を求めさせたら、
> P # the best 3 periods
[1] 16.42332566845 9.80575319285 2.13958610840
と出てきました。最初の16.4233256..は私の極小時刻を入れて求められた2.7367288のきっちり６倍となりました。
あの１月２０日のたまたま観測したデータがなかったら、もっと長い周期の変光星と思ったことでしょう。
CK Aurという周期不明の変光星があり、それを試みたのですがそれについてはまた。
EASolverを少し改良
先日ブログに書いたEASolverを少し改良してみました。
devianceはいわばO-Cの自乗和なのですが、仮定している周期が長ければO-Cが当然大きい値になり得るし、仮定している周期を短くすれば小さくなるわけです。
極端な話、周期を0.001日位に仮定したら自乗和もほぼそれ以下になるわけですから、仮定している周期の自乗とO-Cの自乗和との比(P^2)/devを取って、その大小で周期を評価したほうが良いかと考えました。（大きいほど確からしい）
実際のスクリプトは下の方に書いておきます。（極小のデータの数、ｍが大きいと当然devも増えるわけで、(m*P^2)/devでそれぞれの周期を評価するのがもっと良いかもしれません。）
このように改良したスクリプトで、CW Oriのホフマイスターのデータを使って実際にやってみました。（この星の周期の改訂を2014年にVSXに届けました。）
これがホフマイスターの極小観測データです。
2425202.61 13.5
2425245.49 13.5
2425298.47 13.5
2425322.33 13.5
2425346.35 13.5
2425676.33 13.5
改良したスクリプトで周期を1－20日の範囲で探させました。
ベスト３の周期です。
ホフマイスターが出した周期は4.78505で、今でもGCVSにそう記載されており、彼がこの周期を正しいと思ったのも当然のようです。
私が観測データやASASデータなどから得た実際の周期は1.70993でした。
一発で正解、とはいかないようですが、このように候補をかなり絞れそうです。
下のグラフは改良前のスクリプトで作ったものです。1日－5日の範囲で探させました。
下はペランソのEASolverで延々と時間をかけて探させた場合です。
改良前のスクリプトの結果によく似ています。
３つの候補は出てきますが、順位が全く逆です。
改良したスクリプトです。ところどころ前のスクリプトと違う文字を使用していますので、全部載せます。
# mini: vector with m elements of minima(HJD) in ascending order
# Period to search between P1: minimum period & P2: maximum period
P1 <- 1
P2 <- 20
options(digits=12)
m <- length(mini)
d <- mini - mini[1]
d[m] # distance between the first & the last minimum in day
N1 <- floor(d[m] / P2)
N2 <- ceiling(d[m] / P1)
n <- N2 - N1
pow <- vector("numeric",n+1)
for (i in 1:(n+1)) {
P <- d[m] / (N1 + i - 1)
N <- round(d / P)
res <- lm(d ~ N)
dev <- deviance(res)
pow[i] <- (P^2) / dev
}
plot(d[m]/(N1+(1:(n+1))-1),pow, pch=16,col=3,grid(col=4),xlab="P")
max(pow)
N0 <- N1 + which.max(pow) - 1
N0 # N for the last minimum
P0 <- d[m] / N0
P0 # found Period
N <- round(d /P0)
res <- lm(mini ~ N)
coef(res) # improved epoch & Period
residuals(res) # residuals for each minimum
ord <- order(pow, decreasing=T)
P <- d[m] / (N1 + ord[1:3] -1)
P # the best 3 periods
ヨハネスケプラー
先週レーゲンスブルクに行きました。そこで何か見るものはないかと、旅行案内書を眺めていたら、なんとケプラーはこの街で1630年に亡くなっているのです。彼が住んで亡くなった家が当時のままに最近修復されてケプラー博物館となっていることがわかり尋ねていったのですが、残念ながら行った日は閉館日で中には入れませんでした。
ちょっと汚れているが建物にあった記念碑。
しまっているガラス扉越しに中を覗いてみた。
一方レーゲンスブルクの駅のすぐ前に緑地と小さな教会があります。
その教会の入り口脇にケプラーのことを書いた記念碑がありました。
「ここ、かつてのザンクトペーター墓地に、天文学者、宇宙の調和者、キリスト教統一運動の創始者ヨハネス・ケプラー（1571-1630）が彼の妻スザンネとともに眠っている。」
と書かれています。
彼が埋葬されたあとで戦争によって墓地があらされてしまい、どこに誰が埋葬されているか全くわからなくなってしまったそうです。そしてその墓地は今では緑地公園となっていて、私が訪れたときは一面凍った雪で覆われていました。
この雪の下のどこかにケプラーは眠っているわけですね。
KO GemとOQ Gem Mdyさんの観測
５日間の小旅行から昨晩戻りました。
その間にMdyさんがKO GemとOQ Gemの観測をしてくださいました。
ありがとうございました！
KO Gemは第二極小。ほぼ計算通りでした。先日のKisさんの第一極小を入れた位相図にピンク色でMdyさんの第二極小の観測を書き加えました。光度、フェーズ共にピッタリのようです！
OQ Gemの方、観測報告に書いておられましたが、
あとちょっと、というところでしたね！
まだ食の前のようでした。
今後もどうかよろしくお願いします！
以前書きましたが、今仮定している周期2.736..日で計算した第二極小の時に一度観測したのですが、全く変光が見られなかったので、ひょっとしたら本当の周期はこの２倍かもしれません。その場合は、今回のMdyさんが観測なさろうとした極小のように予報のnが奇数の極小が第二極小にあたるので、その観測が一層重要さを帯びます。（Mdyさんの場合はn=5でした。）日本で次のそのようなチャンスは
date UT n
14 15.2362 9
１４日となります。
KO Gem とOQ Gemの2月の極小予報
両星の要素がほぼわかったようなので、両星の2月の極小予報を書いておきます。
KO Gem。第二極小は以前に書きましたが、かなり浅いようですがあるはずです。
OQ Gen。こちらの第二極小は僕の観測ではかかりませんでした。
今晩スイスでなら第一極小を観測できるはずなのですが、雨で駄目です。11日後の2月11日にまたベストチャンスが来ます。
明日から日曜まで小旅行に出かけますので、ブログ書けないかもしれませんがあしからず。
OQ Gem その後
この星を晴れればしつこく観測していますが、極小は１月20日に観測できたっきりです。明けても暮れても（暮れても暮れてもが本当でしょうが）食外でほとんど変光しません。自作のEASlolverで見つけた周期で位相図を作ると、幸い食外のデータが極小付近に来ることはありません。（見つけた周期がデータと矛盾しないということ）しかし、フェーズ0.5のすぐそばのデータが減光を全く示していないのは気になります。
こういう場合はセバスチャンによれば周期が２倍かもしれません。（第一極小の減光幅が大変に大きい場合は第二極小がほとんど受からない可能性も多いようですが。）
２９日の朝にまた極小があるはずでした。スイスではもう夜が開けてしまっていますが、その極小への減光が受かるかもしれないと２８日夜半から観測しました。はっきりとしたことは言えないのですが、その減光の開始が受かった感じもします。
次の位相図で水色のデータがその観測データです。
減光開始に見えますか？
この2.7367..日という周期は２倍がほぼ5.5日なので11日おきにほぼ同じ時刻にまた極小が来ます。極小が受かった1月20日の11日後、明日３１日の夜はスイス時間で２３時５０分ごろに極小があり最高のチャンスなのですが、天気予報は雨！
KO Gem 続き
２６日にmei/nekoさんが観測してくださった極小を今までのSWASPデータ、ホフマイスターの唯一の極小などと組み合わせ位相図を作るとこのようになりました。mei/nekoさんの観測はIcフィルターによるものですが、SWASPデータにほとんどぴったり乗ります。
SWASPデータとmei/nekoさんの極小を重ねるとホフマイスターの極小が右にずれます。（mei/nekoさんの極小を元期にし、今までの周期2.8601607では短すぎるので少し長くして2.86017にします。）
ホフマイスターの極小がフェーズ0.0に来るようにするにはこれでは周期が長すぎで、2.8601607にしなければなりません。
そうするとSWASPデータが左に外れます。ホンの僅かですが、周期が変化したのかもしれません。
