球面調和関数の値をPython3/matplotlibで３D表示する(2)

ネットで見つけた、球面調和関数１５式でグラフを作る。
コードは、こちら

```harmony2.py
import math
from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

def sgn(x):
  if (x >= 0.0):
    return 1.0
  else:
    return -1.0

def f(x,y,z,cs):
  r = math.sqrt(x*x + y*y + z*z)
  th = math.acos(z/r)
  ph = sgn(y)*math.acos(x/math.sqrt(x*x + y*y))
  if cs == 0:
    return (- x*y - y*z - z*x, 'f=-xy-yz-zx')
  if cs == 1:
    return (x*x + y*y - 2*z*z, 'f=x^2 + y^2 - 2z^2')
  if cs == 2:
    return (x*y*z, 'f=xyz')
  if cs == 3:
    return (z*(x*x - y*y), 'f=z(x^2-y^2)')
  if cs == 4:
    return (z*(5*z*z - 3*r*r), 'f=z(5z^2 - 3r^2)')
  if cs == 5:
    return ((-3*x*x*z - 3*y*y*z + 2*z*z*z)/(r*r*r), 'f=(-3xxz-3yyz+2z^3)/r^3')
  if cs == 6:
    return (35*z*z*z*z - 30*z*z + 3, 'f=35*z^4 - 30z^2 + 3')
  if cs == 7:
    return (x*z*(7*z*z - 3), '1-1-0.5-1-1 formation')
  if cs == 8:
    return ((x*x - y*y)*(7*z*z - 1), '4-4-4 formation')
  if cs == 9:
    return (x*z*(x*x - 3*y*y), '6-6 formation')
  if cs == 10:
    return (x*x*x*x - 6*x*x*y*y + y*y*y*y, '0-8-0 formation')
  if cs == 11:
    return (math.cos(2*ph)*(math.sin(th)**2)*(3*(math.cos(th)**3) - math.cos(th)), '4-4-4-4 formation')
  if cs == 12:
    return (math.cos(3*ph)*(math.sin(th)**3)*(9*(math.cos(th)**2) - 1), '6-6-6 formation')
  if cs == 13:
    return (math.cos(3*ph)*(math.sin(th)**3)*(11*(math.cos(th)**3) - 3*math.cos(th)), '6-6-6-6 formation')
  if cs == 14:
    return (math.cos(5*ph)*(math.sin(th)**5)*(math.cos(th)), '10-10 formation')

for cs in range(0,15):
  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,title) = f(x,y,z,cs)
      udata[c] = u
      c = c + 1
      s = s + u
  ax.set_aspect('equal')
  ax.scatter3D(xdata, ydata, zdata, c=udata, cmap=cm.coolwarm);
  plt.title(title)
  plt.pause(2)
plt.pause(200)
```

<table><tr>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/8aa1cbc9-4bd8-a712-8fbd-8fe6b10b2a44.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/d199a69f-2394-bb96-a5a0-62df3668a4a0.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/3eb0ce4d-2a39-fc9b-96e1-49698bfefc67.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/6e5b4ae0-f129-940c-9312-7f4f6f2b15d2.png"></td>
</tr>
<tr>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/594cbcd0-737a-9fa7-3092-8821a07d8eb4.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/92fb076a-1ded-638f-130a-de60831be696.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/d15ee88e-7cb7-a590-37bb-497a1fcf6b9f.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/37546009-eeaa-e4d5-ee86-6054c485089c.png"></td>
</tr>
<tr>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/681dd08e-958f-ed8b-717d-eb4c22b23362.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/a4c89a42-c70c-76b4-b497-bb8a38e4b4b9.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/8fa2bc61-206a-3c40-3083-0baedf792330.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/2f8a9cd2-fec9-3ee8-dfac-ff3308bde27a.png"></td>
</tr>
<tr>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/94702944-9010-05d7-d67e-16f756a2509b.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/8951ec7a-ddde-605f-68d5-a8adcbabf02a.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/c76ba21b-ccfc-f2f4-604e-35202ab5185d.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/204041/2f8a9cd2-fec9-3ee8-dfac-ff3308bde27a.png"></td>
</tr>
</table>

それぞれの図において、球面上の点を均等にサンプルして色表示しています。赤がプラス。青がマイナス。
赤と青の濃さを合計するとゼロになります。
これから、電子軌道が導き出せますが、私には説明できません。