meineko’s blog

元つくばの某独立行政法人勤務の植物屋です。最近は、ほぼ、突発天体の話題です。

Rで周期解析(自動解析編)

天文 R

Katさんが、htmlを自動生成してくれるプログラムを書いてくれたので、テストしました。
VSOLJの観測データを元にしたミラ型変光星の光度曲線集です。
https://meineko.com/pgm/mira_pgm_test.html
なお、自動解析(探索範囲が固定)なので、おかしな結果のものも混ざっていますが、ほとんど、良好です。
光度曲線の形はいくつかのタイプに分けられそうですね。周期と形の関係はどうでしょう？周期以外の要因は？スペクトル型とか？


以下の様なm.typeというファイルを作っておいて

M
M*

liststar | seltype m.typ > mira.lst

を実行すると、vsolj.dbfにデータのあるミラ型変光星のリスト(mira.lst)ができます。
なお、実行には、vartype.dbfが要ります。
で、以下を実行すると

for (i in 1:nrow(obj)) {
    star <- obj$V1[i]
    d <- readdb(star)
    if (!is.null(d) && nrow(d) > 1000) {
        pergrmhtml(star)
    }
}

リストを元に、次々と周期解析をして、png画像とhtmlファイルを作ってくれます。
なお、データ数が1000より多い星に絞って解析をしています。


あらかじめリストを作っておくのは、効率化のためと、データのない星を呼ぶとエラーで止まるからですが、エラー処理をうまくやれば、他の方法で作った星のリストを食わせるだけでもよいかもしれません。
#データのない星ばかりだと無駄に時間がかかりますが。

なお、pergrmhtmlは、

pergrmhtml <- function(star,start=100,end=500,step=2000) {
    d <- readdb(star)
    pgm <- pergrm(d,start,end,step)
    p <- max(pgm)
    f1 <- sprintf("%s_per.png",star)
    f2 <- sprintf("%s_phase.png",star)
    fout <- sprintf("%s.html",star)
    png(f1)
    plot(pgm)
    dev.off()
    png(f2)
    plot(phave(d,0,p,100))
    dev.off()
    cat("<html>\n",file=fout)
    cat("<h2>\n",file=fout)
    cat("<center>\n",file=fout,append=TRUE)
    cat(sprintf("<img src=%s>\n<br>\n",f1),file=fout,append=TRUE)
    cat(sprintf("Period=%.1f\n<br>\n",p),file=fout,append=TRUE)
    cat(sprintf("<img src=%s>\n<br>\n",f2),file=fout,append=TRUE)
    cat("</center>\n",file=fout,append=TRUE)
    cat("</h2>\n",file=fout,append=TRUE)
    cat("</html>\n",file=fout,append=TRUE)
}

です。


次々と改良版を出されるので、fix出来ないのですが、以下、本日段階でのバージョン、自分用のメモです。
jdfil.awkの改良版

BEGIN {
}

{
    if (substr($0,1,1) != "#" && index($2,"<") == 0) {
        mag = $2;
        gsub("[<A-Za-z:]","",mag)
        sys = $2;
        gsub("[<0-9.:/-]","",sys)
        if (sys == "")
            sys = "-";
        if (index(mag,".") == 0)
            mag = sprintf("%.1f",mag*0.1)
        printf("%s %s %s %s\n",$1,mag,sys,$3)
    }
}

lfit <- function(dat) {
  res <- lm(dat$V2 ~ dat$V1)
  dat$V3 <- dat$V2
  dat$V2 <- res$residuals
  return(dat)
}


phase <- function(dat,epoch,period)
{
    if (class(period) == "minper") {
        p <- period$minper[1]
    }
    else {
        p <- period
    }
    f <- dat
    f$mag <- dat$V2
    f$ph <- ((dat$V1 - epoch) %% p)/p
    if (!is.null(dat$err)) {
        f$err <- dat$err
    }
    f1 <- subset(f,f$ph < 0.5)
    f1$ph <- f1$ph + 1
    f2 <- subset(f,f$ph >= 0.5)
    f2$ph <- f2$ph - 1
    f <- rbind(f2,f,f1)
    class(f) <- "phase"
    return(f)
}

phave <- function(dat,epoch,period,ndiv=60)
{
    if (class(period) == "minper") {
        p <- period$minper[1]
    }
    else {
        p <- period
    }
    dat$ph <- ((dat$V1 - epoch) %% p)/p
    br <- seq(0,1,length=ndiv+1)
    dat$bin <- factor(cut(dat$ph,breaks=br,labels=FALSE))
    ph <- as.vector(tapply(dat$ph,dat$bin,mean))
    num <- as.vector(tapply(dat$V2,dat$bin,length))
    mag <- as.vector(tapply(dat$V2,dat$bin,mean))
    err <- as.vector(tapply(dat$V2,dat$bin,sd))
    err <- ifelse(num > 1, err/sqrt(num-1), NA)
    f <- data.frame(ph,mag,err,num)
    f1 <- subset(f,f$ph < 0.5)
    f1$ph <- f1$ph + 1
    f2 <- subset(f,f$ph >= 0.5)
    f2$ph <- f2$ph - 1
    f <- rbind(f2,f,f1)
    row.names(f) <- 1:nrow(f)
    class(f) <- c("phave",class(f))
    return(f)
}

plot.phave <- function(a,xlim=NULL,ylim=NULL,xlab="Phase",ylab="Mag",axes=TRUE)
{
    if (is.null(ylim)) {
        ymax <- max(a$mag)
        ymin <- min(a$mag)
        y1 <- ymin - (ymax-ymin)*0.1
        y2 <- ymax + (ymax-ymin)*0.1
        yl <- c(y2,y1)
    }
    else {
        yl <- ylim
    }
    if (is.null(xlim)) {
        xlim <- c(-0.5,1.5)
    }
    plot(x=a$ph, y=a$mag, xlim = xlim, ylim = yl, pch = 19, xlab = xlab, ylab = ylab,axes=axes)
    arrows(a$ph, a$mag-a$err, a$ph, a$mag+a$err, length=.05,angle=90,code=3)
}

phave.save <- function(a,outfile)
{
    fout <- file(outfile,"w")
    for (i in 1:nrow(a)) {
        s <- sprintf("%8.5f %8.5f %8.5f %4d\n",
                     a$ph[i],a$mag[i],a$err[i],a$num[i])
        cat(file=fout,s)
    }
    close(fout)
}

yerr.plot <- function(x,y,e,ylim=NULL,xlab="x",ylab="y")
{
    plot(x=x, y=y, pch = 19, ylim=ylim, xlab = xlab, ylab = ylab)
    arrows(x, y-e, x, y+e, length=.05,angle=90,code=3)
}

plot.phase <- function(a,xlim=NULL,ylim=NULL,xlab="Phase",ylab="Mag",axes=TRUE)
{
    if (is.null(ylim)) {
        ymax <- max(a$mag)
        ymin <- min(a$mag)
        y1 <- ymin - (ymax-ymin)*0.1
        y2 <- ymax + (ymax-ymin)*0.1
        yl <- c(y2,y1)
    }
    else {
        yl <- ylim
    }
    if (is.null(xlim)) {
        xlim <- c(-0.5,1.5)
    }
    plot(x=a$ph, y=a$mag, xlim = xlim, ylim = yl, pch = 19, xlab = xlab, ylab = ylab,axes=axes)
}

phave.points <- function(a,xlab="Phase",ylab="Mag")
{
    points(x=a$ph, y=a$mag, pch = 19, xlab = xlab, ylab = ylab)
    arrows(a$ph, a$mag-a$err, a$ph, a$mag+a$err, length=.05,angle=90,code=3)
}


#readdb <- function(obj,start="19000101",end="20991231") {
#    cmd <- sprintf("vcut %s %s-%s | pack | tojd | awk -f ~/bin/jdfil.awk",obj,start,end)
#    read.table(pipe(cmd))
#}


readdb <- function(obj,start="19000101",end="20991231") {
    cmd <- sprintf("vcut %s %s-%s | pack | tojd | gawk -f ./jdfil.awk > temp.txt",obj,start,end)
    system(cmd)
    read.table("temp.txt")
}

pergrm <- function(d,start,end,div) {
    d$V1 <- d$V1 - mean(d$V1)
    d$V2 <- d$V2 - mean(d$V2)
    p <- seq(start,end,length=div+1)
    pw <- as.numeric(lapply(p, function(tp) {
        ph <- ((d$V1/tp) %% 1)*pi*2
        (sum(d$V2*sin(ph)))^2 + (sum(d$V2*cos(ph)))^2}))
    r <- data.frame(period=p,power=pw)
    class(r) <- c("Pergrm",class(r))
    return(r)
}

plot.Pergrm <- function(pgm,ylim=NULL,xlab="Period",ylab="Power",...) {
    if (is.null(ylim)) {
        ymax <- max(pgm$power)
        ymin <- min(pgm$power)
        y1 <- ymin - (ymax-ymin)*0.05
        y2 <- max(1,ymax + (ymax-ymin)*0.05)
        yl <- c(y1,y2)
    }
    plot(x=pgm$period, y=pgm$power, xaxs="i", ylim = yl, typ="l",
         xlab = xlab, ylab = ylab, ...)
}

max.Pergrm <- function(pgm,...) {
    n <- which.max(pgm$power)
    if (n <= 3 || n >= nrow(pgm)-2) {
        return(pgm$period[n])
    } else {
        a <- pgm[(n-2):(n+2),]
        ma <- mean(a$period)
        a$period <- a$period - ma
        m <- lm(a$power ~ poly(a$period,2,raw=TRUE))
        ce <- as.vector(coef(m))
        a <- ce[3]
        b <- ce[2]
        c <- ce[1]
        return(ma-b/(2*a))
    }
}

min.Pergrm <- function(pgm,...) {
    max.Pergrm(pgm)
}
    
fouriercomp <- function(d,p) {
    ph <- ((d$V1/p) %% 1)*pi*2
    ss <- sin(ph)
    cc <- cos(ph)
    m <- lm(d$V2 ~ ss + cc)
    a <- m$coefficients["ss"]
    b <- m$coefficients["cc"]
    return(c(a,b))
}

subfourier <- function(d,p,comp,gain) {
    ph <- ((d$V1/p) %% 1)*pi*2
    ss <- sin(ph)
    cc <- cos(ph)
    d$V2 <- d$V2 - (ss*comp[1] + cc*comp[2])*gain
    return(d)
}

clean <- function(d,start,end,div,step,gain=0.3) {
    px <- numeric(step)
    py <- numeric(step)
    for (i in 1:step) {
        pgm <- pergrm(d,start,end,div)
        p <- max(pgm)
        comp <- fouriercomp(d,p)
        d <- subfourier(d,p,comp,gain)
        px[i] <- p
        py[i] <- sqrt(sum(comp2))*gain
    }
    return(list(px=px,py=py))
}

plot.Clean <- function(cl,bw="nrd0",...) {
    wt <- cl$py/sum(cl$py)
    plot(density(cl$px,weights = wt,bw=bw),...)
}

pergrmhtml <- function(star,start=100,end=500,step=2000) {
    d <- readdb(star)
    pgm <- pergrm(d,start,end,step)
    p <- max(pgm)
    f1 <- sprintf("%s_per.png",star)
    f2 <- sprintf("%s_phase.png",star)
    fout <- sprintf("%s.html",star)
    png(f1)
    plot(pgm)
    dev.off()
    png(f2)
    plot(phave(d,0,p,100))
    dev.off()
    cat("<html>\n",file=fout)
    cat("<h2>\n",file=fout)
    cat("<center>\n",file=fout,append=TRUE)
    cat(sprintf("<img src=%s>\n<br>\n",f1),file=fout,append=TRUE)
    cat(sprintf("Period=%.1f\n<br>\n",p),file=fout,append=TRUE)
    cat(sprintf("<img src=%s>\n<br>\n",f2),file=fout,append=TRUE)
    cat("</center>\n",file=fout,append=TRUE)
    cat("</center>\n",file=fout,append=TRUE)
}