太陽高度を計算する

地球は太陽の周りを回る。
地軸が特定の角度だけ傾いた状態を維持したまま、地球は太陽の周りを回る。
その地軸を回転軸として地球は回る。
一周回ると一日。
地上の点の天頂方向は、地球の中心を原点として、地軸ベクトルを回転軸として回転する。

太陽の仰角(夜間は地平線の下に隠れているが、負の仰角を持つ)は、この天頂方向ベクトルと、太陽-地球ベクトルとが成す角を90度から引いた角度になる。

これをRで計算してみる。

# 緯度 degree
ido <- 35
ido.pi <- 35 / 180 * pi
# date of year
date <- 1 # 元日
# date of year, fraction
date.y <- date / (365 + 0.25)
date.y.pi <- date.y * 2 * pi
# 地軸の傾き(太陽に向かって)
psi <- 23.4
psi.pi <- 23.4 / 180 * pi
# 時刻
t <- 0 # 0 ... 24 hours


# 地球の回転軸単位ベクトルは、季節・時刻によらず一定
v.earth <- c(-cos(psi.pi),0,sin(psi.pi))
# 太陽光が地球に向かう(平行)光線単位ベクトルは季節によって変動、時刻によらず一定
v.solar <- c(cos(date.y.pi),sin(date.y.pi),0)

# 着目地点の天の頂上方向単位ベクトルは季節によらず、時刻によって変動
# そのベクトルは、あるベクトルを、地球の回転軸に関して時刻に応じた角度だけ回転させたもの
# この計算は四元数を用いると便利
# v.pointをある地点の、計算しやすい、特定の方向ベクトルとする
v.point <- c(-cos(ido.pi - psi.pi),0,sin(ido.pi - psi.pi))
# 以下がその計算関数
library(onion)

my.tencho.time <- function(v.earth,v.point,time){
	time.pi <- time/24 * 2 * pi
	v.earth.q <- Hi * v.earth[1] + Hj * v.earth[2] + Hk * v.earth[3]
	rotation.q <- cos(time.pi/2) + sin(time.pi/2) * v.earth.q
	rotated <- rotate(matrix(v.point,ncol=3),rotation.q)
	return(rotated)
}

# 太陽の高度は、天の頂上方向ベクトルと太陽光線ベクトルとの成す角で決まる
my.angle.solar <- function(tencho.v,solar.v){
	ip <- sum(tencho.v * solar.v) / (sqrt(sum(tencho.v^2)) * sqrt(sum(solar.v^2)))
	ang <- acos(ip) / pi * 180
	return(ang)
}

# ある日のある時刻の太陽高度を計算する関数
my.sun.angle <- function(ido,date,time,psi=23.4){
	# 地軸の傾き(太陽に向かって)
	#psi <- 23.4
	psi.pi <- 23.4 / 180 * pi
	# ベクトルの成分を求めるときには、地軸の傾きそのものではなく
	# pi/2 - psi.pi を用いる
	v.earth <- c(-cos(pi/2-psi.pi),0,sin(pi/2-psi.pi))
	# date of year, fraction
	date.y <- date / (365 + 0.25)
	date.y.pi <- date.y * 2 * pi
	
	ido.pi <- ido / 180 * pi
	# v.pointをある地点の、計算しやすい、特定の方向ベクトルとする
	# その方向は、北半球では夏至のときの地球の位置と太陽の位置を
	# 結ぶ線の方向に「地点」があるときの地点方向ベクトルを計算しやすい方向とする
	
	dir <- psi.pi
	v.point <- c(-cos(ido.pi - psi.pi),0,sin(ido.pi - psi.pi))

	tencho.v <- my.tencho.time(v.earth, v.point, time)
	
	# 太陽光が地球に向かう(平行)光線単位ベクトルは季節によって変動、時刻によらず一定
	solar.v <- c(cos(date.y.pi),sin(date.y.pi),0)

	ang <- 90 - my.angle.solar(tencho.v,solar.v)
	
	return(ang)
}

ido <- 35
# おおよそ、夏至、秋分、冬至、春分
dates <- (0:3) / 4 * 365 + 1
# １日。ただし、季節によって、このtimesの値のどこが「南中か」にはずれがあるので
# 太陽高度が一番高くなっているtimes 値を、南中時刻とみなすべき
times <- seq(from=0,to=24,length=200)
dates.times <- expand.grid(dates,times)

angs <- rep(0,length(dates.times[,1]))

for(i in 1:length(dates.times[,1])){
	date <- dates.times[i,1]
	time <- dates.times[i,2]
	angs[i] <- my.sun.angle(ido,date,time)
}

#plot(times,angs,type="l")

matplot(times,matrix(angs,ncol=length(dates),byrow=TRUE),type="l",xlab="時間",ylab="太陽高度")
abline(h=0) # 水平線

###############
# ある日の24時間の太陽高度を計算する関数

my.sun.angle.24 <- function(ido,date,n.times=24*10,psi=23.4){
	angs <- rep(0,n.times)
	times <- seq(from=0,to=24,length=n.times)
	for(i in 1:n.times){
		angs[i] <- my.sun.angle(ido,date,times[i],psi=psi)
	}
	idmin <- which.min(angs)
	ids <- c(idmin:n.times,1:(idmin-1))
	return(angs[ids])
}


ido <- 35
# おおよそ、夏至、秋分、冬至、春分
dates <- (0:3) / 4 * 365 + 1
# １日。ただし、季節によって、このtimesの値のどこが「南中か」にはずれがあるので
# 太陽高度が一番高くなっているtimes 値を、南中時刻とみなすべき

n.times <- 240

angs <- matrix(0,ncol=length(dates),nrow=n.times)


for(i in 1:length(dates)){
	date <- dates[i]
	angs[,i] <- my.sun.angle.24(ido,date,n.times)
}

#plot(times,angs,type="l")
times <- seq(from=0,to=24,length=n.times)
matplot(times,angs,type="l",xlab="時間",ylab="太陽高度")
abline(h=0) # 水平線
abline(v=c(0,12,24))