氫原子相關(guān)matlab代碼。

用于確認(rèn)躍遷類(lèi)型的代碼補(bǔ)充到后續(xù)專(zhuān)欄文章。

%生成波函數(shù)

% matlab 自帶的球坐標(biāo)、直角坐標(biāo)轉(zhuǎn)換與計(jì)算等值面的函數(shù)好像不太配套

% 所以就額外寫(xiě)了坐標(biāo)轉(zhuǎn)換的函數(shù)

rho = @(x,y,z) sqrt(x.^2+y.^2+z.^2);% 由x,y,z給出球坐標(biāo)中的r

% 球坐標(biāo)中的theta、phi也被定義為函數(shù)。

L = 20; % 立方形區(qū)域邊長(zhǎng)

Num = 201;% 單邊長(zhǎng)格點(diǎn)數(shù)，總格點(diǎn)數(shù)為 Num^3

x=-L:(2*L/(Num-1)):L;

Nx = length(x);?

[X,Y,Z] = meshgrid(x);% 三維格點(diǎn)

% 變成球坐標(biāo)

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); % 生成波函數(shù)

% state 2

n2 = 2;

l2 = 0;

m2 = 0;

E2 = -1/(8*pi*n2^2);

psi2 = Hstate(R,T,P,n2,l2,m2);

%驗(yàn)證歸一化

% 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>的歸一化波函數(shù)

Nx = size(R);?

Nx = Nx(1);

% 球諧函數(shù)

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

%徑向波函數(shù)

% Rd = laguerreL(n-l-1,2*l+1,2*R/n);%matlab自帶函數(shù)太慢

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)% 球坐標(biāo)的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) % 球坐標(biāo)的phi

ph = atan2(y,x);

end

function L = mlaguerre(n,p,x) % 根據(jù)網(wǎng)上的建議，直接用級(jí)數(shù)法計(jì)算Laguerre多項(xiàng)式

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

% 畫(huà)圖

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; % 時(shí)間間隔

Tt = 10000; % 最大時(shí)長(zhǎng)

for t=0:dt:Tt

%? ? ?t 時(shí)刻的波函數(shù)

? ? 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

????

????% 計(jì)算并繪制等值面

? ? fv =...

? ? ? ? isosurface(X,Y,Z,abs(psi).^2*(nm^6),...

? ? ? ? 0.006, real(exp(1i*angle(psi))) ); % 等值面的值為0.006（可更改），按照?exp(1i*angle(psi)) 的實(shí)部著色

? ?Sf = patch(ax,...

? ? ? ? fv,...

? ? ? ? 'FaceColor','interp','EdgeColor','none','facealpha',0.3);

????% 3維圖參數(shù)

? ? 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