2012-01-27

美しくない

R 球球面三角形角座標

こちらで、多次元の球面三角形の面積をモンテカルロで算出している
角座標から普通の球の球面三角形の面積を出すとしたらどうなるかをやってみる
結論から言うと、すごく汚い式の積分になるようなので、「方針が間違っている」と思う
簡単のために、 $(1,0,0),(cos(s),0,sin(s)),(cos(t),sin(t),0)$ の３点が囲む「直角球面三角形」をの面積を考えることにする
$(1,0,0)$ から点U $(cos(u),sin(u),0)$ へと $u$ をパラメタとして動くことを考える
このとき点Uを通って、xy平面に垂直な大円が三角形を切るその弧の長さを考える
三角形の面積は弧のuに関する積分になる(はず)
この弧の長さは、大円が三角形の斜辺と交わる点の仰角に相当する値になるから、この交点の座標を角座標で表せれば( $u$ で仰角を表せれば)よさそうだ
交点の座標を $(cos(u)cos(v),sin(u)cos(v),sin(v))$ と表わすと、「仰角」は $v$
交点が、斜辺上にあることから、２点 $(cos(s),0,sin(s)),(cos(t),sin(t),0)$ を通る面上にあることになるが、この面の法線ベクトル $(x,y,z)$ は $x cos(s)+z sin(s)=0$ を満足し、 $x cos(t) + y sin(t)=0$ を満足している
これを利用してを条件としていることからの関係を導くと
- $tan(v)=tan(s)\times(cos(u)-sin(u)\frac{1}{tan(t)})$ となるらしい
したがって $\int_{0}^t v du$ が三角形の面積なのではないかと思うけれども、いかにも、美しくない
三角形が球面上で作る角からスタートすると、いかにもきれいになることからして、角座標から書き起こすのは得策ではないようだ

t<-seq(from=0,to=1,length=100)*2*pi

L1<-cbind(cos(t),sin(t),rep(0,length(t)))
L2<-cbind(cos(t),rep(0,length(t)),sin(t))
L3<-cbind(rep(0,length(t)),cos(t),sin(t))


CategoryVector <-
function (nc = 3) 
{
    df <- nc - 1
    d <- df + 1
    diagval <- 1:d
    diagval <- sqrt((df + 1)/df) * sqrt((df - diagval + 1)/(df - 
        diagval + 2))
    others <- -diagval/(df - (0:(d - 1)))
    m <- matrix(rep(others, df + 1), nrow = df + 1, byrow = TRUE)
    diag(m) <- diagval
    m[upper.tri(m)] <- 0
    as.matrix(m[, 1:df])
}
RandomSphere <-
function (df = 3, r = 1, n = 100) 
{
    rs <- matrix(rnorm(df * n), nrow = n)
    rs/sqrt(apply(rs^2, 1, sum)) * r
}

PerpFoot<-function(V,A,P=NULL){
	if(is.null(P))P<-rep(0,length(A))
	
	x<-solve(V%*%t(V))%*%V%*%(A-P)
	B<-P+t(V)%*%x
	return(list(X=B,k=x))
}

R.sp<-RandomSphere(3,1,100000)

tmpcol<-rep(0,length(R.sp[,1]))
s1<-pi/5
s2<-pi/3
V<-matrix(c(1,0,0,cos(s1),sin(s1),0,cos(s2),0,sin(s2)),3,3)
V<-t(V)
for(i in 1:length(R.sp[,1])){
	foot<-PerpFoot(V,R.sp[i,])
	posi<-foot$k>0
	if(prod(posi))tmpcol[i]<-2
	#tmpcol[i]<-sum(2^(0:2)*posi)
}
library(rgl)

Pts<-rbind(L1,L2,L3)

Pts2<-rbind(R.sp[which(tmpcol==2),],Pts)

col<-c(tmpcol[which(tmpcol==2)],rep(1,length(Pts[,1])))

s<-seq(from=0,to=s1,length=100)
p<-atan(tan(s2)*(cos(s)-sin(s)/tan(s1)))
A<-cbind(cos(s)*cos(p),sin(s)*cos(p),sin(p))

Pts3<-rbind(Pts2,A)
col<-c(col,rep(3,length(A[,1])))


plot3d(Pts3,col=col)