% Visualise the drift velocity of the particle

clear all

Data = dlmread('paperfigure8a');
%Data = Data';
binWidth = Data(1,1); % in microns
rows = Data(1,2);
columns = Data(1,3);

% Unpack data
start = 2;

Occ1 = Data(start:start+rows-1,:);
vMag1 = Data(start+rows:start+2*rows-1,:);
vX1 = Data(start+2*rows:start+3*rows-1,:);
vY1 = Data(start+3*rows:start+4*rows-1,:);

Occ2 = Data(start+4*rows:start+5*rows-1,:);
vMag2 = Data(start+5*rows:start+6*rows-1,:);
vX2 = Data(start+6*rows:start+7*rows-1,:);
vY2 = Data(start+7*rows:start+8*rows-1,:);


maxOcc = max(max(Occ2));
maxSpeed = max(max(vMag2));

% maxOcc2 = max(max(Occ2));
% for n = 1:rows
%     for m = 1:columns
%         if Occ2(n,m)>0
%             Occ2(n,m) = maxOcc2 - Occ2(n,m);
%         end
%     end
% end

%plot magnitude data

X = 0:binWidth:((rows-1)*binWidth);
Y = 0:binWidth:((columns-1)*binWidth);
Z = zeros(columns,rows);
[X, Y] = meshgrid(X,Y);

% Create figure
figure0 = figure('Color',[1 1 1]);
colormap(jet(512));
%gray(512)]
% Create axes
axes0 = axes('Parent',figure0, 'DataAspectRatio',[1 1 1]);
box(axes0,'on');
hold(axes0,'all');

% Create surface
h(1) = surface('Parent',axes0,'ZData',Z,'YData',Y,'XData',X,...
    'LineStyle','none',...
    'FaceColor','flat',...
    'CData',Occ2');
axis tight

% Create figure
figure1 = figure('Color',[1 1 1]);
colormap(jet(512));
%gray(512)]
% Create axes
axes1 = axes('Parent',figure1, 'DataAspectRatio',[1 1 1]);
box(axes1,'on');
hold(axes1,'all');

% Create surface
h(1) = surface('Parent',axes1,'ZData',Z,'YData',Y,'XData',X,...
    'LineStyle','none',...
    'FaceColor','flat',...
    'CData',vMag2');
axis tight

reductionFactor =3;

rowsq = floor(rows/reductionFactor);
columnsq = floor(columns/reductionFactor);
Xvecq = zeros(rowsq,columnsq);
Yvecq = zeros(rowsq,columnsq);
Xq = zeros(columnsq,rowsq);
Yq = zeros(columnsq,rowsq);

s=floor(reductionFactor/2);
for n=1:rowsq
    for m=1:columnsq
        Xvecq(n,m) = vX2(n*reductionFactor-s,m*reductionFactor-s);
        Yvecq(n,m) = vY2(n*reductionFactor-s,m*reductionFactor-s);
        Xq(m,n) = X(m*reductionFactor-s,n*reductionFactor-s);
        Yq(m,n) = Y(m*reductionFactor-s,n*reductionFactor-s);
    end
end

quiver(Xq,Yq,Xvecq',Yvecq','Color',[0 0 0],'AutoScaleFactor',1.3);