太陽の軌道

天球、天の赤道・黄道、地上からの天体の運動軌道に関して、三次元単位ベクトルを軸としたベクトルの回転としてとらえ、その計算を四元数を用いて行う話。

# 黄道と日の出・日の入り・南中
# 黄道は、北緯23.4度地点での天の赤道が
# その地点で見上げた空に作る軌道と同じ軌道を
# 天球上に描いたものなので
# 上記の関数 my.seizaban() の引数 x に23.4 / 180 * pi を
# 与えた場合と似通った処理により
# 太陽の天球上の座標を求めることができる
library(onion)
# day は太陽暦の1月1日からの日数
my.sun <- function(day){
	
	# 冬至を基準として１年で360度を回るとしたときの角度
	# 冬至を12/22としておこう
	angle.geshi <- (day + 9)/365 * 2 * pi
	# 北緯 23.4度での天の北極方向単位ベクトルの四元数表現
	theta <- 23.4/180 * pi
	N <- -cos(theta) * Hi + sin(theta) * Hk
	# 天の赤道の南中方向単位ベクトルの四元数表現
	P <- -cos(theta+pi/2) * Hi + sin(theta+pi/2) * Hk
	
	# N 周りにPを角度angle.geshiして得られる単位方向ベクトルは
	Nrot <- cos(angle.geshi/2) + sin(angle.geshi/2) * N
	# を用いて
	
	res <- Nrot * P * Conj(Nrot) 
	return(res)
}

# 一年のある日 day の赤緯は、my.sun()が返す単位ベクトルと
# 単位ベクトル(-1,0,0)　(天球における天の北極に相当） とのなす角

my.sun.sekii <- function(day){
	sun.v <- my.sun(day)
	# 天の北極単位ベクトルと太陽方向単位ベクトルとの内積
	ip <- (-1) * i(sun.v)
	# 天の北極を０度とした角度
	ang <- acos(ip)
	# 赤緯
	sekii <- -ang + pi/2
	return(sekii)
}

# ある赤緯の天体を、ある北緯地点で観測すると、どのような軌道で見えるかを返す関数

my.star.location <- function(hokui,sekii,jikoku=seq(from=0,to=24,length=24*60+1)){
	# 天の北極方向ベクトル
	N <- -cos(hokui) * Hi + sin(hokui) * Hk
	theta <- pi/2 - sekii
	theta2 <- hokui - theta
	# 南中の１２時間前に相当する方向ベクトル
	P <- -cos(theta2) * Hi + sin(theta2) * Hk
	
	# jikoku は hour単位
	# これを角度に変える
	jikoku.angle <- jikoku / 24 * 2 * pi
	
	Rot.q <- cos(jikoku.angle/2) + sin(jikoku.angle/2) * N
	
	loc.q <- Rot.q * P * Conj(Rot.q)
	
	res <- list(jikoku=jikoku,loc = as.matrix(loc.q))
	return(res)
}

# 緯度（北緯）と、一年の日と時刻から、太陽の軌道を得る

my.sun.loc <- function(hokui,day,jikoku=seq(from=0,to=24,length=24*60+1)){
	sekii <- my.sun.sekii(day)
	out <- my.star.location(hokui,sekii,jikoku=jikoku)
	loc <- t(out$loc[2:4,])
	res <- list(jikoku=jikoku,loc=loc)
	return(res)
}
# 地平を描き、天体の軌道を重ねる
library(rgl)

x <- seq(from=0,to=2*pi,length=100)
chihei <- cbind(cos(x),sin(x),rep(0,length(x)))
plot3d(chihei)

# 北緯３５度
hokui <- 35 / 180 * pi
# 春分・秋分の日の太陽の軌道
# 太陽は天の赤道上にある＝赤緯０
sekii <- 0 / 180 * pi
out <- my.star.location(hokui,sekii)
plot3d(t(out$loc[2:4,]),add=TRUE,col=2)

# 冬至は、12/22頃、それは、元日を基準にして -9 日
# 冬至のとき、太陽は天の赤道からもっとも南寄りにある

out1 <- my.sun.loc(hokui,-9)
plot3d(out1$loc,add=TRUE,col=3)

# 夏至は冬至に半年 365/2 日分　を加える
out2 <- my.sun.loc(hokui,-9+365/2)
plot3d(out2$loc,add=TRUE,col=4)