二維無質(zhì)量 Dirac 型粒子的相關(guān)代碼
二維無質(zhì)量 Dirac 型粒子相關(guān)模擬。
L = 20; % Interval Length
N = 600; % 設(shè)置 N*N 格點
x = linspace(-L/2,L/2-L/N,N)'; % Coordinate vector
dx = x(2) - x(1); % Coordinate step
NT = 50000; % No. of time steps
TF = 50; T = linspace(0,TF,NT); % Time vector
dT = T(2)-T(1); % Time step
%**************************
% 設(shè)置勢能
% 可按需調(diào)整
[X,Y]=meshgrid(x,x);
V=150;?
lth = 2;
Ub = V*heaviside(X).*heaviside(lth-X);
% Ub = V*exp(-X.^2/lth^2);
UVh = exp(-1i*Ub*dT/2);
% 動量相關(guān)
v0 = 10; % 無質(zhì)量Dirac粒子的速度
e1 = ones(N,1);?
K = [2*pi*(-N/2:-1)/L,2*pi*(0:N/2-1)/L]; % 動量范圍
Kx = kron(K,e1);
Ky = Kx';
Ka = sqrt(Kx.^2+Ky.^2);
ThK = angle(Kx+1i*Ky);
%**************************
% 初態(tài)設(shè)置
ao = 1.8;% 控制動量空間波包大小
% 動量中心
ko = 12;
Th = 0;
kxo = ko*cos(Th);
kyo = ko*sin(Th);
% 位矢中心
xo = -2.5;
yo = 0;
% first component?in k space
psi01 = exp(-((Kx-kxo).^2+100*(Ky-kyo).^2)/(ao^2))...
? ? .*exp(-1i*(xo*Kx+yo*Ky)).*(cos(ThK)-1i*sin(ThK))/sqrt(2);
% second component
psi02 = exp(-((Kx-kxo).^2+100*(Ky-kyo).^2)/(ao^2))...
? ? .*exp(-1i*(xo*Kx+yo*Ky)).*1/sqrt(2);
% 位矢空間波函數(shù)，未歸一
psi0 = [reshape( fftshift(ifft2(fftshift(psi01))),N*N,1 )...
? ? ;reshape( fftshift(ifft2( fftshift(psi02) )),N*N,1)];
psi0 = psi0/sqrt(psi0'*psi0); %歸一化
% 初態(tài)的雙分量位矢空間波函數(shù)
psi21 = reshape(psi0(1:end/2),N,N);
psi22 = reshape(psi0(0.5*end+1:end),[N,N]);
%**************************
% 時間演化和畫圖
%建立figure和axes
fig=figure;
fig.Position = [0 0 2000 1000];
ax = axes(fig);
h = 20e-5; % z軸相關(guān)尺度
xlim manual
ylim manual
zlim manual
%關(guān)于顏色
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);
%相位平面位置
Z0 = -1.5*h*ones(N,N);
for t = 1:NT
%時間演化，Strang 方法
? ? psi21 = full(UVh.*psi21);
? ? psi22 = full(UVh.*psi22);
? ? psik1 = fftshift(fft2(psi21));
? ? psik2 = fftshift(fft2(psi22));
? ??
? ? psik10 = psik1;
? ? psik1 = cos(Ka*v0*dT).*psik1-1i*exp(-1i*ThK).*sin(Ka*v0*dT).*psik2;
? ? psik2 = cos(Ka*v0*dT).*psik2-1i*exp(1i*ThK).*sin(Ka*v0*dT).*psik10;
? ??
? ? psi21 = full(ifft2(ifftshift(psik1)));
? ? psi22 = full(ifft2(ifftshift(psik2)));
? ??
? ? psi21 = full(UVh.*psi21);
? ? psi22 = full(UVh.*psi22);
%以上是一步（時間間隔為 dT）的時間演化
? ? % 每隔10步畫圖
? ? if rem(t,10)
? ? ? ? if t~=1
? ? ? ? ? ? continue
? ? ? ? end
? ? end
? ? rho1 = (abs(psi21).^2); % 第一個分量的密度
% rgb顏色，常數(shù)項是本底顏色，根號項凸顯小振幅，一次和平方項給大振幅產(chǎn)生顏色漸變
? ? CO1(:,:,1) = 0.2*ones(N,N)+10*ones(N,N).*sqrt(rho1/h)./((1+50*rho1/h));
? ? CO1(:,:,2) = 15*ones(N,N).*rho1/h;
? ? CO1(:,:,3) = 0.75*ones(N,N)+ones(N,N).*rho1.^2/h^2;
? ? s01 = surf(ax,x,x,rho1+0.1*h,CO1,'FaceAlpha',0.4);?
? ? shading interp
? ? s01.EdgeColor = 'none';
? ??
? ? zlim(ax,[-1.5*h 2.1*h]);
? ? xlim(ax,[-0.5*L 0.5*L]);
? ? ylim(ax,[-0.5*L 0.5*L]);
? ??
? ? hold on
? ? rho2 = (abs(psi22).^2); % 第二個分量的密度
? ? CO2(:,:,1) = 0.2*ones(N,N)+10*ones(N,N).*sqrt(rho2/h)./((1+50*rho2/h));
? ? CO2(:,:,2) = 15*ones(N,N).*rho2/h;
? ? CO2(:,:,3) = 1*ones(N,N)+ones(N,N).*rho2.^2/h^2;
? ? s02 = surf(ax,x,x,-rho2-0.1*h,1-CO2,'FaceAlpha',0.4);?
? ? shading interp
? ? s02.EdgeColor = 'none';
? ??
? ? text(ax,-0.4*L,0.4*L,2.1*h,...
? ? ? ? ['$T=',num2str(dT*t,'%.2f'),'$'],...
? ? ? ? 'Interpreter','latex','FontSize',20); %時間
? ??
? ? text(ax,0.35*L,0.35*L,2.05*h,'方勢壘高 ','FontSize',20);
? ? text(ax,0.4*L,0.4*L,1.75*h,...
? ? ? ? ['$V=',num2str(V,'%.1f'),'$'],...
? ? ? ? 'Interpreter','latex','FontSize',20);
? ? ? ??
%勢能偽色圖
? ? imagesc(ax,[-0.5*L,0.5*L],[-0.5*L,0.5*L],...
? ? ? ? h*Ub,'AlphaData',0.1);
? ??
%顯示相位，閾值是0.0001*h，可按需調(diào)整
? ? ? ? phase1 = ...
? ? ? ? h*angle(psi21).*heaviside(rho1-0.0001*h);
? ? ? ? phase2 = ...
? ? ? ? h*angle(psi22).*heaviside(rho2-0.0001*h);
? ? Ph1 = surf(ax,x,x,-Z0,phase1,'FaceAlpha',0.3);
? ? shading interp
? ? Ph1.EdgeColor = 'none';
? ??
? ? Ph2 = surf(ax,x,x,Z0,phase2,'FaceAlpha',0.3);
? ? shading interp
? ? Ph2.EdgeColor = 'none';
? ? hold off
? ? % 顏色軸和三維視圖設(shè)置
? ? colormap(hsv)
? ? caxis([-pi*h pi*h])
? ? pbaspect([1 1 0.5]);
? ? campos([-15 -32 2.2*h])
? ? camtarget([0 0 0.3*h])
? ? camva(20)
%立刻作圖
? ? drawnow
end
