eclminで極小時刻を求める
今日たまたま古いコンピューターを動かしたのでそこに残っているVSOLJの古いメールを眺めたら、eclminというのをKatさんが書いておられました。
以前はRの使い方が全くわからなかったのでほとんど読んでいなかったが、今なら使えるかもしれないと、試してみました。
見事に動きます。
eclmin(dat, center, width)
centerに大雑把な極小の時刻をいれ、幅は最初は少し大きめに取って動かすと、min(ecl)で極小時刻が出るのでそれを新たにcenterにいれてwidthを少し狭くするのが良いのかなと理解しています。
今月18日のOQ Gemの観測データを入れて極小時刻を求めたら、7803.3728(3)と出ました。ペランソでやったのは7803.3731と出ていました。
一回目はwidthに0.02を、二回目には0.01を入れてみました。
さーこれで、ペランソのお世話にならないでもすべてRで処理できます。
まだ試していないのですが、極小近くのデータに穴が空いていても動くのでしょうか？
追記
Ngaさんが求めてくださった極小時刻は7803.3729(3)でほとんどeclminで求めた値7803.3728(3)と同じでした！
PFilter 周期をふるいにかける
先週、OQ Gemの二回目の極小を観測するまでに幾つもの食外の観測がありました。（Mdyさんも何度か手伝って下さいました。）
そこで、PFinderが出した良さそうな周期を使って何十も位相図を作り、今まで観測した食外のデータが極小付近に来ないかどうかで周期をふるい分けました。かなり大変な仕事なので、それをスクリプトでできないか考えて、PFilterというのを作りました。
これは今までの食外の観測の開始時刻と終了時刻をリストにしたマトリックスと候補の周期の集合をインプットして、その周期集合をふるいにかけるものです。
細かいことはまた書きます。
結果を載せますと、まず周期集合
これをふるいにかけました。
ほとんど全滅で困り、相当のストレスでした。
そうこうしていたら１８日夜の観測で極小が引っかかったのです。
これでいろんな基礎データが大きく変わりました。まず極小データが一つ増え、周期の可能性が、僕の二回の極小の間隔28.8688日の1/1、1/1.5、1/2、1/2.5、、、1/14.5という２８個だけのの数列になりました。（周期が2日未満は不可とわかっています。）
PFinderと同じような考えでその中から良さそうな周期を選択するスクリプトPSelectorを作り、周期を選ばせたら上位から３つが
となりました。これは周期集合がたった28個なので瞬時に答えが出ます。
このスクリプトは手作業で周期7.2160日を見つけた後で作ったのですが、この結果をみてだいぶ安心しました。
その計算を載せておきます。
> div <- seq(1,14.5,0.5)
> P <- 28.8688/div
>
> n <- length(P)
> pow <- vector("numeric",n)
>
> for (i in 1:n) {
+ N <- round(d / P[i])
+ res <- lm(d ~ N)
+ P[i] <- coef(res)[2]
+ dev <- deviance(res)
+ pow[i] <- (P[i]^2)* m / dev
+ }
>
> max(pow)
[1] 281.268739107
> P0 <- P[which.max(pow)]
> N <- round(d /P0)
> res <- lm(d ~ N)
> coef(res) # epoch & Period
(Intercept) N
0.203144855046 7.216070079084 選ばれた周期
> residuals(res) # residuals for each minimum
1 2 3 4
-0.2031448550418 0.0144779481667 0.9367946040411 -0.1907731464772
5 6
-0.2809371170861 -0.2764174336027 残差です。
>
> ord <- order(pow, decreasing=T)
> P[ord[1:3]] # the best 3 periods
[1] 7.21607007908 4.12424539395 3.20775292472
やはり３番目の極小のずれが0.9367..と大きいです。
３番目の極小データを取り除いて同じことをさせるとPが少し減り7.2160056..となって手仕事で見つけた7.2160と同じになります。
> coef(res) # epoch & Period
(Intercept) N
0.11360649727 7.21600564715 選ばれた周期。
> residuals(res) # residuals for each minimum
1 2 3 4
-0.1136064972650 0.1098796119930 0.0358763680185 -0.0184634469843
5
-0.0136860357622 残差です。
>
> ord <- order(pow, decreasing=T)
> P[ord[1:3]] # the best 3 periods
[1] 7.21600564715 28.86955769417 3.60871702552
PFinder（EASolver改め）
これの続きです
EASolverと名付けましたが、ペランソのソフトにEASolverという名がもう付いているのと区別するためと、なかなか問題を解決(solve)してくれないので、「良さそうな周期の候補を評価付きで見つける」という感じでPFinderにしました。スクリプトも少し不必要なところを削ったりしましたが、中身がかわったわけではありません。
ところで、このPFinderでOQ Gemの極小データ（前回書いたようにHoffmeister、Frank、Kasaiのそれぞれ2回づつの計６つのデータ）を使って周期を探させても今回見つかった周期7.2160日と言うのは下のリストに出てこないのです。がっくり！
このトップのいくつかの周期で位相図を先週末作ってみたのですが僕の２回の極小があまりにも重ならない。例えば
結果論的ですが、原因は昨日の位相図でフェーズ0.0にも0.5からも外れているFrankの極小１つのためです。下図での0.0から右に外れている緑の三角。
これは下のリスト中３番めの極小データです。
これをリストから取り除いてPFinderで周期を探させると
7.2160がちゃんと第３位で登場します！残念ながら１つのデータのせいで、見つかるものも見つからなくなってしまいますね。
OQ Gemの周期を探そう
極小が２回受かったのは大変に幸運で良いのですが、周期をどう探すか考えてあれこれスクリプトを書いたりしてだいぶ苦労しました。でも最終的にはクラシカルな方法で言わば手作業でやったのが確実なようでした。
現在のところ、Hoffmeisterの極小が２回、Frankのが２回、僕のが２回と合計６回の極小観測があります。これがそのリストです。
これらの極小時刻の間隔を求めます。
HoffmeisterとFrankの間が相当空いています。
これらの間隔に僕の極小観測の間隔がいくつ入るか、28.8688で割ります。
５番目は当然１になるはずの数値です。
他のをよく見てみましょう。１番目は22.75回、２番めは484.75回、３番めは24.5回、４番目は139回と思われます。
これらをベクトルにします。
少数点以下に0.75があるのは周期が28.8688では長すぎで、その半分ないし1/４でなければならないことを意味します。ｈを２倍してｇを割ると各部分での周期が出ます。
これを見ると周期が一旦短くなってまた伸びたように見えます。
周期14.43日付近で極小データがどのように位相図にのるか試した結果、14.4320が良さそうとなりました。元期は１８日夜の極小時刻です。赤丸がHoffmeisterで第一極小と第二極小にきれいに載っています。緑の三角はFrankの極小で１つ少し外れます。僕の極小は二つとも当然0.0に載っています。
本当の周期がこの14.4320日か、その半分の7.2160日かは後者の周期で計算して第二極小に当たる時に減光が見られるかどうかが決め手になるかと思います。
前者の周期で作成した位相図です。黄色が18日の観測で、その下に1月20日の極小データ（緑）がありますが重なってほとんど見えません。
後者の周期で作成した位相図です。
この位相図の要素で予報を作りました。
なお来る２月２２日は第二極小なので変光範囲は小さいかもしれません。上記のように、全く変光がなかったら周期は２倍の14.4320が適当と決められるでしょう。あるいはnが奇数の極小データが偶数の極小データと異なっても周期は２倍と確実にいえますね。
とにかく今度は当たりますように！外れたら「狼少年」になってしまって誰も信じてくれなくなってしまう！
２月分
３月分
OQ Gemの極小を再び観測
昨晩は霧が途中で発生することもなくずーっと晴れてくれ、食をほぼ始めから終わりまで観測できました。
この極小や他の極小観測から周期を求めようとしていますが一筋縄では行かない感じです。私の２回の極小観測が重なるように周期を探すと、4.819427日というのが出ましたが、これだとHoffmeisterやFrankの極小観測が0.0にも0.5にもほとんどのりません。
OQ Gem ただ今減光を観測中！
今日の午後またこの食変光星の周期について色々試行錯誤して、8.3日ぐらいというのがみつかり、計算したら今夜副極小があると出ました。今のところきれいに晴れていて観測中で、どんどん減光しています！当たりました！
詳しくはまた明日。
OQ Gem その後
これの続きです。
あの時点ではとりあえずの結論として、周期7.06....が良さそうと書きましたが、よくよく見るとダメでした。1月6日の平坦な観測が復光直前に引っかかります。
７日の周期を考えた時に１月２０日に極小（緑）があって、その２週間前の１月６日の最初の観測（赤）に変光が全く引っかかっていないのはおかしいと思いついたのです。
今、EASolverの続きのようなスクリプトを製作中で、それでもこの7.06...はダメと出たのです。
