在知乎看到了 [王欢]的[用 Python+SciPy 可视化定态氢原子波函数（一）] [王大可]的[利用Mathematica将定态氢原子波函数可视化]
手痒，想用MATLAB也实现一下，David Griffiths 所著的 Introduction to Quantum Mechanics 中: ψ n , l , m ( r , θ , ϕ ) = R n l ( r ) Y l m ( θ , ϕ ) \psi_{n, l, m}(r,\theta, \phi)=R_{nl}(r)Y_{l}^{m}(\theta, \phi) ψn,l,m(r,θ,ϕ)=Rnl(r)Ylm(θ,ϕ) 即，函数由径向波函数及球形谐函数两部分组成，展开一下径向波函数是这样的： ψ n , l , m ( r , θ , ϕ ) = ( 2 n a ) 3 ( n − l − 1 ) ! 2 n ( n + l ) ! e − r / n a ( 2 r n a ) l [ L n − l − 1 2 l + 1 ( 2 r n a ) ] Y l m ( θ , ϕ ) \psi_{n, l, m}(r,\theta, \phi)=\sqrt{\left(\frac{2}{n a}\right)^{3} \frac{(n-l-1) !}{2 n(n+l) !}} e^{-r / n a}\left(\frac{2 r}{n a}\right)^{l}\left[L_{n-l-1}^{2 l+1}\left(\frac{2 r}{n a}\right)\right] Y_{l}^{m}(\theta, \phi) ψn,l,m(r,θ,ϕ)=(na2)32n(n+l)!(n−l−1)! e−r/na(na2r)l[Ln−l−12l+1(na2r)]Ylm(θ,ϕ)
其中 a a a 为玻尔半径，数值约为5.2917721092(17)×10^-11米，本文仅做可视化，方便起见，将 a a a 数值设置为1，并省略归一化系数，因此我们实际上本文绘制的是如下函数的图像：
ψ n , l , m ( r , θ , ϕ ) = e − r / n ( 2 r n ) l [ L n − l − 1 2 l + 1 ( 2 r n ) ] Y l m ( θ , ϕ ) \psi_{n, l, m}(r,\theta, \phi)=e^{-r / n}\left(\frac{2 r}{n}\right)^{l}\left[L_{n-l-1}^{2 l+1}\left(\frac{2 r}{n}\right)\right] Y_{l}^{m}(\theta, \phi) ψn,l,m(r,θ,ϕ)=e−r/n(n2r)l[Ln−l−12l+1(n2r)]Ylm(θ,ϕ)
其中非常棒的是 l a g u e r r e laguerre laguerre 多项式 L n − l − 1 2 l + 1 ( 2 r n ) L_{n-l-1}^{2 l+1}\left(\frac{2 r}{n}\right) Ln−l−12l+1(n2r)，MATLAB中有自带的函数 laguerreL可以使用，但是球谐函数(SphericalHarmonicY) Y l m ( θ , ϕ ) Y_{l}^{m}(\theta, \phi) Ylm(θ,ϕ)，并不是MATLAB的自带函数，不过好在，当 m ≥ 0 m\ge0 m≥0时球谐函数的展开式：
Y
l
m
(
θ
,
ϕ
)
=
(
−
1
)
m
(
2
l
+
1
)
4
π
(
l
−
∣
m
∣
)
!
(
l
+
∣
m
∣
)
!
e
i
m
ϕ
P
l
m
(
cos
θ
)
Y_{l}^{m}(\theta, \phi)=(-1)^m \sqrt{\frac{(2 l+1)}{4 \pi} \frac{(l-|m|) !}{(l+|m|) !}} e^{i m \phi} P_{l}^{m}(\cos \theta)
Ylm(θ,ϕ)=(−1)m4π(2l+1)(l+∣m∣)!(l−∣m∣)!
eimϕPlm(cosθ) 其中
L
e
g
e
n
d
r
e
Legendre
Legendre 函数
P
l
m
(
cos
θ
)
P_{l}^{m}(\cos \theta)
Plm(cosθ): 
虽然，也很复杂。。但是耐不住MATLAB自带这个函数啊，甚至MATLAB还给出了如何通过该函数生成球谐函数的实例，简直太贴心了，爱了爱了：
dx=pi/60;
col=0:dx:pi;
az=0:dx:2*pi;
[phi,theta]=meshgrid(az,col);
% 计算 l=3 的网格上的 P_{l}^{m}(\cos \theta)
l=3;
Plm=legendre(l,cos(theta));
% 由于 legendre 为 m 的所有值计算答案，因此 Plm 会包含一些额外的函数值。
% 提取 m=2 的值并丢弃其余值。
% 使用 reshape 函数将结果定向为与 phi 和 theta 具有相同大小的矩阵。
m=2;
if l~=0
Plm=reshape(Plm(m+1,:,:),size(phi));
end
% 计算Y_3^2的球谐函数值
a=(2*l+1)*factorial(l-m);
b=4*pi*factorial(l+m);
C=sqrt(a/b);
Ylm=C.*Plm.*exp(1i*m*phi);
% 球面坐标转换为笛卡尔坐标并绘图
[Xm,Ym,Zm]=sph2cart(phi,pi/2-theta,abs(real(Ylm)));
surf(Xm,Ym,Zm)
title('$Y_3^2$ spherical harmonic','interpreter','latex')
万事俱备开始编程：
工具函数function psi=HydWave(n,l,m,x,y,z)
% @author:slandarer
% 由坐标点计算向量模长及角度
r=vecnorm([x(:),y(:),z(:)]')';
theta=atan2(vecnorm([x(:),y(:)]')',z(:));
phi=atan2(y(:),x(:));
% 恢复矩阵型状
r=reshape(r,size(x));
theta=reshape(theta,size(x));
phi=reshape(phi,size(x));
% 利用MATLAB自带legendre函数计算球谐函数
Plm=legendre(l,cos(theta));
if l~= 0
Plm=reshape(Plm(m+1,:,:),size(phi));
end
C=sqrt(((2*l+1)*factorial(l-m))/(4*pi*factorial(l+m)));
Ylm=C.*Plm.*exp(1i*m*phi);
% laguerreL函数部分计算
Lag=laguerreL(n-l-1,2*l+1,2.*r./n);
% 整合起来
psi=exp(-r./n).*(2.*r./n).^l.*Lag.*Ylm;
psi=real(conj(psi).*psi);
end
[X,Z]=meshgrid(linspace(-30,30,120));
Y=zeros(size(X));
psi=HydWave(4,2,0,X,Y,Z);
surf(X,Z,psi,'EdgeColor','none')
view(2);caxis([0,2])
function HydWaveDemo
% @author:slandarer
% 构造一个好看的颜色映射实属不易
map=[0.1294 0.0549 0.1725;0.2196 0.1608 0.2902;0.3882 0.1804 0.4941;
0.4392 0.1922 0.4706;0.5333 0.2235 0.4392;0.6471 0.2588 0.3686;
0.7137 0.2745 0.3294;0.7725 0.3059 0.2902;0.8510 0.3725 0.2275;
0.9137 0.4196 0.1804;0.9608 0.5020 0.2000;0.9765 0.5529 0.2078;
0.9804 0.6431 0.2549;0.9843 0.6627 0.2706;0.9765 0.7176 0.3412;
0.9765 0.7686 0.4000;0.9765 0.8118 0.4902;0.9725 0.8510 0.5961;
0.9882 0.9020 0.6667;1.0000 0.9451 0.8431;1.0000 0.9961 0.9804;
1.0000 1.0000 1.0000];
Xi=1:size(map,1);Xq=linspace(1,size(map,1),800);
map=[interp1(Xi,map(:,1),Xq,'linear')',...
interp1(Xi,map(:,2),Xq,'linear')',...
interp1(Xi,map(:,3),Xq,'linear')'];
% 情况信息列表
condList=[1,2,0,0,-10,10,0.01;
2,3,0,0,-20,20,0.01;
6,2,1,0,-12,12,0.03;
7,3,1,0,-20,20,0.08;
8,3,1,1,-20,20,0.08;
11,2,1,1,-12,12,0.03;
12,3,2,0,-23,23,0.35;
13,3,2,1,-23,23,0.35;
14,3,2,2,-23,23,0.35;
16,4,0,0,-36,36,0.008;
17,4,1,0,-30,30,0.2;
18,4,1,1,-30,30,0.1;
19,4,2,0,-30,30,0.95;
20,4,2,1,-30,30,0.95;
21,4,2,2,-30,30,0.95;
22,4,3,0,-35,35,6;
23,4,3,1,-35,35,6;
24,4,3,2,-35,35,6;
25,4,3,3,-35,35,1.8];
fig=gcf;
fig.Color=[0,0,0];
set(fig,'InvertHardCopy','off');
% 文本信息
axh=subplot(5,5,[3,4,5]);
axh.XLim=[0,3];
axh.YLim=[0,1];
axh.XTick=[];
axh.YTick=[];
axh.XColor=[0,0,0];
axh.YColor=[0,0,0];
axh.Color=[0 0 0];
text(axh,.04,0.8,'Hydrogen Electron Orbita','Color',[1,1,1].*.98,...
'FontName','cambria','FontWeight','bold','FontSize',15)
text(axh,.04,0.2,'$\psi_{n, l, m}=e^{-r / n}\left(\frac{2 r}{n}\right)^{l}\left[L_{n-l-1}^{2 l+1}\left(\frac{2 r}{n}\right)\right] Y_{l}^{m}(\theta, \phi)$',...
'Color',[1,1,1].*.98,'FontSize',12,'Interpreter','latex')
% 文本信息2
axh2=subplot(5,5,[9,10]);
axh2.XLim=[0,2];
axh2.YLim=[0,1];
axh2.XTick=[];
axh2.YTick=[];
axh2.XColor=[0,0,0];
axh2.YColor=[0,0,0];
axh2.Color=[0 0 0];
text(axh2,.04,0.8,'Not normalized by:','Color',[1,1,1].*.98,...
'FontName','cambria','FontSize',13)
text(axh2,.04,0.2,'$\sqrt{\left(\frac{2}{n a}\right)^{3} \frac{(n-l-1) !}{2 n(n+l) !}}$','Color',[1,1,1].*.98,...
'FontName','cambria','FontSize',12,'Interpreter','latex')
% 循环绘图
for i=1:size(condList,1)
ax=subplot(5,5,condList(i,1));
% 绘制密度分布
[X,Z]=meshgrid(linspace(condList(i,5),condList(i,6),120));
Y=zeros(size(X));
psi=HydWave(condList(i,2),condList(i,3),condList(i,4),X,Y,Z);
surf(X,Z,psi,'EdgeColor','none')
caxis(ax,[0,condList(i,7)]);colormap(map);view(2)
axis tight
% 坐标区域修饰
ax.XLim=[condList(i,5),condList(i,6)];
ax.YLim=[condList(i,5),condList(i,6)];
ax.DataAspectRatio=[1,1,1];
ax.XTick=[];
ax.YTick=[];
ax.XColor=[1,1,1].*.4;
ax.YColor=[1,1,1].*.4;
ax.LineWidth=1.5;
ax.Box='on';
ax.FontName='cambria';
ax.XLabel.String=['(',num2str(condList(i,2)),',',...
num2str(condList(i,3)),',',...
num2str(condList(i,4)),')'];
drawnow
end
saveas(fig,'result.png')
end
