%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALCULATE THE ENERGY LEVELS OF A 2D QUANTUM RING WITH FINITE WALLS 	    %
%                                                                           %
% The Schroeding equation is solved by means of a finite differences scheme %
%               J. Planelles and J.I. Climente, September 2011              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

% INPUT DATA

Rin_nm=3.0;  	% inner radius (3 nm) - if zero, we have a disk instead of a ring
W_nm=30.0;      % ring section (nm), Rout=Rin+W
%Mz=1;		    % azimuthal angular momentum
Bmin=0.0;       % smallest magnetic field (T)
Bmax=20.0;      % largest magnetic field (T)
Bstep=1.0;      % increment of B (nm)
Lb_nm=10.0;     % length of the barrier (nm)
m=0.05;         % electron effective mass (times m0)
Vin_eV=0.00;    % potential energy inside the QD (eV)
Vout_eV=0.25;   % potential energy outside the QD (eV)
h_nm=0.25;      % grid step (nm)
neig=3;         % number of states to be calculated (starting from ground state)

% CALCULATION
for Mz=-20:1:8

for B_T=Bmin:Bstep:Bmax % open loop of B

% 1) Convert into atomic units
kdist=0.0529177; % kdist nm = 1 bohr radius
kener=27.21138;  % kener eV = 1 hartree
tesla=235054;	 % tesla T  = 1 at.unit B

Rin=Rin_nm/kdist;
W=W_nm/kdist;
Lb=Lb_nm/kdist;
h=h_nm/kdist;
Vin=Vin_eV/kener;
Vout=Vout_eV/kener;
B=B_T/tesla;

% 2) Define grid

Rout=Rin+W;        % outer radius
Lt= Rout + Lb;     % total length
n=round(Lt/h) + 1; % number of grid points

% 3) Calculate potential V(x)

for i=1:n
    x(i)= (i-1)*h; % position at i-th grid point
    
    V(i)=Vout;

    if (x(i) > Rin -1e-12) && (x(i) < Rout+1e-12)
        V(i)=Vin;  % inside the QD
    end;

end;
    
% 4) Define coefficients of f(i), f(i+/-1)

for i=2:n
    c(i)=-1/(2*m)*(-2/h^2 - Mz^2/x(i)^2 - B*Mz) + (B*x(i))^2/(8*m) + V(i); 	% coef of f(i)
    a(i)=-1/(2*m)*(1/h^2 - 1/(x(i)*2*h));					% coef of f(i+1)
    b(i)=-1/(2*m)*(1/h^2 + 1/(x(i)*2*h));					% coef of f(i-1)
end;

% if Mz=0, amend BC 
if (Mz==0)
	c(2)=c(2)+a(2);
end


% 5) Build the Hamiltonian matrix

% define the three diagonals (vectors)
dc=c(2:n-1)';   % main diagonal f(i)
di=a(3:n)';     % subdiagonal   f(i-1)
ds=b(1:n-2)';   % superdiagonal f(i+1)

% introduce them into a sparse, tridiagonal matrix
matrix=spdiags([di,dc,ds],[-1,0,1],n-2,n-2);

% 6) Diagonalize the matrix
[P,D]=eigs(matrix,neig,'sm');
eners=sort(diag(D))*kener; % grep energy -diagonal of D-, sort it and convert it to eV

% WRITE RESULTS TO DISK

% Write energy to external file
file=fopen('energies_B.dat','a'); % open file to append results

fprintf(file,'%10.1f',B*tesla);   % write QD size (nm)
for i=1:neig
    fprintf(file,'%10.3f',eners(i)*1000); % write eigenvalues (meV)
end;
fprintf(file,'\n');

fclose(file);

end; % close loop of B
end; % close loop of Mz
