スクリプトでミスが
この記事の中のスクリプトで最後の行にカッコのミスがありましたので訂正します。
list(coef(res),residuals(res),P[order[1:3]])となっていますが
このorderの後1:3のカッコは()でなければなりません。
P[order(1:3)]が正しいです。すみません。
カッコは小さくてカギカッコなのか丸いのかよく見えない！
AC CMaの食
２月２６日夜に、日本でMdyさんがこの星の食を観測してくださいました！
１９４２年来、食は初観測でしょうか？この星は南にあって観測できる時間は短めで高度は低いところを観測してくださってありがとうございます！
以前ブログに載せたASASデータで作った位相図を見ると食の底はASASには写っていないので、ひょっとしたら皆既食が続いているのかと思っていましたら、そうでした。
ASASの位相図にMdyさんの観測を重ねた位相図を作りました。
すごい減光です。食外と皆既食の光度差は２等級以上あります。今これを作っていて何か間違っていないかと思ったのですが、VSXにも写真等級で12.8-15とありますから、これであっているのでしょう。Mdyさんの観測では、30分で約0.5等暗くなってます。
AC CMa の３月の極小予報
日本だと５日、１５日、２２日あたりが主極小のチャンスでしょうか？
KO Gemを観測しています
mei/nekoさんが５２年ぶりの極小を観測してくださった星です。ようやくスイスでも極小を観測できるチャンスが来ました。今日は雲もなく昨晩より良く安定した観測ができそうです。極小時刻はスイス時間で２２時２０分頃のはずですが、夕方に観測し始めたときにはすでに減光し始めていて、どんどん暗くなってます。
OQ Gem減光確認できず
昨晩２６時過ぎに極小があるはずで、あれば周期が7.2160日が確定できたのですが、減光が始まる前に完全に曇られました。残念です。。
雲の隙間から一枚でも二枚でも写らないかと１時半まで粘ってみましたが全くチャンスなし。オートガイドを切ってその後は寝ました。今朝撮影したフレームを見たけど一枚も変光星は写っていなかった。（明るい星は少しは写っているけど、１４等台は駄目です。）次のチャンスは３月末。。。
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の周期を探してみました。（この星について昨年８月に調査してブログにもすでに少し書きました。）
これが今までわかっている４回の極小観測時刻です。
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]])
}
