雙縫干涉相關(guān)代碼
相當(dāng)于在粒子位矢波函數(shù)上張量積了兩個(gè)雙態(tài)空間，方法還是 time splitting 方法。
很多地方?jīng)]怎么精細(xì)處理，比較臃腫，不過(guò)能跑就行
L = 15; % Interval Length
N = 400; % 設(shè)置 N*N 格點(diǎn)
x = linspace(-L/2,L/2,N)'; % Coordinate vector
dx = x(2) - x(1); % Coordinate step
NT = 20000; % No. of time steps
TF = 10; T = linspace(0,TF,NT); % Time vector
dT = T(2)-T(1); % Time step
%**************************
% 二維單位矩陣和 Pauli x 矩陣
I2 = speye(2,2);
sigmax = sparse(2,2);
sigmax(1,2) = 1; sigmax(2,1) = 1;
%**************************
% 位矢空間相關(guān)算符
e1 = ones(N,1);
Rx = kron(x',e1);
Ry = Rx';
InArea = heaviside(-2-Rx);
InArea = reshape(InArea,1,N*N); % 篩出 Rx<-2 范圍內(nèi)的振幅
%**************************
% 勢(shì)能
V = 1000;
Ud = V*exp(-25*(Rx+2).^2)...
? ?.*heaviside(0.75*pi-Ry).*heaviside(Ry+0.75*pi)...
? ? ? ?.*cos(Ry*2).^2+...
? ? ? ?V*exp(-25*(Rx-0.5*L).^2)+V*exp(-25*(Rx+0.5*L).^2)+...
? ? ? ?V*exp(-25*(Ry-0.5*L).^2)+V*exp(-25*(Ry+0.5*L).^2);% 雙縫與邊界
% 演化算符相關(guān)
EUdx = kron(...
? ?reshape(exp(-1i*0.5*Ud*dT),1,N*N),...
? ?[1;1;1;1]);
? ?
%**************************
% probe
g = 90; % coupling
% coupling region
UpR = g*exp(-10*((Rx-0.03*L+2).^2+(Ry-0.25*pi).^2)); % right
UpL = g*exp(-10*((Rx-0.03*L+2).^2+(Ry+0.25*pi).^2)); % left
% 演化算符相關(guān)
SUpRx = sin(0.5*UpR*dT);
SUpRx = kron(reshape(SUpRx,1,N*N),[1;1;1;1]);
SUpRd = kron(I2,sigmax);
CUpRx = cos(0.5*UpR*dT);
CUpRx = kron(reshape(CUpRx,1,N*N),[1;1;1;1]);
SUpLx = sin(0.5*UpL*dT);
SUpLx = kron(reshape(SUpLx,1,N*N),[1;1;1;1]);
SUpLd = kron(sigmax,I2);
CUpLx = cos(0.5*UpL*dT);
CUpLx = kron(reshape(CUpLx,1,N*N),[1;1;1;1]);
%繪制探測(cè)區(qū)域的邊界
circxR = cos(linspace(0,2*pi,50))/sqrt(10)+0.03*L-2;
circyR = sin(linspace(0,2*pi,50))/sqrt(10)+0.25*pi;
circzR = zeros(50);
circxL = cos(linspace(0,2*pi,50))/sqrt(10)+0.03*L-2;
circyL = sin(linspace(0,2*pi,50))/sqrt(10)-0.25*pi;
circzL = zeros(50);
%**************************
% 動(dòng)量相關(guān)算符
e1 = ones(N,1);
K = [2*pi*(-N/2:-1)/L,2*pi*(0:N/2-1)/L]; % 動(dòng)量范圍
Kx = kron(K,e1);
Ky = Kx';
Ks = (Kx.^2+Ky.^2);
% 時(shí)間演化相關(guān)
UMmtk = exp(-0.5*1i*dT*(Ks));
%**************************
% 初態(tài)
a = 1;% 控制波包大小
ko = 30;% 動(dòng)量中心
% 位置中心
xo = -4;
yo = 0;
psi0 = exp(-((Rx-xo).^2+(Ry-yo).^2)/(a^2)).*exp(1i*Rx*ko); % Gauss 波包
psi0 = psi0/sqrt(sum(sum(conj(psi0).*psi0))); % 歸一化
psi2 = kron(reshape(psi0,1,N*N),[1;0;0;0]); % 化為方便處理的形式
%**************************
% 時(shí)間演化與繪圖
%繪圖相關(guān)設(shè)置
fig=figure;
fig.Position = [0 0 2000 1000];
ax = axes(fig);
h = 0.5e-3;
xlim manual
ylim manual
zlim manual
CO1=zeros(N,N,3);
CO1(:,:,1) = 0.2*ones(N,N);
CO1(:,:,3) = 0.75*ones(N,N);
CO2 = CO1;
CO2(:,:,3) = 1*ones(N,N);
for t = 1:NT
? ?% Strang 方法計(jì)算時(shí)間演化
? ?% 依次作用 right probe, left probe, potential 相關(guān)算符
? ?% （不過(guò)三者其實(shí)對(duì)易，具體順序不重要）
? ?psi2 = CUpRx.*psi2-1i*SUpRx.*(SUpRd*psi2);
? ?psi2 = CUpLx.*psi2-1i*SUpLx.*(SUpLd*psi2);
? ?psi2 = EUdx.*psi2;
? ?% 計(jì)算動(dòng)量波函數(shù)，十分臃腫，放棄思考
? ?% 并作用動(dòng)量相關(guān)算子
? ?psik1 = fftshift(fft2(reshape(psi2(1,:),N,N)));
? ?psik1 = UMmtk.*psik1;
? ?psik2 = fftshift(fft2(reshape(psi2(2,:),N,N)));
? ?psik2 = UMmtk.*psik2;
? ?psik3 = fftshift(fft2(reshape(psi2(3,:),N,N)));
? ?psik3 = UMmtk.*psik3;
? ?psik4 = fftshift(fft2(reshape(psi2(4,:),N,N)));
? ?psik4 = UMmtk.*psik4;
? ?% 再逆變換回位矢波函數(shù)
? ?psi2(1,:) = reshape(ifft2(ifftshift(psik1)),1,N*N);
? ?psi2(2,:) = reshape(ifft2(ifftshift(psik2)),1,N*N);
? ?psi2(3,:) = reshape(ifft2(ifftshift(psik3)),1,N*N);
? ?psi2(4,:) = reshape(ifft2(ifftshift(psik4)),1,N*N);
% 再次依次作用 right probe, left probe, potential 相關(guān)算符
? ?psi2 = CUpRx.*psi2-1i*SUpRx.*(SUpRd*psi2);
? ?psi2 = CUpLx.*psi2-1i*SUpLx.*(SUpLd*psi2);
? ?psi2 = EUdx.*psi2;
? ?
? ?if rem(t,10) % 控制畫(huà)圖頻率
? ? ? ?if t~=1
? ? ? ? ? ?continue
? ? ? ?end
? ?end
? ?
? ?rho = reshape(sum(abs(psi2).^2,1),N,N); % 總位矢分布密度
? ?% 曲面繪制與著色
? ?CO1(:,:,1) = 0.2*ones(N,N)+10*ones(N,N).*sqrt(rho/h)./((1+50*rho/h));
? ?CO1(:,:,2) = 15*ones(N,N).*rho/h;
? ?CO1(:,:,3) = 0.75*ones(N,N)+ones(N,N).*rho.^2/h^2;
? ?s01 = surf(ax,x,x,rho+0.01*h,CO1,'FaceAlpha',0.4);
? ?shading interp
? ?s01.EdgeColor = 'none';
? ?
? ?zlim(ax,[-0.015*h 3*h]);
? ?xlim(ax,[-0.5*L 0.5*L]);
? ?ylim(ax,[-0.5*L 0.5*L]);
? ?
? ?hold on
? %時(shí)間
? ?text(ax,-0.4*L,0.4*L,2.1*h,...
? ? ? ?['$T=',num2str(dT*t,'%.3f'),'$'],...
? ? ? ?'Interpreter','latex','FontSize',20);
? ?% 兩個(gè)探針在 2*2=4 個(gè)狀態(tài)上的分布
? ?None = sum(abs(psi2(1,:)).^2);
? ?Left = sum(abs(psi2(2,:)).^2);
? ?Right = sum(abs(psi2(3,:)).^2);
? ?Both = sum(abs(psi2(4,:)).^2);
? ?text(ax,-0.3*L,-0.5*L,2.3*h,...
? ? ? ?['$None:\ ',num2str(None,'%.3f'),'$',newline,...
? ? ? ?'$Left:\ ',num2str(Left,'%.3f'),'$',newline, ...
? ? ? ?'$Right:\ ',num2str(Right,'%.3f'),'$',newline, ...
? ? ? ?'$Both:\ ',num2str(Both,'%.3f'),'$'],...
? ? ? ?'Interpreter','latex','FontSize',20);
? ?% 耦合強(qiáng)度
? ?text(ax,-0.4*L,0.4*L,1.8*h,...
? ? ? ?['$g =',num2str(g,'%.0f'),'$'],...
? ? ? ?'Interpreter','latex','FontSize',20);
? ?% 直接計(jì)算反射率
? ?Rfl = sum(abs(InArea.*psi2(1,:)).^2)+...
? ?sum(abs(InArea.*psi2(2,:)).^2)+...
? ?sum(abs(InArea.*psi2(3,:)).^2)+...
? ?sum(abs(InArea.*psi2(4,:)).^2);
? ?
? ?text(ax,-0.3*L,-0.5*L,3*h,...
? ? ? ?['Reflect:$\ ',num2str(Rfl,'%.3f'),'$'],...
? ? ? ?'Interpreter','latex','FontSize',20);
? ?% 用顏色顯示勢(shì)能
? ?imagesc(ax,[-0.5*L,0.5*L],[-0.5*L,0.5*L],...
? ? ? ?10*h*Ud/V,'AlphaData',0.15);
? ?% 標(biāo)注探測(cè)區(qū)域
? ?plot3(ax,circxR,circyR,circzR,'r');
? ?plot3(ax,circxL,circyL,circzL,'b');
? ?
? ?hold off
? ?% 繪圖參數(shù)
? ?colormap(hsv)
? ?caxis([-pi*h pi*h])
? ?pbaspect([1 1 0.5]);
? ?campos([25 -10 5*h])
? ?camtarget([0 0 1.5*h])
? ?camva(23)
? ?
? ?drawnow
end
