Homec4science
Diffusion Computational physics (master) Browse

compphys/exo3/landau.m
33bba04159d6master

Last Change
Raw File

landau.m
View Options

global a
global B
global hop
global N
hop = -1; % hopping integral
N = 50;
B = 0.05;
L = 1.0;
a = L / (N+1);
%bc = "periodic"
bc = "dirichlet"
if abs(B) > eps %&& bc == "periodic"
a = sqrt(4 * pi / (B * N));
end
% print 'a' parameter
a
% size(H) = [N^2, N^2]
H = zeros(N^2);
% bind values to the hamiltonian
for i = 1:N
for j = 1:N
% debug
%fprintf("Index: %d", grid_k(i,j))
% diagonal
H(grid_k(i,j), grid_k(i,j)) = 4;
% off diagonal
if i > 1
H(grid_k(i,j), grid_k(i-1,j)) = hop_int_D(j);
end
if j > 1
H(grid_k(i,j), grid_k(i,j-1)) = hop_int_L(i);
end
if j < N
H(grid_k(i,j), grid_k(i,j+1)) = hop_int_R(i);
end
if i < N
H(grid_k(i,j), grid_k(i+1,j)) = hop_int_U(j);
end
end
end
if bc == "dirichlet"
% nothing to do
elseif bc == "periodic"
% for each x distretization (j): find y boundary
for j = 1:N
H(grid_k(1,j), grid_k(N,j)) = hop_int_D(j); % i = 1
H(grid_k(N,j), grid_k(1,j)) = hop_int_U(j); % i = N
end
% for each y distretization (i): find x boundary
for i = 1:N
H(grid_k(i,1), grid_k(i,N)) = hop_int_L(i); % j = 1
H(grid_k(i,N), grid_k(i,1)) = hop_int_R(i); % j = N
end
end
% normalize
H = H / a^2;
% H debugging
%diag(H)'
%diag(H,-1)'
%diag(H,1)'
% solving
[psi, E] = eig(H);
E = diag(E)
%energies = E'
[~, index] = min(E);
%index = 2;
E_0 = E(index)
psi_0 = abs(psi(:,index)).^2;
adds_on = [];
%adds_on = [adds_on 135, 239, 341] % stable
%adds_on = [adds_on 97, 194, 296, 394, 494] % 1/2 transition
adds_on = [adds_on 85, 187, 282, 679, 2097, 2203] % 1/4 transition
%adds_on = [adds_on 1250] % central point
%adds_on = [76, 171, 268] % bottom discontinuity
%adds_on = [104, 208, 312] % top discontinuity
%eig_v = psi_0'
% output data to a file
out = fopen(sprintf('output/landau_%s_B%.3f_N%d.out', bc, B, N),'w');
for k = 1:(N^2)
fprintf(out, '%d %g\n', k, E(k));
end
fclose(out);
% plot lowest level
figure
xlabel('x [m]')
ylabel('y [m]')
pcolor(x(1:N), y(1:N), reshape(psi_0, N, N))
%contourf(x(1:N), y(1:N), reshape(psi_0, N, N))
for i = 1:length(adds_on)
fprintf('Printing index %d\n', adds_on(i))
figure
xlabel('x [m]')
ylabel('y [m]')
pcolor(x(1:N), y(1:N), reshape(abs(psi(:,adds_on(i))).^2, N, N))
end
% plot eigenvalues
figure
xlabel('Index')
ylabel('Eigenvalue')
plot(1:(N^2), E)
% functions
function k = grid_k(i, j)
global N
k = (i-1) * N + j;
end
function pos = x(j)
global a
pos = (j-1) * a;
end
function pos = y(i)
global a
pos = (i-1) * a;
end
function g = hop_int_R(i)
global hop
global B
global a
g = hop * exp(-1i * B * y(i) * a / 2);
end
function g = hop_int_L(i)
global hop
global B
global a
g = hop * exp(1i * B * y(i) * a / 2);
end
function g = hop_int_U(j)
global hop
global B
global a
g = hop * exp(1i * B * x(j) * a / 2);
end
function g = hop_int_D(j)
global hop
global B
global a
g = hop * exp(-1i * B * x(j) * a / 2);
end
c4science · Help