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

muttenz's blog

スイス星空だより

EASolverを少し改良

先日ブログに書いたEASolverを少し改良してみました。

devianceはいわばO-Cの自乗和なのですが、仮定している周期が長ければO-Cが当然大きい値になり得るし、仮定している周期を短くすれば小さくなるわけです。

極端な話、周期を0.001日位に仮定したら自乗和もほぼそれ以下になるわけですから、仮定している周期の自乗とO-Cの自乗和との比(P^2)/devを取って、その大小で周期を評価したほうが良いかと考えました。(大きいほど確からしい)

実際のスクリプトは下の方に書いておきます。(極小のデータの数、mが大きいと当然devも増えるわけで、(m*P^2)/devでそれぞれの周期を評価するのがもっと良いかもしれません。)

このように改良したスクリプトで、CW Oriのホフマイスターのデータを使って実際にやってみました。(この星の周期の改訂を2014年にVSXに届けました。)

これがホフマイスターの極小観測データです。

2425202.61 13.5
2425245.49 13.5
2425298.47 13.5
2425322.33 13.5
2425346.35 13.5
2425676.33 13.5

改良したスクリプトで周期を1-20日の範囲で探させました。

f:id:muttenz:20170208180708p:plain

ベスト3の周期です。

f:id:muttenz:20170208181459j:plain

ホフマイスターが出した周期は4.78505で、今でもGCVSにそう記載されており、彼がこの周期を正しいと思ったのも当然のようです。

私が観測データやASASデータなどから得た実際の周期は1.70993でした。

一発で正解、とはいかないようですが、このように候補をかなり絞れそうです。

下のグラフは改良前のスクリプトで作ったものです。1日-5日の範囲で探させました。

f:id:muttenz:20170208185603p:plain

下はペランソのEASolverで延々と時間をかけて探させた場合です。

f:id:muttenz:20170208181621j:plain

改良前のスクリプトの結果によく似ています。

3つの候補は出てきますが、順位が全く逆です。

 

改良したスクリプトです。ところどころ前のスクリプトと違う文字を使用していますので、全部載せます。

# mini: vector with m elements of minima(HJD) in ascending order
# Period to search between P1: minimum period & P2: maximum period


P1 <- 1
P2 <- 20
options(digits=12)
m <- length(mini)
d <- mini - mini[1]
d[m] # distance between the first & the last minimum in day

N1 <- floor(d[m] / P2)
N2 <- ceiling(d[m] / P1)
n <- N2 - N1
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)
dev <- deviance(res)
pow[i] <- (P^2) / dev
}

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

max(pow)
N0 <- N1 + which.max(pow) - 1
N0 # N for the last minimum
P0 <- d[m] / N0
P0 # found Period
N <- round(d /P0)
res <- lm(mini ~ N)
coef(res) # improved epoch & Period
residuals(res) # residuals for each minimum

ord <- order(pow, decreasing=T)
P <- d[m] / (N1 + ord[1:3] -1)
P # the best 3 periods