原子核四极形变的 python 可视化

把原子核当作液滴，讨论较小的形变，经过一系列处理，可以用公式处理：

\[R = R_0 \left\{ 1 + \beta \sqrt{\frac{5}{16\pi}} (\cos\gamma (3\cos^2\theta -1) + \sqrt{3}\sin\gamma\sin^2\theta\cos2\phi) \right\}. \]

import numpy as np
import matplotlib.pyplot as plt
from scipy import special
import mpl_toolkits.mplot3d.axes3d as axes3d

n = 15

thetas, phis = np.linspace(0, np.pi, n), np.linspace(0, 2*np.pi, n)
THETA, PHI = np.meshgrid(thetas, phis)

beta = 0.5
gamma = np.pi/6

def R(theta, phi):
    temp1 = 3 * np.cos(theta) * np.cos(theta) - 1
    temp2 = np.sin(theta) * np.sin(theta) * np.cos(2*phi)
    temp = np.cos(gamma) * temp1 + np.sqrt(3) * np.sin(gamma) * temp2
    return 1 + beta * np.sqrt(5/16/np.pi) * temp

R = R(THETA, PHI)
X = R * np.sin(THETA) * np.cos(PHI)
Y = R * np.sin(THETA) * np.sin(PHI)
Z = R * np.cos(THETA)
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')
plot = ax.plot_surface(
    X, Y, Z, rstride=1, cstride=1, cmap=plt.get_cmap('jet'),
    linewidth=0, antialiased=False, alpha=0.5 )
#plt.title(r"$\sum_m|Y^m_4|^2$")
# below are codes copied from stackoverflow, to make the scaling correct
max_range = np.array([X.max()-X.min(), Y.max()-Y.min(), Z.max()-Z.min()]).max() / 2.0
mid_x = (X.max()+X.min()) * 0.5
mid_y = (Y.max()+Y.min()) * 0.5
mid_z = (Z.max()+Z.min()) * 0.5
ax.set_xlim(mid_x - max_range, mid_x + max_range)
ax.set_ylim(mid_y - max_range, mid_y + max_range)
ax.set_zlim(mid_z - max_range, mid_z + max_range)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$z$')
ax.set_aspect('equal')
#help(ax)
ax.view_init(elev=30,azim=0) #调节视角，elev指向上(z方向)旋转的角度，azim指xy平面内旋转的角度
#plt.legend(r"$|Y_{lm}|^2$")
plt.show()

长椭球是一个维度特别长，另外两个维度有旋转对称性：

扁椭球是一个维度特别短，另外两个维度有旋转对称性：

三轴形变是三个轴都不同，分为长、中、短：