clear all; close all; clc;

%global Rf Lf Wf g % If you want to run an empty channel, activate this instead of the next line
global posx posy g num Rf Lf Wf % Use with scatterers in the channel

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.

Lf = L/2; Wf = W/2; % These are now in \ell.
R = 60; Rf = R/2; % Length (along x) of reservoirs on either side of the disordered channel.

% Generate noise realisation
num = round(ff*W*L/4);
posx = randi(Lf,[1,num]);
posy = randi(Wf,[num,1]);
%load blah.mat % OR load a pre-saved one if you like

g = 5; % Scatterer height in units of E_0, as defined in the paper.

model = createpde(1);

shape = [3;4;-Rf;Lf+Rf;Lf+Rf;-Rf;0;0;Wf;Wf];
geom = decsg(shape);
geometryFromEdges(model,geom);
applyBoundaryCondition(model,'dirichlet','Edge',[1,2,3,4],'u',0);

% If you want to run an empty channel, use this line:
    %specifyCoefficients(model,'m',0,'d',-1i,'c',1,'a',0,'f',0);
% If you want to evolve in V, use this line:
    specifyCoefficients(model,'m',0,'d',-1i,'c',1,'a',@ScatPot,'f',0); % Go into 'ScatPot.m' and set width of scatterers!!!

generateMesh(model,'Hmax',1); % This is an important resolution parameter. For all simulations in the paper, OK to use 1. If you go to high energies (by increasing k_0), you may wish to reduce Hmax to 1/2 and compare to a run with 1 (obviously use the same noise realisation for this comparison!). Large systems with Hmax=1/2 will take forever or will max out your RAM (just a friendly warning).
setInitialConditions(model,@Gaussian); % Open 'Gaussian.m' and set parameters for the initial condition!!!

time = (0:1:20); % Time vector in t_0. It should start at zero. Feel free to change the step size or final evolution time as you need to. Note that solver uses its own step size for solving -- this is purely for evaluating and plotting.
result = solvepde(model,time);

step = 0.1; % Again, this step size here is for evaluation purposes only. The solver does its own thing.
x = (-Rf:step:Lf+Rf);
y = (0:step:Wf);

[X,Y] = meshgrid(x,y);

den1d = zeros(length(time),length(x)); % We'll be wanting to calculate and store the longitudinal density.

for k=1:length(time)
    psie = interpolateSolution(result,[X(:),Y(:)]',k);
    psie = reshape(psie,size(X));
    
	% Enable this line to see a 2D movie, which you have to flick through by hitting a key every time in pauses.
    %surface(X,Y,abs(psie).^2); shading interp; colormap('jet'); grid on; pause; close all;
        
    if k==1
        total = trapz(x,trapz(y,abs(psie).^2)); % The initial condition is not normalised, so we do that here.
    end
    
    den1d(k,:) = trapz(y,abs(psie).^2/total);
end

% Now plot a 1D movie, which you have to flick through by hitting a key every time in pauses.
for k=1:length(time)
    k
    plot(x,den1d(k,:),'bo'); grid on; pause;
end



