% This first half (till line 76) should be done once per noise realisation and then switched off. After that you can run the second half for several e-states in this same realisation.
if 1
% Load desired workspace and save noise realisation in working folder
clear all; close all; clc;
L = 50; W = 50; ff = 0.06; % 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.
count = 1; % Noise realisation label
folder = strcat('L',int2str(L),'W',int2str(W),'ff',num2str(ff),'P');
folder = strrep(folder,'.','_');
cd SimCity; % pick folder
cd(folder);
load(strcat('wrkspc',int2str(count),'.mat'));
cd ../..

% To use H with V, include this line:
save('blah.mat','posx','posy');
% To use H with W_E, include this line INSTEAD:
%save('blah.mat','info');

% Now run exact diagonalisation and save e-vals and e-vecs.
% To use H with V, include this line:
DiagH; % Open this script first, and comment out lines 1 & 3. Make sure that scatterer positions (posx and posy) are loaded, not generated anew.
% To use H with W_E, include this line INSTEAD:
%DiagHwithWE; % Open this script first, and comment out lines 1 & 3.

close all;
evecs = V;
evals = E;
save('exactresults.mat','evecs','evals','jones','Nx','Ny','many');

% Reload desired workspace
clear all; close all; clc;
L = 50; W = 50; ff = 0.06; % Because we've just cleared the workspace, you have to set these params to match the above load section
count = 1;
folder = strcat('L',int2str(L),'W',int2str(W),'ff',num2str(ff),'P');
folder = strrep(folder,'.','_');
cd SimCity; % pick folder
cd(folder);
load(strcat('wrkspc',int2str(count),'.mat'));
cd ../..
load exactresults.mat

% Define x & y for plotting
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;

% Plot e-vecs and note which e-vecs and which domains to examine (I've provided my selection in the HowTo file...)
for n=1:many
evec = reshape(evecs(:,jones(n)),[Nx+2,Ny+2]).';
figure;
surface(x,y,log(abs(evec))); colormap('jet'); shading interp; grid on; hold on;

for k=1:length(C)
    traj = C{k,1};
    c = find(traj(3,:)<1/evals(n));
    plot3(traj(1,c),traj(2,c),traj(3,c),'k.');
    traj = C{k,2};
    c = find(traj(3,:)<1/evals(n));
    plot3(traj(1,c),traj(2,c),traj(3,c),'w.');
end

for k=1:size(maxima,1)
    plot(maxima(k,1),maxima(k,2),'ko');
    text(maxima(k,1),maxima(k,2),strcat('\leftarrow ',int2str(k)),'FontSize',20);
end

pause; close all;
end

end

if 1
eignum = 1; % Pick eigenstate number. Turn off the above section after running it once, turn this on, and then you can run it several times for diff eigenstates without clearning workspace or redoing the above.

evec = reshape(evecs(:,jones(eignum)),[Nx+2,Ny+2]).';
approx = chebfun2(fliplr(flipud(abs(evec))),[0 Lf 0 Wf]);
evals(eignum) % Prints out eigenvalue of the eigenvector you picked

% integrate abs(psi) over the domains
psi_int = zeros(length(maxima),1);

for n=1:length(maxima)
    [r,c] = find(owners==n);
    walls = [];
    for k=1:length(r)
        traj = C{r(k),1};
        walls = [walls, traj];
        traj = C{r(k),2};
        walls = [walls, traj];
    end
    walls(3,:)=[];
    r = find(owners2==n);
    for k=1:length(r)
        traj = ChannelWalls{r(k)};
        walls = [walls, traj];
    end
        
    box = [round(min(walls(1,:))/e)*e round(max(walls(1,:))/e)*e round(min(walls(2,:))/e)*e round(max(walls(2,:))/e)*e];
    x = (box(1):e:box(2));
    y = (box(3):e:box(4));
    Psimat = zeros(length(y),length(x));
    
    walls = round(walls/e)*e;
    IntBnds = zeros(length(y),2);
    for k=1:length(y)
        r = find(round(walls(2,:)/e)==round(y(k)/e));
        if ~isempty(r)
            blah = min(walls(1,r));
            IntBnds(k,1) = find(round(x/e)==round(blah/e));
            blah = max(walls(1,r));
            IntBnds(k,2) = find(round(x/e)==round(blah/e));
        end
    end
    
    [r,c] = find(IntBnds==0);
    for k=1:length(r)
        IntBnds(r(k),c(k)) = (IntBnds(r(k)-1,c(k)));
    end
    
    xk = (0:e:Lf);
    yk = (0:e:Wf);
    xc = intersect(round(xk/e),round(x/e))*e;
    yc = intersect(round(yk/e),round(y/e))*e;
        
    [Xc,Yc] = meshgrid(xc,yc);
    approxe = feval(approx,Xc,Yc);
    Psimat(round((yc-box(3))/e+1),round((xc-box(1))/e+1)) = approxe;
    blah = zeros(size(Psimat));
    for k=1:length(y)
        blah(k,IntBnds(k,1):IntBnds(k,2)) = Psimat(k,IntBnds(k,1):IntBnds(k,2));
    end
    Psimat = blah;
    clear blah
        
    psi_int(n) = trapz(x,trapz(y,Psimat,1));
end

weight = psi_int./area; % Calculate average of |\psi| on each domain

% Plot for visual inspection by a human
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;

surface(x,y,log(abs(evec))); colormap('jet'); shading interp; grid on; hold on;

for k=1:length(C)
    traj = C{k,1};
    c = find(traj(3,:)<1/evals(eignum));
    plot3(traj(1,c),traj(2,c),traj(3,c),'k.');
    traj = C{k,2};
    c = find(traj(3,:)<1/evals(eignum));
    plot3(traj(1,c),traj(2,c),traj(3,c),'w.');
end

plot(maxima(:,1),maxima(:,2),'bo','MarkerSize',6,'MarkerFace','b');
end
