双缝干涉相关代码

相当于在粒子位矢波函数上张量积了两个双态空间，方法还是 time splitting 方法。

很多地方没怎么精细处理，比较臃肿，不过能跑就行

L = 15; % Interval Length

N = 400; % 设置 N*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;

%**************************

% 位矢空间相关算符

e1 = ones(N,1);

Rx = kron(x',e1);

Ry = Rx';

InArea = heaviside(-2-Rx);

InArea = reshape(InArea,1,N*N); % 筛出 Rx<-2 范围内的振幅

%**************************

% 势能

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);% 双缝与边界

% 演化算符相关

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

% 演化算符相关

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]);

%绘制探测区域的边界

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);

%**************************

% 动量相关算符

e1 = ones(N,1);

K = [2*pi*(-N/2:-1)/L,2*pi*(0:N/2-1)/L]; % 动量范围

Kx = kron(K,e1);

Ky = Kx';

Ks = (Kx.^2+Ky.^2);

% 时间演化相关

UMmtk = exp(-0.5*1i*dT*(Ks));

%**************************

% 初态

a = 1;% 控制波包大小

ko = 30;% 动量中心

% 位置中心

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]); % 化为方便处理的形式

%**************************

% 时间演化与绘图

%绘图相关设置

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 方法计算时间演化

   % 依次作用 right probe, left probe, potential 相关算符

   % （不过三者其实对易，具体顺序不重要）

   psi2 = CUpRx.*psi2-1i*SUpRx.*(SUpRd*psi2);

   psi2 = CUpLx.*psi2-1i*SUpLx.*(SUpLd*psi2);

   psi2 = EUdx.*psi2;

   % 计算动量波函数，十分臃肿，放弃思考

   % 并作用动量相关算子

   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;

   % 再逆变换回位矢波函数

   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 相关算符

   psi2 = CUpRx.*psi2-1i*SUpRx.*(SUpRd*psi2);

   psi2 = CUpLx.*psi2-1i*SUpLx.*(SUpLd*psi2);

   psi2 = EUdx.*psi2;

   

   if rem(t,10) % 控制画图频率

       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

  %时间

   text(ax,-0.4*L,0.4*L,2.1*h,...

       ['$T=',num2str(dT*t,'%.3f'),'$'],...

       'Interpreter','latex','FontSize',20);

   % 两个探针在 2*2=4 个状态上的分布

   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);

   % 耦合强度

   text(ax,-0.4*L,0.4*L,1.8*h,...

       ['$g =',num2str(g,'%.0f'),'$'],...

       'Interpreter','latex','FontSize',20);

   % 直接计算反射率

   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);

   % 用颜色显示势能

   imagesc(ax,[-0.5*L,0.5*L],[-0.5*L,0.5*L],...

       10*h*Ud/V,'AlphaData',0.15);

   % 标注探测区域

   plot3(ax,circxR,circyR,circzR,'r');

   plot3(ax,circxL,circyL,circzL,'b');

   

   hold off

   % 绘图参数

   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