VSOLJのデータからAB Aurの光度曲線を描いてみたのですが、ばらつきもあってよくわかりませんでした。
そこで、10日毎くらいの平均光度曲線を描いてみようと思ったのですが、googleっても、離散的で欠測値も多いデータの期間を区切っての平均の求め方はよくわかりませんでした。
ならばと、Katさんのphaveを読んで以下のようなのをでっち上げたました。
エラーバーの計算もphaveをパクリです。
#brの計算、もう少し真面目のやったほうが良かったかも？
apply(tapply)って便利ですよね。ループを記述しなくていいの。
d <- read.csv("aurab.jd", header=F)
br <- (max(d$V1)-min(d$V1))/10
bin <- factor(cut(d$V1, breaks=br, label=NULL))
jd <- as.vector(tapply(d$V1, bin, mean))
mag <- as.vector(tapply(d$V2, bin, mean))
num <- as.vector(tapply(d$V2, bin, length))
err <- as.vector(tapply(d$V2, bin, sd))
err <- ifelse(num > 1, err/sqrt(num-1), NA)
plot(jd, mag, xlim=c(55000, 56000), ylim=c(8, 6), xlab="jd-2400000", ylab="mag.")
arrows(jd, mag+err, jd, mag-err, length=0.05, angle=90, code=3)
ちなみに、平均光度曲線の黒丸は、Iakさんのデジカメ測光とわたしのCCDでの観測を後から上書きしたものです。
ちゃんと操作内容を理解していないから、嘘やってたら、ごめんなさい。
結果をみて、よさそうだったので、それ以上追求していません。
#errがNAの時に、arrowsが、「長さが0だ」と警告を返すので、本当は、ちゃんとエラー処理をしないといけません。
