周期を見つけるための３つのスクリプト
先日から書いていますが、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]])
}
