3次元図形についてのちょっとしたメモ

my.Legendre <- function(k){
	Pm <- list()
	Pm[[1]] <- c(1)
	Pm[[2]] <- c(0,1)
	if(k <= 1){
		return(Pm)
	}
	for(i in 3:(k+1)){
		z <- i-1
		x1 <- c(0,Pm[[i-1]])
		x2 <- Pm[[i-2]]
		L.x1 <- length(x1)
		L.x2 <- length(x2)
		L.max <- max(L.x1,L.x2)
		x1 <- c(x1,rep(0,L.max-L.x1))
		x2 <- c(x2,rep(0,L.max-L.x2))
		Pm[[i]] <- 1/(z)*((2*z-1)*x1 - (z-1)*x2)
	}
	Pm
}
    • さらに、この多項式係数ベクトルに値をあたえて、関数値を出す関数も作っておく
my.poly.calc <- function(Pm,x){
	n <- length(Pm)
	ret <- rep(0,length(x))
	for(i in 1:n){
		ret <- ret + Pm[i] * x^(i-1)
	}
	ret
}
k <- 20
Pms <- my.Legendre(k)
ys <- lapply(Pms,my.poly.calc,x)
plot(x,ys[[1]],type="l",ylim=c(-1,1))
for(i in 2:length(ys)){
	points(x,ys[[i]],type="l",col=i)
}

  • このルジャンドル多項式を使って、球面調和関数を記述する
    • 球面座標は極座標表示で表す。経度を0 \le \theta < 2\piとし、緯度は、極を0として0 \le \phi \le \piで表す
    • x = r \cos(\theta)\sin(\phi)
    • y=r\cos(\theta)\sin(\phi)
    • z=r\cos(\phi)
    • 半径を\theta,\phiの関数(r(\theta,\phi))として表せばよい
      • Wolfram MathWorldを基本に使用(\theta,\phiがこの記事と逆であることに注意)
      • ごちゃごちゃと書いてあるけれど、球面調和関数は以下のように表せる
      • Y_m^n (\theta,\phi)=\sqrt{\frac{2m+1}{4\pi}\frac{(m-n)!}{(m+n)!}}P_m^n(\cos(\phi))e^{i n \theta}
      • 3つの構成要素からなっている
        • \sqrt{\frac{2m+1}{4\pi}\frac{(m-n)!}{(m+n)!}}は、関数の値がそれなりになるような「調整係数」
        • P_m^n(\cos(\phi))。m,nは整数でP_m^n(\cos(\phi))は、associated Legendre 多項式と呼ばれるもので、それはm,nで式が決まり、\cos(\phi)はその引数
        • e^{i n \theta}=\cos(n\theta) + i \sin(n\theta)であるが、これは、経度\thetaに関して、n=1ならただの円だし、大きくなれば周回が増えるように定めたもの。要するに経度について、半径を周期的に変化させる成分。実際に半径を考えるときには、実数部分を考えればよい
      • 問題はP_m^n(\cos(\phi))。これは緯度についての関数。実は、ルジャンドルm次多項式のn階微分と関係している。ルジャンドル多項式-1 \le x \le 1の範囲でうまく閉じた関数であるから、そのxに、-1 〜 1をとる値\cos(\phi)を(これはz座標に対応する)。このように、この成分は、緯度に関して周期性を持たせるための成分。
      • 実際P_m^n(x) = (-1)^n(1-x^2)^{n/2} \frac{d^n}{dx^n}P_m(x)、ただしP_m(x)はm次ルジャンドル多項式
      • これは、m次ルジャンドル多項式のn階微分と、(-1)^n(1-x^2)^{n/2}との積になっている
      • さらにx=\cos(\phi)を考慮すると(-1)^n(1-\cos(\phi)^2)^{n/2}=(-1)^n(\sin(\phi))^{n/2}となる。ここで、球面調和関数が「極座標」と親和性が高いことがわかる
      • ルジャンドル多項式のn階微分をするので、多項式係数を微分して変化する関数を作っておく
my.differential <- function(x,k){
	n <- length(x)
	if(k==0){
		ret <- x
	}else if(n < (k+1)){
		ret <- c(0)
	}else{
		tmp <- x[(k+1):n]
		a <- 1:(n-k)
		b <- (k):(n-1)
		ret <- factorial(b)/factorial(a) * tmp
	}
	ret
}
      • 球面多項式を足し合わせることを考慮した関数を作る。ここでは、0,1,...,m次の多項式の、0,1,...k階微分の重みを(m+1)x(m+1)行列に納めることとする
      • また、経度に関して位相変化をさせたいので、それについても(m+1)x(m+1)行列に納めることにする。
my.spherical.harm <- function(A,B=matrix(0,ncol=ncol(A),nrow=nrow(A)),Pm=my.Legendre(length(A[,1])),theta=seq(from=0,to=1,length=100)*2*pi,phi=seq(from=0,to=1,length=100)*pi,normalization=TRUE){
# 角度の座標
	tp <- expand.grid(theta,phi)
	z <- rep(0,length(tp[,1]))
	M <- ncol(A)
# 重みづけ行列でのループ
	for(i in 1:M){
		for(j in 1:i){
			tmp1 <- (A[i,j] * cos((j-1)*tp[,1] + B[i,j]) ) * (-1)^(j-1)*(1-cos(tp[,2])^2)^((j-1)/2)
# ルジャンドル多項式のj-1階微分の係数ベクトル
			tmp2 <- my.differential(Pm[[i]],j-1)
			tmp3 <- cos(tp[,2])
			tmp4 <- rep(0,length(z))
			for(k in 1:length(tmp2)){
				tmp4 <- tmp4 + tmp3^(k-1) * tmp2[k]
			}
			if(normalization){
				tmp5 <- sqrt((2*(i-1)+1)/(4*pi)*factorial(i-j)/factorial(i+j-2))
			}else{
				tmp5 <- 1
			}
			
			z <- z + tmp1 * tmp4 * tmp5
		}
	}
	X <- z * cos(tp[,1]) * sin(tp[,2])
	Y <- z * sin(tp[,1]) * sin(tp[,2])
	Z <- z * cos(tp[,2])
	X. <- cos(tp[,1]) * sin(tp[,2])
	Y. <- sin(tp[,1]) * sin(tp[,2])
	Z. <- cos(tp[,2])
	return(list(r=z,tp=tp,X=cbind(X,Y,Z),X.=cbind(X.,Y.,Z.)))
}
n <- 6
A <- matrix(0,n,n)
ret <- matrix(list(),n,n)
for(i in 1:n){
	for(j in 1:i){
		A. <- A
		A.[i,j] <- 1
		ret[i,j] <- list(my.spherical.harm(A.))

	}
}

tmp <- ret[4,1][[1]]
out <- rbind(tmp$X,rep(min(tmp$X),3),rep(max(tmp$X),3))
plot3d(out,type="l",col=rainbow(128))

col1 <- col2 <- rep(0,length(tmp$r))
col1[which(tmp$r>=0)] <- tmp$r[which(tmp$r>=0)]
col2[which(tmp$r<0)] <- -tmp$r[which(tmp$r<0)]
col1 <- col1/max(col1)
col2 <- col2/max(col2)
plot3d(tmp$X.,col=rgb(1-col1,1-col2,0),size=10)
  • 表示の仕方
    • 半径1の球面の各点の半径がいくつになるのかを決めているが、その視覚表示にはいくつかある
    • 1つは、半径の値を色であらわすもの
    • もう一つは、半径の二乗によって、球面上の点の位置を拡大するもの
    • この2つが、こちらにある

    • もう一つは、実際の半径を使って点の座標を定める方法
      • これをすると、半径が負の場合が出る
      • 今回のこの記事を書いている理由は、球面調和関数を使って、球体をいびつにすることなので、この方法がよい
n <- 6
A <- matrix(0,n,n)
ret <- matrix(list(),n,n)
for(i in 1:n){
	for(j in 1:i){
		A. <- A
		A.[i,j] <- 1
		ret[i,j] <- list(my.spherical.harm(A.))

	}
}
# 半径そのまま
tmp <- ret[4,1][[1]]
out <- rbind(tmp$X,rep(min(tmp$X),3),rep(max(tmp$X),3))
plot3d(out,type="l",col=rainbow(128))
# 半径の2乗で
out.. <- rbind(tmp$X*tmp$r,rep(min(tmp$X*tmp$r),3),rep(max(tmp$X*tmp$r),3))
plot3d(out..,type="l",col=rainbow(128))
# 半径を色で
col1 <- col2 <- rep(0,length(tmp$r))
col1[which(tmp$r>=0)] <- tmp$r[which(tmp$r>=0)]
col2[which(tmp$r<0)] <- -tmp$r[which(tmp$r<0)]
col1 <- col1/max(col1)
col2 <- col2/max(col2)
plot3d(tmp$X.,col=rgb(1-col1,1-col2,0),size=10)
      • ただし、基本となる正球の上に、関数を足す必要があるので、以下のようにする
n <- 6
A <- matrix(0,n,n)
ret <- matrix(list(),n,n)
for(i in 1:n){
	for(j in 1:i){
		A. <- A
		A.[1,1] <- 5 # 正円の分
		A.[i,j] <- 1
		ret[i,j] <- list(my.spherical.harm(A.))

	}
}

tmp <- ret[4,1][[1]]
out <- rbind(tmp$X,rep(min(tmp$X),3),rep(max(tmp$X),3))
plot3d(out,type="l",col=rainbow(128))