More than 3 years have passed since last update.

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

Last updated at 2021-10-13Posted at 2021-10-12

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

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)