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

L = 15; % Interval Length

N = 200; % 設(shè)置 N*N 格點(diǎn)

% x = linspace(-L/2,L/2-L/(N),N)'; % Coordinate vector， 周期性邊界時(shí)

x = linspace(-L/2,L/2,N)'; % Coordinate vector，硬邊界時(shí)

dx = x(2) - x(1); % Coordinate step

% operators in 1D:

e = ones(N,1);?

% Dif1 = spdiags([-e 0*e e -e e],[-1 0 1 N-1 -N+1],N,N)/(2*dx);% 有周期性邊界

% Lap1 = spdiags([e -2*e e e e],[-1 0 1 N-1 -N+1],N,N)/dx^2;

Dif1 = spdiags([-e 0*e e],[-1 0 1],N,N)/(2*dx);% 硬邊界

Lap1 = spdiags([e -2*e e],[-1 0 1],N,N)/dx^2;?

I1 = speye(N,N);

% operators in 2D:

Lap2 = kron(Lap1,I1)+kron(I1,Lap1); % 2D Laplacian

PDx = kron(Dif1,I1); % partial diff of x

PDy = kron(I1,Dif1);% partial diff of y

xy = spdiags(x,0,N,N);

X = kron(xy,I1); % x position operator

Y = kron(I1,xy); % y position operator

R2 = kron(xy.^2,I1)+kron(I1,xy.^2); % R^2=x^2+y^2

L3 = -1i*(X*PDy-Y*PDx);% angular momentum, z向

B = 4; % 磁感應(yīng)強(qiáng)度

H = 0.5*(-Lap2+B*L3+0.25*B*B*R2); % 哈密頓算子

NT = 50000; % No. of time steps

TF = 50; T = linspace(0,TF,NT); % Time vector

dT = T(2)-T(1); % Time step

% E = expm(-1i*H*dT); % Evolution

I = speye(size(H));

A = -1i*H*dT;

E = I;

for m=1:11

? ? E = E+(A^m/factorial(m));

end

% spy(E)

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

e1 = ones(N,1);

xf = kron(x,e1);

yf = kron(e1,x);

NLL = 4;% Landau 能級量子數(shù)，與z角動量量子數(shù)

a = 3;% 控制波包大小

d = 0.5*L;% 平移距離，平移波函數(shù)時(shí)使用

% 以原點(diǎn)為中心的旋轉(zhuǎn)對稱的朗道能級波函數(shù)

psi0 = ((xf+1i*yf).^NLL).*exp(-0.25*B*(xf.^2+yf.^2));?

% 中心按照 d 在 x 方向平移的朗道能級波函數(shù)

% psi0 = exp(-1i*0.5*B*yf*d).*(((xf-d)+1i*yf).^NLL).*exp(-0.25*B*((xf-d).^2+yf.^2));

% psi0 = psi0.*exp(-((xf-0.3*L).^2+(yf).^2)/a); % 用于截取波包

psi0 = psi0/sqrt(psi0'*psi0); %歸一化

EM=psi0'*H*psi0 % 能量期望

LM=psi0'*L3*psi0; % Lz 期望

% 理論預(yù)測的能量本征值是 B*(2*NLL+1)/2

ES=sqrt(real(psi0'*H*H*psi0-EM.^2)) % 能量標(biāo)準(zhǔn)差

LS=sqrt(real(psi0'*L3*L3*psi0-LM.^2)); % Lz 標(biāo)準(zhǔn)差

% 計(jì)算密度

rho = abs(psi0).^2;?

rho1 = reshape(rho,[N,N]); % 變成二維數(shù)組，方便畫圖

% 在密度足夠大的范圍內(nèi)計(jì)算相位，閾值設(shè)為0.05*h

h=1e-3;

?phase = reshape(angle(psi0).*heaviside(abs(psi0).^2-0.05*h),[N,N]);