赤道座標を銀河座標に変換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")