clear all; close all; clc;

L = 50; W = 50; ff = 0.1; % L,W are system dimensions. For archaic reasons (that you don't need to worry about), they are in units of HALF \ell, the unit of distance used in the paper. So set them to TWICE the value you want in units of \ell. ff is the fill factor, called 'f' in the paper. The definition of this parameter is the same as in the paper, so it is easy to set.

load blah.mat; % you should run ValleyNetwork.m first for the noise realisation of interest (which you should save to use in DiagH.m, too). Save the cell array 'info' as blah.mat, and load it here. It contains 'u' (the localisation landscape), from which we get W_E. Note that due to laziness on my part, this code will only work properly for 1x1 partitions in ValleyNetwork.m (i.e. Px=Py=1, see that code).

box = 2; % do not change! converts from archaic units of length (\ell/2) to paper-consistent (\ell).
Lf = L/box; Wf = W/box; % these are now in \ell.
Nx = L; % choose number of Cheb points in x
Ny = W; % choose number of Cheb points in y
num = round(ff*W*L/4); % total number of scatterers. The 4 here is hard-coded (don't change) due to the use of L and W and the unit thing above.

g = 21.33; % scatterer height in units of E_0, as defined in the paper.

many = 30; % how many e-vals and e-vecs you want
where = 'sm'; % where to look on the enegry scale ('sm' means smallest, keep it that way)
how = 'bo'; % how to plot energy spectrum
BCtype = 'Dirichlet'; % 'Dirichlet' & 'Neumann' are allowed options for x.
getevecs = 1; % do you want e-vecs as well as e-vals? (1=yes, 0=no)
fulldiag = 0; % would you like to use full matrices rather than sparse? (not recommended. 1=yes, 0=no)

% construct x (longitudinal dim) laplacian
Lx = chebop(@(u) -diff(u,2),[0,Lf],'dirichlet');
[Lx,IDx] = MyDisc(Lx,Nx); % uses 'MyDisc', which I got from chebfun and slightly modified for my purposes.
Lx(1:2,:) = []; IDx(1:2,:) = [];
% construct y (transverse dim) laplacian
Ly = chebop(@(u) -diff(u,2),[0,Wf],'dirichlet');
[Ly,IDy] = MyDisc(Ly,Ny);
Ly(1:2,:) = []; IDy(1:2,:) = [];

% tensor product (note order!)
Lap = kron(Ly,IDx) + kron(IDy,Lx);
ID = kron(IDy,IDx);

% create Cheb discretized vecs
nx = (0:Nx+1);
x = cos(pi*nx/(Nx+1));
x = x + 1;
x = x*Lf/2;
ny = (0:Ny+1);
y = cos(pi*ny/(Ny+1));
y = y + 1;
y = y*Wf/2;

if ff>0
    [X,Y] = meshgrid(x,y);
    
    u=info{1,2}; % so this is where we pull up 'u' from LLT and set the potential = W_E
    Scat = 1./feval(u,X,Y);
       
    Scat = sparse(Scat);

    blah = Scat.';
    blah = blah(:);
    ChebScat = diag(blah);
    ChebScat = sparse(ChebScat);
else
    ChebScat = 0;
end

DownSampled = sparse(ID*ChebScat);

top = sparse(zeros(Nx+2,(Nx+2)*(Ny+2)));
for j=1:Nx+2
    top(j,j) = 1;
end

bot = sparse(zeros(Nx+2,(Nx+2)*(Ny+2)));
for j=1:Nx+2
    bot(j,(Nx+2)*(Ny+1)+j) = 1;
end

if strcmp(BCtype,'Dirichlet')
    
left = sparse(zeros(Ny,(Nx+2)*(Ny+2)));
for j=1:Ny
    left(j,j*(Nx+2)+1) = 1;
end

right = sparse(zeros(Ny,(Nx+2)*(Ny+2)));
for j=2:Ny+1
    right(j-1,j*(Nx+2)) = 1;
end

elseif strcmp(BCtype,'Neumann')
    
Dx = zeros(Nx+2,Nx+2);
zilch = zeros(Nx+2,1);
for n=1:Nx+2
uni = zilch; uni(n) = 1;
Dx(:,n) = fchd(uni).';
end
Dx = -2/Lf*Dx;
Iy = speye(Ny+2);
NBC = sparse(kron(Iy,Dx));

left = sparse(zeros(Ny,(Nx+2)*(Ny+2)));
for j=1:Ny
    left(j,:) = NBC(j*(Nx+2)+1,:);
end

right = sparse(zeros(Ny,(Nx+2)*(Ny+2)));
for j=2:Ny+1
    right(j-1,:) = NBC(j*(Nx+2),:);
end
end

BCs = sparse([top; left; right; bot]);
empty = sparse(zeros(size(BCs)));

DownSampled = [empty; DownSampled];

Lap = [BCs; Lap];
ID = [empty; ID];

H = Lap + DownSampled; % Hamiltonian in units of E_0

if getevecs
    [V,E] = eigs(H,ID,many,where);
elseif ~getevecs && ~fulldiag
    [~,E] = eigs(H,ID,many,where);
elseif ~getevecs && fulldiag
    [~,E] = eig(full(H),full(ID));
end

[E,jones] = sort(diag(E));
E = E(1:many);
indy = (1:1:many);
plot(indy,E,how); grid on; hold on; % plots the eigenvalues in increasing order. Units are E_0.


