2010-09-27

分割方法を点で表わす

三角関数球球面三角形 R

２次元平面上の単位円上に点をとる。第1象限の点の角座標は、第１象限の弧の２分割を表している
そのアナロジーで、３次元空間にある単位球の表面の第１象限を考える。そこに１点をとり、３つのデカルト正規直交座標単位ベクトルの先端への大円にて、第１象限の部分球面を分割することとする
関連記事はこちらやこちらで書いた通り、デカルト座標 $(x_i),i=1,2,3$ と、大円の弧の長さ $\theta_i,i=1,2,3$ とには $\cos(\theta_i)=x_i$ なる関係がある
２次元単位円の場合には、『面積』と『弧の長さ』が同じなので、『円周という面積(長さだが)の分割』は、弧の長さそのものであり、それは、角度そのものである
３次元のとき、球表面の分割と $x_i,\theta_i$ とは、どういう関係になるのか、ということを問題にしている
$x_i$ と $\theta_i$ との関係は、 $f(x_i)=\theta_i$ という１対１の関係が３個あった
球表面積を３分割(球面三角形の対辺を底辺とする球面三角形を、ある点に対応する領域と考えた、分割)を考える
この部分面積と $(x_i)$ との関係は、3変数関数を３個用いて定められる
球面上の点に接する面と軸との交点がであることや、その交点ととを結んだベクトルが、と直交することなどを使い、球面三角形の面積が、その内角を用いるとで表わされることを使うと、に対応する部分の面積は
- $A_1 = \tan^{-1}(x_1/x_3)+\tan^{-1}(x_1/x_2) + \cos^{-1}(\frac{-1}{\sqrt{(\frac{1}{x_2^2}-1)(\frac{1}{x_3^2}-1)}})-\pi$
- $A_2 = \tan^{-1}(x_2/x_1)+\tan^{-1}(x_2/x_3) + \cos^{-1}(\frac{-1}{\sqrt{(\frac{1}{x_3^2}-1)(\frac{1}{x_1^2}-1)}}-\pi)$
- $A_3 = \tan^{-1}(x_3/x_2)+\tan^{-1}(x_3/x_1) + \cos^{-1}(\frac{-1}{\sqrt{(\frac{1}{x_1^2}-1)(\frac{1}{x_2^2}-1)}})-\pi$
ちょっとまだ複雑で、高次元に一般化しやすいのかどうか分かりにくいけれども・・・
整理していない覚書のＲのコードは以下の通り

RandomSphere<-function (df = 3, r = 1, n = 100) 
{
    rs <- matrix(rnorm(df * n), nrow = n)
    rs/sqrt(apply(rs^2, 1, sum)) * r
}


k<-3 # dimension 
cs<-abs(RandomSphere(df=k,n=1)) # random point coordinate in k-dimensional space
sum(cs^2) # it is on the surface of unit sphere in k-dimensional space

thetas<-acos(cs) # distance on the sphere from the point to regular base unit vectors


x<-cs[1]
y<-cs[2]
z<-cs[3]

A1<-acos(y/sqrt(y^2+z^2))
A2<-acos(z/sqrt(y^2+z^2))
B1<-acos(z/sqrt(x^2+z^2))
B2<-acos(x/sqrt(x^2+z^2))
C1<-acos(x/sqrt(x^2+y^2))
C2<-acos(y/sqrt(x^2+y^2))

A1<-atan(z/y)
A2<-atan(y/z)
B1<-atan(x/z)
B2<-atan(z/x)
C1<-atan(y/x)
C2<-atan(x/y)

P<-acos(-1/sqrt((1/x^2-1)*(1/y^2-1)))
Q<-acos(-1/sqrt((1/y^2-1)*(1/z^2-1)))
R<-acos(-1/sqrt((1/z^2-1)*(1/x^2-1)))

A1+B2+P-pi
B1+C2+Q-pi
C1+A2+R-pi



A1
A2
B1
B2
C1
C2
P
Q
R
A1+A2
B1+B2
C1+C2
P+Q+R



##################

#３D第一象限にある3点が作る立体角を求める
KyuumenSankaku<-function(v1,v2,v3){
 t12<-acos(sum(v1*v2))
 t23<-acos(sum(v2*v3))
 t31<-acos(sum(v3*v1))

 s1<-acos((cos(t23)-cos(t31)*cos(t12))/(sin(t12)*sin(t31)))
 s2<-acos((cos(t31)-cos(t12)*cos(t23))/(sin(t12)*sin(t23)))
 s3<-acos((cos(t12)-cos(t31)*cos(t23))/(sin(t31)*sin(t23)))

 area<-s1+s2+s3-pi

# 返り値のarea2は１のときに立体角がpi/2
 list(vs=matrix(c(v1,v2,v3),ncol=3,byrow=TRUE),t=c(t12,t23,t31),s=c(s1,s2,s3),area=area,area2=area/(pi/2))

}
#3D第一象限にある点が(1,0,0),(0,1,0),(0,0,1)と作る立体角(KyuumenSankaku()の返り値のarea2を返す
Rittaikaku3<-function(x){
 v1<-c(1,0,0)
 v2<-c(0,1,0)
 v3<-c(0,0,1)
 a12<-KyuumenSankaku(x,v1,v2)$area2
 a23<-KyuumenSankaku(x,v2,v3)$area2
 a31<-KyuumenSankaku(x,v3,v1)$area2
 c(a23,a31,a12)
}

# df 次元の半径rの球面上の点をランダムに均等にn個発生させる関数
RandomSphere<-
function (df = 3, r = 1, n = 100) 
{
    rs <- matrix(rnorm(df * n), nrow = n) # 多次元正規乱数を用いる
    rs/sqrt(apply(rs^2, 1, sum)) * r # 半径rにする
}
# 第一象限ランダム点の発生
Npoints<-10000
xs<-abs(RandomSphere(df=3,r=1,n=Npoints)) # ３成分の絶対値を取って第一象限化
thetas<-xs # ３立体角を格納する
for(i in 1:Npoints){
 thetas[i,]<-Rittaikaku3(xs[i,]) # 発生させた点について逐一、立体角を計算し格納
}

# xs と thetasとのコプロット
plot(as.data.frame(cbind(xs,thetas)))


par(mfcol=c(1,2)) # 描図領域を１行2列に分割
# 第一図
# thetaとxとの関係は、三角関数に似た形
# ただし、三角関数によって、下限が作られている
plot(thetas[,1],xs[,1]) 
par(new=TRUE) # 重ねて描図
plot(thetas[,1],sin(thetas[,1]*pi/2),col="red") #三角関数

# 第二図
# theta 3成分から、ディリクレ正三角形を作り、そこに、x第一成分の値でグラデーションをしてプロット
# 緩い曲線の等高線が見える
plot(thetas[,1],sqrt(3)/2*(thetas[,2]-thetas[,3]),col=gray(xs[,1]))