3D球面調和関数を、Python3/matplotlibの３Dスキャタ表示してみる

調和関数が、平均値の性質を満たしているか、Python3/matplotlib 3Dスキャタ表示をして試してみる。
平均値の性質を簡単に説明すると、r=1球面上で調和関数の値を積分すると0になる。という事である。
これだと、f(x,y,z)=2　のような調和関数には、当てはまらないので、f(0,0,0)=0の調和関数に限定して考える。

曲座標系の関数はわかりにくので、デカルト座標での関数を考える。

この論文に、デカルト座標の調和関数がリストされている。
http://lucille.sourceforge.net/blog/images/shpoly.pdf

$\sqrt{\frac{3}{4\pi}}$ のような係数は、今回は無視して１とする。
赤の＊は、重複しているので省略します。
$y_{0,0}$は、白一色になるので、省略

$y_{3,-2}=xyz$

GIFアニメーションを生成するコードは、こちら

harmony_animation.py

import math
from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib.animation as animation

def y3_2(x,y,z):
  return x*y*z

fig = plt.figure()
ax = plt.axes(projection='3d')

s = 0.0
c = 0
N = 50
zdata = np.linspace(0,0,2*N*N)
xdata = np.linspace(0,0,2*N*N)
ydata = np.linspace(0,0,2*N*N)
udata = np.linspace(0,0,2*N*N)

for j in range(0,N):
  ps = math.pi*(j+0.5)/N
  w = int(N*math.sin(ps))
  for i in range(-1*w,w):
    th = math.pi*(i+0.5)/w
    z = math.cos(ps)
    x = math.sin(ps)*math.cos(th)
    y = math.sin(ps)*math.sin(th)
    xdata[c] = x
    ydata[c] = y
    zdata[c] = z
    u = y3_2(x,y,z) #replace here!! y53,y63,y65
    udata[c] = u
    c = c+1
    s = s + u

ax.set_aspect('equal')
ax.scatter3D(xdata, ydata, zdata, c=udata, cmap=cm.coolwarm);
#plt.show()
def rotate(angle):
    ax.view_init(azim=angle)

print("Making animation")
rot_animation = animation.FuncAnimation(fig, rotate, frames=np.arange(0, 362, 2), interval=100)
rot_animation.save('rotation.gif', dpi=80, writer='imagemagick')