EASolver 食変光星の周期を見つける

有名な変光星のデータ解析のソフト「Peranso(Period analysis software)」にはEASolverというソフトが入っていて、食変光星、中でもEAのようなタイプの極小が幾つか観測されながら周期が未知の場合に、周期を探すことができます。

昨日書いたように、これで計算するとしばしば大変に長い時間がかかります。（OQ Gemで周期１日から７日まで探させたら３０分）しかも精度もイマイチの感じです。

昨日思い立って、Rでこの目的のスクリプトを書けないだろうかと試行錯誤してみました。当初は周期を少しずつ変化させて、すべての極小時刻をうまく1つの周期で結びつけることを試みたのですが、周期を荒く変化させたのでは精度が落ちるし、ひょっとしたら本当の周期をちょうど飛び越えてしまう危険もあります。だからと言って、周期を非常に細かいステップで変えていくと、ペランソと同じようにいつまでたっても結果が出ません。

そこで思いついたのは極小時刻を早い順に並べて最初と最後の時刻の間隔d[m]をある整数nで割って仮の周期Pを出し、その周期で途中にある極小がそれぞれほぼ何番目Nの極小に当たるかを出すようにしました。そしてその計算上の極小時刻と実際の極小時刻の差、いわばO-Cの自乗和devを出して、この自乗和devが最小となる整数を見つけるようにしてみました。これだといわば飛び飛びの数値で計算をさせる感じです。

昨日書いたように、これでOQ Gemの今までの極小観測を使って周期を出したら2.736...というのが見つかりました。

周期2日から20日の間で探させた場合、約13秒で出てきます。

Interceptのしたの左の数値が元期、右の数値が周期となります。

そのさらにしたにあるのがその要素で計算した極小時刻Cと実際のデータOのO-Cになります。

一応作ってみたスクリプトを載せてみます。

スクリプトを書くのはあまりうまくないのですが、そこはご容赦ください。

どうしてかdev[1]が０になってしまいます。そこを避けて通るのにちょっとエレガントではない行があります。

何か改良やこれはダメとかありましたらどうかご教示ください。

これで見つかった周期が確実とはならないことも、いつか周期を改訂したCW Oriのケースでわかりました。（dev0の値がホンのわずか大きいのが本当の周期でした。）なのでとりあえず候補を３つまであげられるように最後の幾つかの行を付け加えました。

ベクトルminiやP1、P2はそれぞれのスクリプトを動かす前に入れるところです。

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

P1 <- 2
P2 <- 20

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

N1 <- floor(d[m] / P2)
N2 <- ceiling(d[m] / P1)
k <- N2 - N1
dev <- vector("numeric",k+1) # deviance

for (n in 1:k+1) {
P <- d[m] / (N1 + n - 1)
N <- round(d / P)
res <- lm(d ~ N)
dev[n] <- deviance(res)
}

dev0 <- dev[-1]
min(dev0)
N0 <- N1 + which.min(dev0)
N0 # N for the last minimum
P0 <- d[m] / N0 # found Period

N <- round(d /P0)
res <- lm(mini ~ N)
coef(res) # improved epoch & Period
residuals(res) # residuals for each minimum

ord <- order(dev0)
P <- d[m] / (ord[1:3] + N1)
P # the best 3 periods