氢原子波函数相关代码
氢原子相关matlab代码。
用于确认跃迁类型的代码补充到后续专栏文章。
%生成波函数
% matlab 自带的球坐标、直角坐标转换与计算等值面的函数好像不太配套
% 所以就额外写了坐标转换的函数
rho = @(x,y,z) sqrt(x.^2+y.^2+z.^2);% 由x,y,z给出球坐标中的r
% 球坐标中的theta、phi也被定义为函数。
L = 20; % 立方形区域边长
Num = 201;% 单边长格点数，总格点数为 Num^3
x=-L:(2*L/(Num-1)):L;
Nx = length(x);
[X,Y,Z] = meshgrid(x);% 三维格点
% 变成球坐标
R = rho(X,Y,Z);
T = theta(X,Y,Z);
P = phi(X,Y,Z);
% state 1
n1 = 3;
l1 = 2;
m1 = 2;
E1 = -1/(8*pi*n1^2); % energy
psi1 = Hstate(R,T,P,n1,l1,m1); % 生成波函数
% state 2
n2 = 2;
l2 = 0;
m2 = 0;
E2 = -1/(8*pi*n2^2);
psi2 = Hstate(R,T,P,n2,l2,m2);
%验证归一化
% Dim = Num^3;
% psi01=reshape(psi1,[Dim,1]);
% psi02=reshape(psi2,[Dim,1]);
% NormalL = 1/(2*L/(Num-1))^3
% psi01'*psi01/NormalL
% psi02'*psi02/NormalL
function psi = Hstate(R,T,P,n,l,m) %生成|nlm>的归一化波函数
Nx = size(R);
Nx = Nx(1);
% 球谐函数
C1 = sqrt(((2*l+1)*factorial(l-abs(m)))/(4*pi*factorial(l+abs(m))));
lgd = legendre(l,cos(T));
if l
Ylm = C1*reshape(lgd(abs(m)+1,:,:,:),Nx,Nx,Nx)...
.*exp(1i*m*P);
else
Ylm = C1*lgd...
.*exp(1i*m*P);
end
%径向波函数
% Rd = laguerreL(n-l-1,2*l+1,2*R/n);%matlab自带函数太慢
Rd = mlaguerre(n-l-1,2*l+1,2*R/n).*power(2*R/n,l).*exp(-R/n);
C2 = power(2/n,1.5)*sqrt((factorial(n-l-1))/(2*n*factorial(l+n)));
psi = C2*Rd.*Ylm;
end
function th = theta(x,y,z)% 球坐标的theta
r = sqrt(x.^2+y.^2+z.^2);
c = z./r;
c(isinf(c))=0;
c(isnan(c))=0;
th = acos(c);
end
function ph = phi(x,y,z) % 球坐标的phi
ph = atan2(y,x);
end
function L = mlaguerre(n,p,x) % 根据网上的建议，直接用级数法计算Laguerre多项式
ret = 0;
for i=0:n
sum1 = ret + (power(-1,i)* (factorial(n+p)/(factorial(n-i) * factorial(p+i) * factorial(i)))*power(x,i));
ret=sum1;
end
L=ret;
end
Om = 0.2*abs(E1-E2); % Rabi freq
% 画图
fig = figure;
fig.Position = [0 0 2000 1000];
ax = axes(fig);
grid on
xticks(-L:L/5:L)
yticks(-L:L/5:L)
zticks(-L:L/5:L)
xlim manual
ylim manual
zlim manual
nm = min(n1,n2);
Txt0 = text(ax,-0.8*L,0.8*L,0.8*L,...
['$|',num2str(n1,'%.0f'),num2str(l1,'%.0f'),...
num2str(m1,'%.0f'),'\rangle\ \leftrightarrow\ ',...
'\ |',num2str(n2,'%.0f'),num2str(l2,'%.0f'),...
num2str(m2,'%.0f'),'\rangle$'],...
'Interpreter','latex','FontSize',20);
dt = 10; % 时间间隔
Tt = 10000; % 最大时长
for t=0:dt:Tt
% t 时刻的波函数
psi = cos(0.5*Om*t)*exp(-1i*E1*t)*psi1+...
1i*sin(0.5*Om*t)*exp(-1i*E2*t)*psi2;
if t~=0
delete(Sf);
delete(Txt1);
end
% 计算并绘制等值面
fv =...
isosurface(X,Y,Z,abs(psi).^2*(nm^6),...
0.006, real(exp(1i*angle(psi))) ); % 等值面的值为0.006（可更改），按照 exp(1i*angle(psi)) 的实部着色
Sf = patch(ax,...
fv,...
'FaceColor','interp','EdgeColor','none','facealpha',0.3);
% 3维图参数
view(3)
xlim([-L,L])
ylim([-L,L])
zlim([-L,L])
pbaspect([1 1 1]);
campos([-1.5*L -3*L 0.5*L])
camtarget([0 0 0])
camva(36)
Txt1 = text(-0.8*L,0.8*L,0.6*L,...
['$T=',num2str(t,'%2.1f'),'$'],...
'Interpreter','latex','FontSize',20);
% F = getframe(fig);
if t>5710
% close(v)
break
end
drawnow
end
