星空研究Memo

ここは某天文屋の外部記憶装置である。

JD形式からSTD形式へ変換

2013-08-07 00:01:50 | R言語

Muniwin の測光結果は .txt で吐き出すと時間が JD となっている。
これを VSOLJ-obs に報告しようとなると、時間を JST に変換し
STD 形式にしなければならない。
そこで、Muniwin の測光結果を R で STD 形式に変換するコードを書いてみた。

[追 記 & 訂 正]
某ML でコードの間違いの指摘や改良のアドバイスを頂きました。
それに伴い以下には更新版を掲載します。
# if を使わずにという宿題を頂きましたが... うーむ、結局 ifelse() を使ってしまいました。

==========================================
# Muniwin JD format -> VSOLJ STD one
in_file <- "BOOUZ_20130802.txt"
out_file<- "BOOUZ_20130802.std"
star <- "BOOUZ"
obs <- "Iak"
comp_mag <- 12.14
filter <- "C"
#-------(ここから下は変更しない)-------
#------------------------------
data <- read.table(in_file, skip=2)
 mag <- sprintf("%5.3f%s", data$V2+comp_mag, filter)
# JD -> STD
JD <- data$V1
jd <- JD + 0.5
z  <- floor(jd)
f  <- jd-z
aa <- floor((z-1867216.25)/36524.25)
a  <- floor(z+1+aa-floor((aa/4)))
b  <- a+1524
c  <- floor((b-122.1)/365.25)
k  <- floor(365.25*c)
e  <- floor((b-k)/30.6001)

 D <- floor(b-k-floor(30.6001*e))
 M <- ifelse(e<13.5, M1 <- e-1, M2 <- e-13)
 Y <- ifelse(M>2.5,  Y1 <- c-4716, Y2 <- c-4715)
  hh  <- floor(f*24)
  mm  <- floor((f*24-hh)*60)
  ss  <- round((((f*24-hh)*60)-mm)*60, digit=0)
     mm <- ifelse(ss=="60", mm1<-mm+1, mm2<-mm)
     ss <- ifelse(ss=="60", ss1<-0, ss2<-ss)
h9 <- hh + 9

nM  <- nchar(M)
MMM <- ifelse (nM == "1",
        MM <- paste("0", M, sep=""),M)
nD  <- nchar(D)
DDD <- ifelse (nD == "1",
        DD <- paste("0", D, sep=""),D)
nmm <- nchar(mm)
mmm <- ifelse (nmm == "1",
        m <- paste("0", mm, sep=""),mm)
nss <- nchar(ss)
sss <- ifelse (nss == "1",
        s <- paste("0", ss, sep=""),ss)
JST <- paste(Y, MMM, DDD, h9, mmm, sss, sep="")
#------- .std で吐き出す
 STD <- data.frame(star,JST,mag,obs)
 write.table(STD, file=out_file, row.names=F,col.names=F,quote=F)

最新の画像もっと見る

post a comment