球面調和関数をR でやる - 驚異のアニヲタ社会復帰の予備

MikuHatsune2016-09-18

球面調和関数をrPython でやっていたが、R だけで完結するほうが、rPython に投げてやりとりするより圧倒的に高速かつ安定なので、頑張ってR でやった。
ゴールとしては、
$f(\theta, \phi)=\sum_{l}^{\infty}\sum_{m=-l}^lc_l^m Y_l^m$
となる係数 $c_l^m$ を求めることだが、これはルジャンドル陪関数が直交するという性質を利用して
[tex:c_l^m=] (内積)
を求めることでなんとかなる。
degree $l$ は理想的には無限大まで行けば、 $f(\theta, \phi)$ は完璧に級数展開されるが、現実的には30か40くらいまでいけば、たいていの形は再現されるっぽいようなことが論文によく書いてあるから、 $L_{max}=25$ くらいまでがんばった。
　
orthopolynom というパッケージで、ルジャンドル関数までは用意してくれるが、 $m$ 階の微分は自分で用意しないといけない。
実際に3次元物体を扱うには実数部分が必要だから、sph1 関数で最終的に必要な実部もしくは虚部だけ取ってくる。

library(orthopolynom)
l <- 5
m <- 2
theta <- pi/7
phi <- pi/8
p_l_m <- function(l, m){ 
  legp <- legendre.polynomials(l)[l+ 1]
  for(mm in 1:m){
    legp <- polynomial.derivatives(legp)
  }
  polynomial.coefficients(legp)[[1]]
}
p_l_m <- function(x, l, m){
  legp <- legendre.polynomials(l)[l+1]
  polynomial.values(eval(parse(text=paste0(paste(rep("polynomial.derivatives(", abs(m)), collapse=""), "legp", paste(rep(")", abs(m)), collapse=""))))[1], x)[[1]]
}

sph <- function(theta, phi, l, m){
  x <- cos(theta)
  y <- cos(m*phi) + 1i*sin(m*phi)
  #Plm <- (-1)^m * (1-x^2)^(m/2)*p_l_m(x, l, m)
  Plm <- sin(theta)^abs(m)*p_l_m(x, l, m)
  #res <- (-1)^((m+abs(m))/2)
  res <- (-1)^((m+abs(m))/2)*Plm*sqrt((2*l+1)/4/pi*factorial(l-abs(m))/factorial(l+abs(m)))*y
  return(res)
}

sph1 <- function(theta, phi, l, m){
  s <- sph(theta, phi, l, m)
  ifelse(rep(m, length(s)) >= 0, Re(s), Im(s))
}

$Y_0^0=\frac{1}{\sqrt{4\pi}}$ だから確かめると、合ってる。

sph1(0.5, 0.5, 0, 0)

[1] 0.2820948

1/sqrt(4*pi)

[1] 0.2820948

　
単位球面上に均一な点を発生させて、 $f(\theta,\phi)=1$ ができるかどうか検算してみる。gss というオブジェクトが単位球で、SphericalCubature パッケージのrect2polar 関数で極座標に変換できる。

library(rgl)
library(geometry)
library(SphericalCubature)
# フィボナッチ数を返す関数
FN <- function(n){
	# 負の数のフィボナッチ数は、正のそれに符号をつけたもの
	if(n<0){
		return(both.fn.3(-n)*(-1)^(n+1))
	}
	# 0,1,2はフィボナッチ数発生の種なので指定する
	if(n==0){
		return(0)
	}else if (n==1 | n==2){
		return(1)
	}else{# 3以上は、計算する
		tmp <- rep(0,n)
		tmp[1] <- tmp[2] <- 1
		for(i in 3:n){
			tmp[i] <- tmp[i-1]+tmp[i-2]
		}
		return(tmp[length(tmp)])
	}
}

# 球面らせん法
# f1は黄金比
fib.lattice.S2 <-function(N,f1=(sqrt(5)+1)/2,k=1){
	f2 <- f1-1 # 黄金比はx^2-x-1=0の１つの根。もう１つの根を取り出す
	# 点の数2N+1
	P <- 2*N+1
	# -N から Nまでの整数列
	i <- seq(from=-N, to=N,by=k)
	# zの座標は-1 から 1の均等割り
	theta <- asin(2*i/(2*N+1))
	z <- sin(theta)
	# そのようなz座標に対応する、(x,y)座標を取り出すにあたり
	# 黄金比から作ったf2に応じてxy平面での角座標を作る
	phi <- 2*pi*i*f2
	# その角を使って、z座標に対応した半径の座標とする
	x <- cos(theta)*cos(phi)
	y <- cos(theta)*sin(phi)
	return(cbind(x,y,z))
}
# FN(n)個の点をつくるために、引数をFN(n)/2とする
gss <- fib.lattice.S2(FN(18)/2)
rad <- rect2polar
rad <- t(rect2polar(t(gss))$phi)

# l はLmax まで、m は -l, ..., 0, ..., l までひたすら計算する。
Lmax <- 25
ylm <- function(theta, phi, Lmax) mapply(function(l) mapply(function(m) sph1(theta, phi, l, m), -l:l), 0:Lmax)

# theta, phi に対応して ylm を計算する
# 計算結果は、m+1 番目のリストに、(theta, phi) の組の数だけの行と
# -l, ..., 0, ..., l の列の行列
w <- ylm(rad[,1], rad[,2], Lmax)

# 係数は内積で得られる
coef <- function(vec, Y) mapply(function(z) c(vec %*% z), Y)
# 単位球なので f(theta, phi)=1 である
a <- coef(rep(1, nrow(rad)), w)

# あるm について、-l ~ l まですべて足す
# その後、m についても足す
r <- 4*pi*rowSums(sapply(mapply(function(z1, z2) sweep(z1, 2, z2, "*"), w, a), rowSums)

# ひとつの頂点が担当している球面上の領域は、均一なので頂点数で割ればいいはず
r <- r/nrow(gss)

hist(r)