読者です 読者をやめる 読者になる 読者になる

muttenz's blog

スイス星空だより

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

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

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

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

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

muttenz.hatenablog.com

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

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

この4回のデータをpfindにいれて周期1日から20日の範囲で探させました。

f:id:muttenz:20170225053214j:plain

f:id:muttenz:20170225053552j:plain

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

3.096745でした。たった4回の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]])
}