clear all; close all; clc;

numruns = 20; % Set number of realisation to average over.
meany = zeros(101,1351); % Number of rows = length(time) from tdepsims.m. Number of columns = length(x) from tdepsims. You need to look these up by loading one realisation and checking how long these vectors are before you run this.

% The different realisations to be averaged are assumed to be in the working directory. You may need to copy this script into the subfolder containing the data and run it there.
for count=1:numruns
    load(strcat(int2str(count),'.mat'));
    meany = meany + den1d;    
end

meany = meany/numruns;

% Now we plot movies for you to look at (hit any key to move from frame to frame), first on a linear scale, then on a log scale. Red vertical lines show the boundaries between the reservoirs and the channel.

for k=1:length(time)
    k
    plot(x,meany(k,:),'bo'); grid on; hold on;
    plot([0,0],[min(meany(k,:)) max(meany(k,:))],'r-');
    plot([Lf,Lf],[min(meany(k,:)) max(meany(k,:))],'r-');
    hold off; pause;
end

for k=1:length(time)
    k
    plot(x(10:end-10),log(meany(k,10:end-10)),'bo'); grid on; hold on;
    plot([0,0],[min(log(meany(k,10:end-10))) max(log(meany(k,10:end-10)))],'r-');
    plot([Lf,Lf],[min(log(meany(k,10:end-10))) max(log(meany(k,10:end-10)))],'r-');
    hold off; pause;
end

