2014-12-03

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

昨日の記事の続き
球面調和関数の和で表す
- ルジャンドル多項式というのがある
  - $x=\pm 1$ の特異点を持つ多項式群である。ある微分方程式の解を構成する
- この多項式は $(n+1)P_{n+1}(x)=(2n+1)xP_n(x)-nP_{n-1}(x)$ , $P_0(x)=1,P_1(x)=x$ という漸化式であることが知られているので、これを使って、ルジャンドル多項式の係数を出す関数を作っておく

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)$
- 半径をの関数()として表せばよい
  - 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}$
  - ３つの構成要素からなっている
    - $\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)

表示の仕方
- 半径１の球面の各点の半径がいくつになるのかを決めているが、その視覚表示にはいくつかある
- １つは、半径の値を色であらわすもの
- もう一つは、半径の二乗によって、球面上の点の位置を拡大するもの
- この２つが、こちらにある

図

- もう一つは、実際の半径を使って点の座標を定める方法
  - これをすると、半径が負の場合が出る
  - 今回のこの記事を書いている理由は、球面調和関数を使って、球体をいびつにすることなので、この方法がよい

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))
# 半径の２乗で
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))