JD形式からSTD形式へ変換

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)

R: 最近の作図

R でグラフを描くさい、見栄えなどにこだわりだすと
気が付けばもの凄い時間が経っていたりする（ｗ
出来上がったコードはコピペすれば色々な局面で応用が利くので
最近使っているものを（部分的に）メモしておくなり。 

私の場合、スペクトルのグラフを描く機会が多い。
ライトカーブも連続測光した場合はクイックルック用にグラフを作りますが、
こちらはあまり見栄えなどはこだわっていない。

スペクトルのグラフ（気合いれた版）は T Pyx 論文で描いたような図が
最近のブームになっていて、ここ一年近く研究会の集録原稿などでも応用して使ってます。

以下は某新星の一晩分のスペクトルを作図したさいのコード。
フォントはデフォルトではなく、 Times 系を使うのがお気に入りです。


=====

　　setEPS()                     　　　　                 # EPS で出力
　　postscript("nova_hogehoge.eps",          # 出力ファイル名
                family="Times",                       # フォント
               height = 6, width = 7)              # グラフの大きさ            

　　par(tcl=0.4,lwd=2,ps=11)                     # 目盛の長さ、枠線、文字サイズなどの指定
　　plot(data, type="l", axes=F,                  # 外枠や軸はあとでプロットするのでaxes=F
　　　　　cex.lab=1.2,lab=c(7,7,8),
　　　xlab=expression(wavelength (ring(A))), 　# 横軸の波長Å
　　　ylab="log10(flux) + const.", col="red", lwd=1,xlim=c(4000,8500),ylim=c(0,2.5))

# gnuplot みたいに上下左右に目盛り線が欲しい場合は以下の通り。 

　　axis(1, lwd.ticks=2)                               # 軸を描く
　　axis(2, lwd.ticks=2,las=1)                      # las=1 は縦軸の数字を横書きにする
　　axis(3, lwd.ticks=2,lab=F)                      
　　axis(4, lwd.ticks=2,lab=F)
　　box()                                                   # 外枠を描く

# 以下、スペクトル線の元素を示す場合に使う。
# text (横軸位置, 縦軸位置, "出力文字", ...) を使う。
# srt=90 は文字を90°回転させている。 

　　text(6563,2.4,expression(paste("H",alpha)),srt=90,cex=1.2)     # Balmer
　　text(4861,2.0,expression(paste("H",beta)),srt=90,cex=1.2)
　　text(4340,2.0,expression(paste("H",gamma)),srt=90,cex=1.2)

　　text(6300,1.5,"[O I]",srt=90,cex=1.2)   # [O I]
　　text(7773,1.5,"O I",srt=90,cex=1.2)     # O I

# 以下、スペクトル１本のときはあまり使わないけど紹介。
# あるスペクトル線を点線などで強調したい場合は abline を使っている。
# 複数のスペクトルをオフセットして並べ、ラインの変化などを見るとき便利。

　abline(v=6563,lty=2,lwd=1)    # H-alpha
　abline(v=7773,lty=2,lwd=1)    # O I

# abline(v=強調したい波長, ...)、となっている。
# つまり指定した横軸の値に対して垂線が引かれる。 
# この点線は黒色で多様すると見栄えが微妙になることもある。
# 灰色など控えめの色を使うと幸せになれるかもしれない。 

　　dev.off()

======


とまぁ、 実際プロットした例図（スペクトル）をここに UP できれば良いのですが、
ブログでは割愛します (^^;
あと T Pyx 論文では補助目盛を使っていて、  
 

=====

　　xm <- seq(3800,9000,by=100)
　　　axis(1,xm,lab=F,tcl=0.4)
　　　axis(3,xm,lab=F,tcl=0.4)
　　ym <- seq(0,8,by=0.2)
　　　axis(2,ym,lab=F,tcl=0.4)
　　　axis(4,ym,lab=F,tcl=0.4)

=====

という感じで刻むことができる。この場合、主目盛の長さは tcl=0.6 くらいあると
良いかもしれない。このあたりは好みの問題でしょうか。

横軸を共有して二つのグラフを合体させたい場合の方法は、一度ブログで紹介しています。 
T Pyx 論文では plot() の前に以下のおまじないを唱えました。

=====
# グラフのレイアウトを行列で考える(一つの図の大きさの割合)

　　　mat <- matrix(c(1,1,1,1,2,2,2,2), # グラフのレイアウトを行列で考える(一つの図の大きさの割合)
　　　　　　　　　　　　　　　8,1,byrow=T)

　　　setEPS()  
　　　postscript("tpyx_spectra.eps",family="Times",
　　　　　　　　　　　　height=9,width =9)

　　　　 par(mfrow=c(2,1))    # 実際のグラフ数を 2x1 として縦に並べる
　　　　 layout(mat)
　　　　 par(mar=c(0,4,0,4))  # 各グラフをぴったり接続するおまじない(上下マージン0)
===== 

 

stellaR を入れてみた

StellaR: a Package to Manage Stellar Evolutio Tracks and Isochrones
http://arxiv.org/abs/1301.3695

VSOLJ ML で紹介されていたので、さっそくパッケージを入れてみました。
とりあえずここから .zip を落としてきて、ローカルからインストールしました。 

詳しい内容はさておき (おいっ) 、動作確認をしてみた。
コマンドは論文の通りにポチポチ入れていくと、気が付くとかっちょいい図が数種類できあがる。 

　showComposition() 

から順を追っていくと良い。
出来上がった図の例は以下の通り。
エラーなども無く、とりあえず問題なく動作したと思う。

R で星図と星座線をかく

R で Bright Star Catalog を読み込んで星をプロットし、
ついでに星座線まで出るようにしてみた。

ありがたいことにカタログと星座線のセットが以下に .csv 形式で置いてあるので使わせて頂くことに。
　星空横丁 (http://hoshizora.yokochou.com/)

星をプロットする場合、星の大きさをどうするか悩んだのですが
とりあえず kawashima さんの方法を使わせて頂きました。

あと星座線は RA が 例えば一方が23時で、もう一方が0時とかになると (例えばペガスス座の四辺形)
そのまま segments() で描くと酷いことになります。
そのようなケースが7件ほどあるので、とりあえず取り除いています。
# なので星座線つきの星図としてはちょっと未完成品orz

以下プロットできた星図です。


クリックしたら原寸画像で表示。

新星などが出たときに、出現場所を大きめの星図でまず確認するので
そのような場合に活用するかもです。

あとは星座境界線なんかも入れてみたいところです。
使いやすそうなデータセットがないか探し中。

以下、R のコード。途中星座線の座標を拾うのに for なんてものを使ってるので遅いです。
常用する場合は星座線の座標は別途保存しておいて
それを読み込むほうが実用的でしょう。

======

#------任意の領域指定
RA_lim  <- c(24,0)
DEC_lim <- c(-70,70)
#------星表の読み込みと星のプロット
data <- read.csv("star.csv",header=F)
line <- read.csv("star_line.csv",header=F)
　ra  <- data$V6 + (data$V7/60) + (data$V8/3600)
　dec <- data$V10 + (data$V11/60) + (data$V12/3600) 
　dec <- ifelse(data$V9==0, dec, dec*-1) # Dec の符号
　mag <- data$V14 + data$V15/100
　HR  <- data$V1
　star<- data.frame(HR,ra,dec,mag)

maxmin  <- summary(mag)
　r <- maxmin[[6]]+1-mag
symbols(ra,dec,circle=r,xlim=RA_lim,ylim=DEC_lim,
bg="black",fg="white",inches=0.06)
#------星座線の座標を拾う
t1<-NULL
t2<-NULL
l1 <- NULL
l2 <- NULL
 for(i in 1:744){
   l1 <- subset(star,star$HR==line[i,2])
   t1 <-rbind(t1,l1)
   l2 <- subset(star,star$HR==line[i,3])
   t2 <-rbind(t2,l2)
#------RA が折り返しそうなlineを拾う
   lr <- data.frame(t1,t2)
   lr1 <- subset(lr,lr$ra>22)
   lr1 <- subset(lr1,lr1$ra.1<3)
   lr2 <- subset(lr,lr$ra.1>22)
   lr2 <- subset(lr2,lr2$ra<3)
   lr0 <- rbind(lr1,lr2)
 }
#------星座線をかく (一部の星座線を取り除いている)
rown <- as.integer(row.names(lr0))
lr <- lr[-which(rownames(lr)%in% rown),]
segments(lr$ra,lr$dec,lr$ra.1,lr$dec.1)

赤道座標を銀河座標に変換2

前回の続きです。

VSOLJ のメーリングリストでご指摘と解説を頂きました (vsolj 11105)。
ありがとうございます。

角度変換に asin とか atan は使っちゃダメで、
というかそれ以前の問題で、私アホなことに定義域とか値域のことも
すっかり忘れていました・・・orz

で、早速前回のコードをいじくりました。
今度はうまくいったように見えます。
次回は適当な星表を読み込んでみて、銀河座標でプロットしてみようと思います。


=======

RA <- 72
DEC<- -10
#-------------
aN <- 282.25
I <- 62.6

 rad <- pi/180
 deg <- 180/pi

RA <- RA*rad
DEC<- DEC*rad
aN <- aN*rad
I  <- I*rad
#------------
x <- cos(DEC)*cos(RA-aN)
y <- sin(DEC)*sin(I) + cos(DEC)*sin(RA-aN)*cos(I)
z <- sin(DEC)*cos(I) - cos(DEC)*sin(RA-aN)*sin(I)

l <- atan2(y,x) *deg +33
l <- ifelse(l<0,360+l,l)
b <- atan2(z,sqrt(x*x+y*y)) *deg

cat("銀経=",sprintf("%3.2f",l),", 銀緯=",sprintf("%3.2f",b),"\n")