pfilter 食外の観測データを使って周期を絞る

周期が未知の変光星を観測しても大体は「変光なし」のデータになってしまいます。今回のOQ Gemでもかなりの食外のデータが得られました。（下参照）それでもある程度食外のデータがたまるとそれを使って周期の可能性を絞ることができます。昨日記載したスクリプトのpfilterがそのためのものです。

スクリプトを読んでお分かりいただけた方も多いと思いますが、一応どう考えたか説明してみます。

変光星の１つの極小が観測されたとしてそれを元期とします。そして今ある周期Pを仮定し、時間軸上で元期を原点として周期２Pのサインカーブを考えます。一方で食外観測の開始時刻と終了時刻をそれぞれo1とo2するとこの観測期間中に極小が来ないための必要十分条件はサインカーブのゼロ点がこの期間中に来ないことです。すると、これはo1とo2の時点でのsinの値が同じ符号(++か--)であることと同じです（ゼロではいけない。Pは当然観測期間より長いのを仮定する）。ただし、今は極小だけを考えましたが、食外の観測ならば減光も増光も受からなかったということなので、食の継続時間の半分をo1から引きo2の方には加えてo1o2の代わりにし期間を前後に増やします。

このような考え方をスクリプトにしました。

サインまで持ち出すのは大げさかとも思いました。例えばフェーズの値でも同じようなことをできそうですが、まあスクリプトのステップとしてはほとんど同じくらいなので思いついたのを書きました。

実際に見てみましょう。

食外の観測の開始、終了時刻です。１２回食外の観測がありました。

Mdyさんのが３回入っています。

2457760.473092 14.177 0.0000
2457760.689499 14.252 0.0000
#
2457775.259450 14.203 0.0000
2457775.619496 14.145 0.0000
#
2457776.302090 14.128 0.0000
2457776.527811 14.174 0.0000
#
2457780.381774 14.192 0.0000
2457780.561295 14.232 0.0000
#
2457781.230634 14.162 0.0000
2457781.369587 14.237 0.0000
#
2457782.480863 14.183 0.0000
2457782.619896 14.275 0.0000
#
2457783.242855 14.289 0.0000
2457783.450784 14.184 0.0000
#
2457788.035484 14.230 0.0000
2457788.072142 14.190 0.0000
#
2457799.006189 14.299 0.0000
2457799.228353 14.470 0.0000
#
2457799.235402 2.611 0.0900
2457799.268080 2.600 0.0310
#
2457800.042613 14.234 0.0000
2457800.184863 14.227 0.0000
#
2457802.464942 2.538 0.0180
2457802.513308 2.550 0.0710

これらのデータをdatに入れ、pfindで見つけた周期のセットをふるいにかけます。

> result <- pfind(mini, 2,20)
> result3　　　　　　　　　pfindで見つけた周期のセット

[1] 2.21904683973 3.20447997799 7.21600572527 3.21991626311 5.81441792209
[6] 2.06590495001 7.13893824660 5.76427706585 2.22643808719 3.58863255498
[11] 3.60800286263 2.21170450419 2.05953961682 2.07230975125 7.29475530380

> pfilter(result3,dat,2457803.3727,0.13)　　0.13が食継続時間の半分。
[1] 7.21600572527 5.81441792209 3.58863255498 3.60800286263

結果として、result3の最初の２つの周期はふるいにかけたらダメとなって、３番目の周期が通りました。めでたしめでたし。

ところで、これを書きながらOQ Gemを観測しています。今夜半過ぎに極小があるはずで、薄曇りですが今のところ観測できる状態です。果たして減光が受かるか？

極小時にはおそらく星没でしょう。

周期を見つけるための３つのスクリプト

先日から書いていますが、OQ Gemの周期を探すのにスクリプトを３つpfind  pfilter  pselctを書いてみました。

全部functionの形にしたのでconsoleから簡単にまとめて呼び出せます。一応下の方に載せておきます。（まだまだ改良出来るかと思います。ご意見がおありでしたらどうかお知らせください。）

今回はもう今までにもブログに書いたPFinderをさらに改名したpfindについて。

pfindはいくつかの既知の極小時刻から、周期をある決まった範囲で探すものです。これで例えばGCVSでもVSXでも周期未知のV1774 Ophの周期を探してみました。（この星について昨年８月に調査してブログにもすでに少し書きました。）

muttenz.hatenablog.com

これが今までわかっている４回の極小観測時刻です。

JD mag from BAN 1968
2435638.835 14.0
2435960.916 14.0
2436728.781 14.0
2436762.818 14.0

この４回のデータをpfindにいれて周期１日から２０日の範囲で探させました。

3.09635日が有力候補として出ました。昨年ASASのデータなどからPDMで得た周期は

3.096745でした。たった４回のdiscreteなデータから周期の候補が得られるのはたいへんありがたいです。（ASAS、KWS、NSVS、SWASPなどのデータからこれらの周期で位相図を作成してみて、どの周期が良さそうか判断することができることがあります。いつもとは行かないのですが。。。）

明日はpfilterについて書きます。

# to find possible periods
# mini: vector with m minima(HJD) in ascending order
# P1: minimum period & P2: maximum period

pfind <- function(mini, P1, P2) {

options(digits=12)
m <- length(mini)
d <- mini - mini[1]

N1 <- floor(d[m] / P2)
N2 <- ceiling(d[m] / P1)
n <- N2 - N1
P <- vector("numeric",n+1)
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)
P[i] <- coef(res)[2]
dev <- deviance(res)
pow[i] <- (P[i] ^ 2) * m / dev
}

plot(d[m]/(N1+(1:(n+1))-1),pow, pch=16,col=3,grid(col=4),xlab="P")

P0 <- P[which.max(pow)] # the best found period
N <- round(d /P0)
res <- lm(mini ~ N)
ord <- order(pow, decreasing=T)
list(coef(res),residuals(res),P[ord[1:15]])
# improved epoch & period
# residuals for each minimum
# the best 15 found periods
}

# to determin if a suposed period is
# acceptable(TRUE) or not(FALSE)
# P: a set of n periods
# obs: data of starts(hjd) and ends(hjd) of
# m observations without variation
# hed: half of eclipse duration (day)

pfilter <- function(P, obs, epo, hed){
pi <- 3.1415926535

OBS <- matrix(dat = obs$V1, byrow=T,nc=2)
OBS[,1] <- OBS[,1] - hed
OBS[,2] <- OBS[,2] + hed

n <- length(P)
det <- vector(length = n)

for (i in 1:n) {
p <- P[i]
ph <- sin(pi * (OBS - epo) / p)
res <- ph[,1] * ph[,2]
det[i] <- all(0 < res)
}

result <- data.frame(P,det)
true <- subset(result, det=="TRUE")
return(true$P)
}

# to select suitable periods from a set of periods
# mini: vector with m minima(HJD) in ascending order
# P: a set of n periods

pselect <- function(P, mini){
options(digits=12)

m <- length(mini)
d <- mini - mini[1]

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
}

P0 <- P[which.max(pow)]
N <- round(d /P0)
res <- lm(d ~ N)
coef(res) # epoch & Period
residuals(res) # residuals for each minimum
ord <- order(pow, decreasing=T)
list(coef(res),residuals(res),P[order[1:3]])
}

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改め）

muttenz.hatenablog.com

これの続きです

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にもほとんどのりません。