function ROIPostProcess(TR,ED,R0_temp,L,tau_in,TM,NT,PI,...
    BROI,BT0,BTF,PROI,PT0,PTF,DROI,DT0,DTF,TN,ProtNo,PhiMov,ROIRetractMov,...
    ColorSurfTraffMov)
% clf

global Trial ExpDuration TotalMass TotalFree PlotTime...
    NodeIndex RaftNodes SurfNodes WaterNodes NodeIndexPrev...
    SurfNodesPrev RaftNodesPrev WaterNodesPrev RaftColor...
    SurfaceColor RedExp Outline TenMil LengthScale ...
    NominalSpringStiffness TrajectoryMu TrajectorySigma  ...
    BackgroundLattice InitialDens SurfFractionThresh nt dt it...
    AvgAntLength NoPauseIt GridSize Nrow Ncol tau_dec ct...
    DomainDensThresh RaftDensThresh NWindow NewNeighDisThresh...
    ParticleSize TotalTime mmPerNode zeta InputInterval...
    RadNeighbors TrajectoryNoise TrajectoryType CrowdingDistFactor...
    ChooseBasedOnDisp MaxDOF PressureRadius parkrate unparkrate...
    SurfNodesHistory Neighbors Nr Ns Ntot gamma alpha delta phi...
    SurfExcess tau Area Circumference WorkLengthScale R0 BulkROI...
    Bulkt0 Bulktf ProtROI Prott0 Prottf DecayROI Decayt0 Decaytf

%% I. Inputs
%% Tracking Inputs
nt = NT;            %No. Successful iterations
xSpeed = 100;           %speedup factor for video
Trial = TR;
R0 = R0_temp;       %Vicsek radius
WorkLengthScale = L;       %free parking curvature
ExpDuration = ED;
InputInterval = PI;
PlotInterval= PI;
FontSize = 20;
BulkROI = BROI;
Bulkt0 = BT0;
Bulktf = BTF;
ProtROI = PROI;
Prott0 = PT0;
Prottf = PTF;
DecayROI = DROI;
Decayt0 = DT0;
Decaytf = DTF;

%% Clear Previous Values
NodeIndexPrev = []; SurfNodesPrev =[]; RaftNodesPrev = []; WaterNodesPrev = [];

%% Unit Conversion and Ant Property Inputs
AntMass = 1e-3; %g/ant
TotalMass = TM;  %g   Based on Experiment 10 - timelapse
TotalFree = floor(TotalMass/AntMass);%3*floor((1+SurfDens/RaftDens)*pi*Rad^2);% %Total number of free ants (not used for SourceLoc ==0,1)
    %if SourceLoc==3, once N_add>TotalFree, surface density will decay
    %until rho_t==SurfFractionThresh, then partilces will be added to middle to
    %maintain rho_t==SurfFractionThresh
AvgAntLength = 2.93;   %mm/body length (1 zeta = 1 ant);
RaftDens = 0.303747695;  %ants/mm^2
RaftDensNorm = 1;   %normalized number of nodes per square node area
mmPerNode = sqrt(1/RaftDens);  %mm/node dimension
SurfDens = 0.072242475; %ants/mm^2
zeta = 1/sqrt(RaftDens);    %nodes per ant body length - normalization parameter & assumed reach of ants in nodes
AntSpeed_zeta = 0.52;    %zeta/s
AntSpeed_mm = AntSpeed_zeta*AvgAntLength;
AntSpeed_node = AntSpeed_mm/zeta;   %ant speed in nodes/s
TenMil = 10/mmPerNode;   %zeta
dt = 1.7/AntSpeed_node;    %s/iteration assuming particles move one body length per iteration 

%% Domain Inputs
NWindow = 10*ceil((0.5*floor(TM/2.75/AntMass)/pi)^(1/2));%10*ceil((0.5*floor(2.2/AntMass)/pi)^(1/2))
GridSize = 2;    %size of a square portion of grid for density-based unparking in nodes-sizes
Nrow = NWindow+2*GridSize;%round(NWindow/(exp(-ExpDuration/tau_dec)),-1); %Size of lattice domain
Ncol = Nrow;
DomainDensThresh = 0; %local density below which particles will be added to edge of domain in nodes/ node length^2
    %Note: Introduced to smooth out grain boundary. DO NOT CHANGE
MaxDOF = 18;   %DOF of movement and connectivity - dist. based after t=0

%% Surface Particle Inputs
InitialDens = 1;%SurfDens/RaftDens; %fraction of surface nodes occupied by active ants (maintained)
SurfFractionThresh = SurfDens/RaftDens; %threshold that surface density will be held at if SourceLoc==3
PauseDuration = 5; %how long ants pause in place in seconds
NoPauseIt = ceil(PauseDuration/dt); %number of iterations particles will park for
TrajectoryType = 1; %Determines mode of movement determination
    %0 for semi-ballistic trajectories with Guassian distribution of
        %turning probability
    %1 for numerical implementation of Vicsek model - ants parking on edge
        %also follow Vicsek model (will not park unless it is compatible 
        %with  their Vicsek-predicted trajectories)
TrajectoryMu = 0;   %mean of Gaussian for directional change
TrajectorySigma = pi/8; %std dev. of Gaussian for directional change
%Inputs if TrajectoryType==1 (Viscek)
RadNeighbors = 1.5;   %Radius of influence for surface particles in Vicsek model
PressureRadius = 1.6;%RadNeighbors;   %Radius over which pressure can be felt on particles
TrajectoryNoise = TN;  %Noise factor for the directional change of a particle
    %in range of 0 (no noise) to 1 (360 degrees of noise)
    %Previous analysis showed that the code is extremely sensitive to this.
CrowdingDistFactor = 1;
ChooseBasedOnDisp = 0;

%% Raft Inputs
NominalSpringStiffness = 0.2;   %Nominal spring stiffness... WHAT ARE UNITS?
    %Note: pick a number that Ncol and Nrow are both divisible by
    %Note: Pick an even number
    %Note: setting smaller will increase computational demand but introduce
        %higher precision over which nodes to remove. Set to 1 for final
        %analysis
NewNeighDisThresh = 2.5; %Max dist. across which new neighbors can form
tau_dec = tau_in; %linear decay rate of raft based on exp. data in s^-1
RaftDensThresh = round(GridSize^2*RaftDensNorm)+1; %local density above which raft particles will unpark in nodes/node length^2
    %NumNodes/GridArea at which an unpark event will occur to maintain
    %nominal raft denstity
    
%% Visualization Options
BackgroundLattice = 0;
RaftColor = 'c';    %raft particle color
SurfaceColor = 'r'; %active particle color
RedExp = 1.15;   %Reduces Marker size by 1/Nrow^RedExp
Outline = 1;    %turn particle outline on
LengthScale = 1;    %1 to turn on a 10mm length bar
ParticleSize = 0.75;

%% Prep Results Output
mkdir('Charts')
mkdir('Movies')
mkdir('Output Data')

if PhiMov==1
    MovieName1 = ['Movies/ROIRaft.T',num2str(Trial)];
    v1 = VideoWriter(MovieName1);
    v1.FrameRate = round(xSpeed/dt);
    open(v1)
end

if ROIRetractMov==1
    MovieName2 = ['Movies/ROIStrainRate.T',num2str(Trial)];
    v2 = VideoWriter(MovieName2);
    v2.FrameRate = round(xSpeed/dt);
    open(v2)
end

if ColorSurfTraffMov==1
    MovieName3 = ['Movies/ROISurfTraffic.T',num2str(Trial)];
    v3 = VideoWriter(MovieName3);
    v3.FrameRate = round(xSpeed/dt);
    open(v3)
end

figure(1)
clf
% figure(2)
% clf

parkrate = [];
unparkrate = [];
PlotTime = [];
gamma = [];
delta = [];
alpha = [];
Nr = [];
Ns = [];
Ntot = [];
Area = [];
Circumference = [];
phi = [];
SurfExcess = []; 
tau = [];
it = [];
TotalTime = [];

PlotCt = 0; %counts number of plots
ct = 0; %counts number of iterations
for it=InputInterval:InputInterval:nt
    ct = ct+1;
    PercentComplete = it/nt*100
    if PercentComplete==100
        break
    end
    
    %% Import Data
    NI_FileName = ['Data/NodeIndex.Trial',num2str(Trial),'.it',num2str(it),'.txt'];
    SN_FileName = ['Data/SurfNodes.Trial',num2str(Trial),'.it',num2str(it),'.txt'];
    Neighb_FileName = ['Data/Neighbors.Trial',num2str(Trial),'.it',num2str(it),'.txt'];
    
    NodeIndex = dlmread(NI_FileName);
    SurfNodes = NodeIndex(NodeIndex(:,3)==1,:);
    RaftNodes = NodeIndex(NodeIndex(:,2)==1,:);
    WaterNodes = NodeIndex(NodeIndex(:,2)==0,:);
    
    SurfNodesHistory = dlmread(SN_FileName);
    
    Neighbors = dlmread(Neighb_FileName);
    
    %Identify Boundary Raft Nodes
    IdentifyEdgeNodes
    
    if it>InputInterval
        CalculateParkingRate;
        CalculateUnparkingRate;
    end
    
    CalculateOutputs
    
    BulkVarphi(ct,:) = CalculateNormalizedVelocity(BulkROI,Bulkt0,Bulktf);
    
    ProtVarphi(ct,:) = CalculateNormalizedVelocity(ProtROI,Prott0,Prottf);
    
    [R,dRdt] = AssignNodeColors;

    TotalTime = it*dt;
    PlotTime(ct,1) = TotalTime;
    
    if ~mod(it,PlotInterval)
        PlotCt = PlotCt+1;
        %Plot Normal Raft
        if PhiMov==1
            figure(1)
            clf
            PlotRaft
            
            F1 = getframe(gcf);
            writeVideo(v1,F1)
            mov1(PlotCt) = F1;
        end
        
        if ROIRetractMov==1
            %Plot Decay
            figure(2)
            clf
            PlotROIDecay(DecayROI)
            
            F2 = getframe(gcf);
            writeVideo(v2,F2)
            mov2(PlotCt) = F2;
            
            RetractionData = [R,dRdt];
            RetractionData = sortrows(RetractionData,1);
            
            if TotalTime>=Decayt0*60 && TotalTime<=Decaytf*60
                mkdir('Retraction Data')
                figure(7)
                clf
                scatter(RetractionData(:,1),RetractionData(:,2),'k');
                
                dlmwrite(['Retraction Data/Trial',num2str(Trial),'.t',...
                    num2str(TotalTime/60,'%.0f'),'.txt'],RetractionData)
            end
        end
        
        if ColorSurfTraffMov==1
            figure(3)
            clf
            PlotColoredSurfTraffic
            
            F3 = getframe(gcf);
            writeVideo(v3,F3)
            mov3(PlotCt) = F3;
        end
        
    end
    
    NodeIndexPrev = NodeIndex;
    SurfNodesPrev = SurfNodes;
    RaftNodesPrev = RaftNodes;
    WaterNodesPrev = WaterNodes;
    
end
if PhiMov==1
    close(v1)
end
if ROIRetractMov==1
    close(v2)
end

if PhiMov==1
    %Plot Bulk ROI
    figure(3)
    clf
    
    PlotRange0 = find(PlotTime(:,1)/60<=Bulkt0,1,'last');
    if isempty(PlotRange0)
        PlotRange0 = 1;
    end
    PlotRangeF = find(PlotTime(:,1)/60>=Bulktf,1,'first');
    NoPts = 100;
    PlotRange = (PlotRange0:floor((PlotRangeF-PlotRange0)/NoPts):PlotRangeF)';
    Time0 = PlotTime(PlotRange0)/60;
    TimeF = PlotTime(PlotRangeF)/60;
    
    PlotBulkVarphi = movmean(BulkVarphi(PlotRange,4),51);
    PlotTimeMean = movmean(PlotTime(PlotRange),51)/60;
    
    scatter(PlotTimeMean,PlotBulkVarphi,'k')
    hold on
    
    xlabel('$t$ (min)','FontSize',FontSize,'Interpreter','latex')
    ylabel('$\varphi_b$','FontSize',FontSize,'Interpreter','latex')
    
    xlim([Time0 TimeF])
    ylim([0 1])
    
    PlotRange0 = find(PlotTime(:,1)/60<=Prott0,1,'last');
    if isempty(PlotRange0)
        PlotRange0 = 1;
    end
    PlotRangeF = find(PlotTime(:,1)/60>=Prottf,1,'first');
    PlotRange = (PlotRange0:PlotRangeF)';
    Time0 = PlotTime(PlotRange0)/60;
    TimeF = PlotTime(PlotRangeF)/60;
    
    PlotProtVarphi = movmean(ProtVarphi(PlotRange,4),51);
    PlotTimeMean = movmean(PlotTime(PlotRange),51)/60;
    
    scatter(PlotTimeMean,PlotProtVarphi,'k^')
    hold on
    
    xlabel('$t$ (min)','FontSize',FontSize,'Interpreter','latex')
    ylabel('$\varphi$','FontSize',FontSize,'Interpreter','latex')
    l = legend({'Bulk','Protrusion'},'FontSize',FontSize/2,'Interpreter','latex',...
        'Location','northwest');
    
    set(gcf,'Color','w')
    % xlim([Time0 TimeF])
    % ylim([0 1])
    
    pbaspect([1 1 1])
    
    saveas(gcf,['Charts/VarphiVsTROI.Trial',num2str(Trial),'.png'])
    
    PlotBulkVarphi(isnan(PlotBulkVarphi),:) = [];
    PlotProtVarphi(isnan(PlotProtVarphi),:) = [];
    
    MeanBulkVarphi = mean(PlotBulkVarphi);
    StdErrorBulkVarphi = mean(PlotBulkVarphi)/sqrt(size(PlotBulkVarphi,1));
    
    MeanProtVarphi = mean(PlotProtVarphi);
    StdErrorProtVarphi = mean(PlotProtVarphi)/sqrt(size(PlotProtVarphi,1));
    
    OutputVarphi = [MeanBulkVarphi StdErrorBulkVarphi MeanProtVarphi StdErrorProtVarphi];
    
    dlmwrite(['AvgVarphi.',num2str(Trial),'.',num2str(ProtNo),'.txt'],OutputVarphi);
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PlotROIDecay(Corners)

global NodeIndex RaftNodes SurfNodes RaftNodeColors NWindow RedExp...
    ParticleSize LengthScale Outline mmPerNode it InputInterval DecayROI

FontSize = 12;
LineWidth = 2;
ShowNodeNo = 0;

BL = Corners(1,:); BR = Corners(2,:);
TR = Corners(3,:); TL = Corners(4,:);
    
%Define Marker Size

BckGrndMarkerSize = 1000*1/NWindow^RedExp;

xs = RaftNodeColors(:,1); 
ys = RaftNodeColors(:,2);
cdatas = zeros(size(xs));

nFree = size(SurfNodes,1);
nRaft = size(RaftNodes,1);
nAdd = nFree+nRaft; %total added plus original raft size
phi_t = nFree/nRaft;

%Plot background lattice
scatter3(xs,ys,cdatas,BckGrndMarkerSize,'k.')
hold on
daspect([1 1 1])
xlim([-NWindow/2 NWindow/2])
ylim([-NWindow/2 NWindow/2])

% Get plot size in pixel units
pos = get(gcf, 'Position'); %// gives x left, y bottom, width, height
FigWidth = pos(3);  %figure width in pixels
MarkerDiam = ParticleSize*FigWidth/NWindow;
MarkerSize = floor((pi/4*MarkerDiam)^2); %in num pixels
if MarkerSize==0
    MarkerSize=1;
end
clf

%Plot a length scale bar
if LengthScale==1
    LengthBarFract = 5;
    LengthBar = NWindow/LengthBarFract;
    LengthBar_mm = round(LengthBar*mmPerNode,-1);
    LengthBar = round(LengthBar_mm/mmPerNode);
    StartBar_x = -LengthBar/2;
    EndBar_x = StartBar_x+LengthBar;
    StartBar_y = -0.9*NWindow/2;
    EndBar_y = StartBar_y;
    
    PlotBarX = [StartBar_x EndBar_x];
    PlotBarY = [StartBar_y EndBar_y];
    
    plot(PlotBarX,PlotBarY,'k','LineWidth',LineWidth)
    hold on
    
    LengthBarText = [num2str(LengthBar_mm,'%.0f'),' mm'];
    text(0,-0.85*NWindow/2,LengthBarText,'FontSize',FontSize,...
        'HorizontalAlignment','center',...
        'Interpreter','Latex')
end

%Plot Raft Background (brown)
h = scatter3(RaftNodes(:,4),RaftNodes(:,5),zeros(size(RaftNodes,1),1),MarkerSize,...
    [0.5 0.25 0],'filled');
if Outline==1
    set(h,'MarkerEdgeColor', 'k');
    hold on
end

%Plot raft lattice that exists within ROI

%use 
    XPlots = RaftNodeColors(:,1);
    YPlots = RaftNodeColors(:,2);
    Colors = RaftNodeColors(:,3);
    
    In = inpolygon(XPlots,YPlots,Corners(:,1),Corners(:,2));
    XPlots = XPlots(In,:);
    YPlots = YPlots(In,:);
    Colors = Colors(In,:);
if it>InputInterval
    XPlots(Colors>1,:) = [];
    YPlots(Colors>1,:) = [];
    Colors(Colors>1,:) = [];
end

h = scatter3(XPlots,YPlots,Colors,MarkerSize,Colors,'filled');
hold on
if Outline==1
    set(h,'MarkerEdgeColor', 'k');
    h.MarkerEdgeAlpha = .8;
    hold on
end
colormap(jet);
c = colorbar;
% label = ylabel(c,'$\dot{R}/\dot{R}_{max}$');
label = ylabel(c,'$dx/dt$ ($\zeta$ s$^{-1}$)');
set(label,'Interpreter','latex')
set(label,'FontSize',FontSize)
caxis([0 0.015])

scatter3(0,0,1,4*MarkerSize,'w+','LineWidth',1.5)

PlotFigureText2(phi_t,nRaft,nFree,nAdd,FontSize);

DBL = DecayROI(1,:);
DBR = DecayROI(2,:);
DTR = DecayROI(3,:);
DTL = DecayROI(4,:);

plot3([DBL(1,1) DBR(1,1)],[DBL(1,2) DBR(1,2)],[100 100],'k','LineWidth',1.5)
hold on
plot3([DBR(1,1) DTR(1,1)],[DBR(1,2) DTR(1,2)],[100 100],'k','LineWidth',1.5)
hold on
plot3([DTR(1,1) DTL(1,1)],[DTR(1,2) DTL(1,2)],[100 100],'k','LineWidth',1.5)
hold on
plot3([DTL(1,1) DBL(1,1)],[DTL(1,2) DBL(1,2)],[100 100],'k','LineWidth',1.5)
hold on

if ShowNodeNo==1
    text(xs,ys,[num2str(NodeIndex(:,1))],'FontSize',5)
end
set(gcf,'color','w');
daspect([1 1 1])
xlim([-NWindow/2 NWindow/2])
ylim([-NWindow/2 NWindow/2])
xticks([])
yticks([])
box on

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PlotDecay

global NodeIndex RaftNodes SurfNodes RaftNodeColors NWindow RedExp...
    ParticleSize LengthScale Outline mmPerNode it InputInterval

FontSize = 12;
LineWidth = 2;
ShowNodeNo = 0;

%Define Marker Size

BckGrndMarkerSize = 1000*1/NWindow^RedExp;

xs = RaftNodeColors(:,1); 
ys = RaftNodeColors(:,2);
cdatas = zeros(size(xs));

nFree = size(SurfNodes,1);
nRaft = size(RaftNodes,1);
nAdd = nFree+nRaft; %total added plus original raft size
phi_t = nFree/nRaft;

%Plot background lattice
scatter3(xs,ys,cdatas,BckGrndMarkerSize,'k.')
hold on
daspect([1 1 1])
xlim([-NWindow/2 NWindow/2])
ylim([-NWindow/2 NWindow/2])

% Get plot size in pixel units
pos = get(gcf, 'Position'); %// gives x left, y bottom, width, height
FigWidth = pos(3);  %figure width in pixels
MarkerDiam = ParticleSize*FigWidth/NWindow;
MarkerSize = floor((pi/4*MarkerDiam)^2); %in num pixels
if MarkerSize==0
    MarkerSize=1;
end
clf

%Plot a length scale bar
if LengthScale==1
    LengthBarFract = 5;
    LengthBar = NWindow/LengthBarFract;
    LengthBar_mm = round(LengthBar*mmPerNode,-1);
    LengthBar = round(LengthBar_mm/mmPerNode);
    StartBar_x = -LengthBar/2;
    EndBar_x = StartBar_x+LengthBar;
    StartBar_y = -0.9*NWindow/2;
    EndBar_y = StartBar_y;
    
    PlotBarX = [StartBar_x EndBar_x];
    PlotBarY = [StartBar_y EndBar_y];
    
    plot(PlotBarX,PlotBarY,'k','LineWidth',LineWidth)
    hold on
    
    LengthBarText = [num2str(LengthBar_mm,'%.0f'),' mm'];
    text(0,-0.85*NWindow/2,LengthBarText,'FontSize',FontSize,...
        'HorizontalAlignment','center',...
        'Interpreter','Latex')
end

%Plot raft lattice
    XPlots = RaftNodeColors(:,1);
    YPlots = RaftNodeColors(:,2);
    Colors = RaftNodeColors(:,3);
if it>InputInterval
    XPlots(Colors>1,:) = [];
    YPlots(Colors>1,:) = [];
    Colors(Colors>1,:) = [];
end

h = scatter3(XPlots,YPlots,Colors,MarkerSize,Colors,'filled');
hold on
if Outline==1
    set(h,'MarkerEdgeColor', 'k');
    hold on
end
colormap(jet);
c = colorbar;
label = ylabel(c,'$\dot{R}/\dot{R}_{max}$');
set(label,'Interpreter','latex')
set(label,'FontSize',FontSize)
caxis([0 1])

scatter3(0,0,1,4*MarkerSize,'w+','LineWidth',1.5)

PlotFigureText2(phi_t,nRaft,nFree,nAdd,FontSize);

if ShowNodeNo==1
    text(xs,ys,[num2str(NodeIndex(:,1))],'FontSize',5)
end
set(gcf,'color','w');
daspect([1 1 1])
xlim([-NWindow/2 NWindow/2])
ylim([-NWindow/2 NWindow/2])
xticks([])
yticks([])
box on

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Dist,Speed] = AssignNodeColors

global RaftNodes RaftNodesPrev RaftNodeColors it InputInterval dt MaxSpeed
  
Dist = 0; Speed = 0;

if it==InputInterval
    X = RaftNodes(:,4);
    Y = RaftNodes(:,5);
    RaftNodeColors = [X Y zeros(size(X,1),1)];
    MaxSpeed = 'NA';
else
    OldRaftNodes = RaftNodesPrev(:,1);
    NewRaftNodes = RaftNodes(:,1);
    
    ConservedRaftNodeInfo = ...
        OldRaftNodes(ismember(OldRaftNodes(:,1),NewRaftNodes(:,1)),:);
    ConservedRaftNodes = unique(ConservedRaftNodeInfo(:,1),'rows');
    
    [~,UniqueRows] = unique(RaftNodes(:,1),'rows');
    [~,UniqueRowsPrev] = unique(RaftNodesPrev(:,1),'rows');
    
    UniqueRaftNodes = RaftNodes(UniqueRows,:);
    UniqueRaftNodesPrev = RaftNodesPrev(UniqueRowsPrev,:);
    
    X = UniqueRaftNodes(ismember(UniqueRaftNodes(:,1),ConservedRaftNodes),4);
    X0 = ...
        UniqueRaftNodesPrev(ismember(UniqueRaftNodesPrev(:,1),ConservedRaftNodes),4);
    Y = UniqueRaftNodes(ismember(UniqueRaftNodes(:,1),ConservedRaftNodes),5);
    Y0 = ...
        UniqueRaftNodesPrev(ismember(UniqueRaftNodesPrev(:,1),ConservedRaftNodes),5);
    
    Disp = [X,Y]-[X0,Y0];
    
    AvgPosX = mean([X,X0],2);
    AvgPosY = mean([Y,Y0],2);
    MeanPos = [AvgPosX,AvgPosY];
    UnitVec = -MeanPos./vecnorm(MeanPos,2,2);
    
    DeltaR = dot(Disp,UnitVec,2);
    
    Dist = vecnorm([X,Y],2,2);
    DeltaT = InputInterval*dt;
    
    Speed = DeltaR/DeltaT;    %zeta/s
    MaxSpeed = max(Speed);
    TempSpeed = sort(Speed);
    MeanSpeed = mean(Speed);
    while MaxSpeed>2.5*MeanSpeed
        TempSpeed(end,:) = [];
        MaxSpeed = max(TempSpeed);
    end
        
%     RaftNodeColors = [X Y Speed/MaxSpeed];
    RaftNodeColors = [X Y Speed];
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PlotRaft 

global NodeIndex SurfNodes RaftNodes RedExp EdgeNodes...
    mmPerNode LengthScale Outline NWindow ParticleSize...
    NewNeighDisThresh SurfNodesHistory BulkROI Bulkt0 Bulktf ProtROI...
    Prott0 Prottf it dt

FontSize = 12;
LineWidth = 2;
PlotBoundNodes = 0;
ShowNodeNo = 0;

SurfPos = [SurfNodes(:,4) SurfNodes(:,5)];
RaftPos = [RaftNodes(:,4) RaftNodes(:,5)];

%Define Marker Size

BckGrndMarkerSize = 1000*1/NWindow^RedExp;

xs = NodeIndex(:,4); ys = NodeIndex(:,5);

nFree = size(SurfNodes,1);
nRaft = size(RaftNodes,1);
nAdd = nFree+nRaft; %total added plus original raft size
phi_t = nFree/nRaft;

%Plot background lattice
scatter(xs,ys,BckGrndMarkerSize,'k.')
hold on
daspect([1 1 1])
xlim([-NWindow/2 NWindow/2])
ylim([-NWindow/2 NWindow/2])

% Get plot size in pixel units
pos = get(gcf, 'Position'); %// gives x left, y bottom, width, height
FigWidth = pos(3);  %figure width in pixels
MarkerDiam = ParticleSize*FigWidth/NWindow;
MarkerSize = floor((pi/4*MarkerDiam)^2); %in num pixels
if MarkerSize==0
    MarkerSize=1;
end
clf

%Plot a length scale bar
if LengthScale==1
    LengthBarFract = 5;
    LengthBar = NWindow/LengthBarFract;
    LengthBar_mm = round(LengthBar*mmPerNode,-1);
    LengthBar = round(LengthBar_mm/mmPerNode);
    StartBar_x = -LengthBar/2;
    EndBar_x = StartBar_x+LengthBar;
    StartBar_y = -0.9*NWindow/2;
    EndBar_y = StartBar_y;
    
    PlotBarX = [StartBar_x EndBar_x];
    PlotBarY = [StartBar_y EndBar_y];
    
    plot(PlotBarX,PlotBarY,'k','LineWidth',LineWidth)
    hold on
    
    LengthBarText = [num2str(LengthBar_mm,'%.0f'),' mm'];
    text(0,-0.85*NWindow/2,LengthBarText,'FontSize',FontSize,...
        'HorizontalAlignment','center',...
        'Interpreter','Latex')
end

%Plot raft lattice
h = scatter(RaftPos(:,1),RaftPos(:,2),MarkerSize,'c','filled');
hold on
if Outline==1
    set(h,'MarkerEdgeColor', 'k');
    hold on
end

%Plot edge nodes
if PlotBoundNodes==1
    %NEED TO WRITE CODE TO FIND EDGE NODES%
        %1. FIND NEIGHBORS
        %2. FIND RAFT NODES WITH WATER NEIGHBORS
    EdgeNodesXY = NodeIndex(ismember(NodeIndex(:,1),EdgeNodes),4:5);
    h = scatter(EdgeNodesXY(:,1),EdgeNodesXY(:,2),MarkerSize,[0 0.5 0.5],'filled');
    hold on
    if Outline==1
        scatter(EdgeNodesXY(:,1),EdgeNodesXY(:,2),MarkerSize,'k')
        hold on
        set(h,'MarkerEdgeColor', 'k');
        hold on
    end
end

% Plot previous surface particles
OldSurfPos = SurfNodesHistory(:,5:6);
NewSurfPos = SurfNodesHistory(:,2:3);

Disp = NewSurfPos-OldSurfPos;
Dist = vecnorm(Disp,2,2);

PlotDist = Dist;
PlotDist(Dist>=0.85*NewNeighDisThresh,:) = 0.85*NewNeighDisThresh;

UnitVec = Disp./Dist;

AbdomenPos = NewSurfPos-UnitVec*1.*PlotDist;
ThoraxPos = NewSurfPos-UnitVec*0.5.*PlotDist;
AbdomenPos(AbdomenPos(:,1)==0,:)=[];
ThoraxPos(ThoraxPos(:,1)==0,:)=[];

h1 = scatter(AbdomenPos(:,1),AbdomenPos(:,2),MarkerSize,'r','filled','MarkerFaceAlpha',0.6);
if Outline==1
    set(h1,'MarkerEdgeColor', 'k');
    hold on
end
h2 = scatter(ThoraxPos(:,1),ThoraxPos(:,2),MarkerSize/2,'r','filled','MarkerFaceAlpha',0.6);
if Outline==1
    set(h2,'MarkerEdgeColor', 'k');
    hold on
end

%Plot current surface particles
h = scatter(SurfPos(:,1),SurfPos(:,2),MarkerSize,'r','filled');
hold on
if Outline==1
    set(h,'MarkerEdgeColor', 'k');     
    hold on
end

%Plot ROIs
Time = it*dt/60;    %min

BBL = BulkROI(1,:);
BBR = BulkROI(2,:);
BTR = BulkROI(3,:);
BTL = BulkROI(4,:);

PBL = ProtROI(1,:);
PBR = ProtROI(2,:);
PTR = ProtROI(3,:);
PTL = ProtROI(4,:);

if Bulkt0<=Time && Time<=Bulktf
    plot([BBL(1,1) BBR(1,1)],[BBL(1,2) BBR(1,2)],'k','LineWidth',1.5)
    hold on
    plot([BBR(1,1) BTR(1,1)],[BBR(1,2) BTR(1,2)],'k','LineWidth',1.5)
    hold on
    plot([BTR(1,1) BTL(1,1)],[BTR(1,2) BTL(1,2)],'k','LineWidth',1.5)
    hold on
    plot([BTL(1,1) BBL(1,1)],[BTL(1,2) BBL(1,2)],'k','LineWidth',1.5)
    hold on
end

if Prott0<=Time && Time<=Prottf
    plot([PBL(1,1) PBR(1,1)],[PBL(1,2) PBR(1,2)],'k','LineWidth',1.5)
    hold on
    plot([PBR(1,1) PTR(1,1)],[PBR(1,2) PTR(1,2)],'k','LineWidth',1.5)
    hold on
    plot([PTR(1,1) PTL(1,1)],[PTR(1,2) PTL(1,2)],'k','LineWidth',1.5)
    hold on
    plot([PTL(1,1) PBL(1,1)],[PTL(1,2) PBL(1,2)],'k','LineWidth',1.5)
    hold on
end

%Plot figure text
PlotFigureText1(phi_t,nRaft,nFree,nAdd,FontSize);

if ShowNodeNo==1
    text(xs,ys,[num2str(NodeIndex(:,1))],'FontSize',5)
end
set(gcf,'color','w');
daspect([1 1 1])
xlim([-NWindow/2 NWindow/2])
ylim([-NWindow/2 NWindow/2])
xticks([])
yticks([])
box on

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PlotDomain

global NodeIndex SurfNodes RaftNodes RedExp...
    Outline NWindow ParticleSize SurfNodesHistory...
    NewNeighDisThresh 

FontSize = 12;
LineWidth = 2;
PlotBoundNodes = 0;
ShowNodeNo = 0;

SurfPos = [SurfNodes(:,4) SurfNodes(:,5)];
RaftPos = [RaftNodes(:,4) RaftNodes(:,5)];

%Define Marker Size

BckGrndMarkerSize = 1000*1/NWindow^RedExp;

xs = NodeIndex(:,4); ys = NodeIndex(:,5);

nFree = size(SurfNodes,1);
nRaft = size(RaftNodes,1);
nAdd = nFree+nRaft; %total added plus original raft size
phi_t = nFree/nRaft;

%Plot background lattice
scatter(xs,ys,BckGrndMarkerSize,'k.')
hold on
daspect([1 1 1])
xlim([-NWindow/2 NWindow/2])
ylim([-NWindow/2 NWindow/2])

% Get plot size in pixel units
pos = get(gcf, 'Position'); %// gives x left, y bottom, width, height
FigWidth = pos(3);  %figure width in pixels
MarkerDiam = ParticleSize*FigWidth/NWindow;
MarkerSize = floor((pi/4*MarkerDiam)^2); %in num pixels
if MarkerSize==0
    MarkerSize=1;
end

%Plot raft lattice
h = scatter(RaftPos(:,1),RaftPos(:,2),MarkerSize,'c','filled');
hold on
if Outline==1
    set(h,'MarkerEdgeColor', 'k');
    hold on
end

% Plot previous surface particles
OldSurfPos = SurfNodesHistory(:,5:6);
NewSurfPos = SurfNodesHistory(:,2:3);

Disp = NewSurfPos-OldSurfPos;
Dist = vecnorm(Disp,2,2);

PlotDist = Dist;
PlotDist(Dist>=0.85*NewNeighDisThresh,:) = 0.85*NewNeighDisThresh;

UnitVec = Disp./Dist;

AbdomenPos = NewSurfPos-UnitVec*1.*PlotDist;
ThoraxPos = NewSurfPos-UnitVec*0.5.*PlotDist;
AbdomenPos(AbdomenPos(:,1)==0,:)=[];
ThoraxPos(ThoraxPos(:,1)==0,:)=[];

h1 = scatter(AbdomenPos(:,1),AbdomenPos(:,2),MarkerSize,'r','filled','MarkerFaceAlpha',0.6);
if Outline==1
    set(h1,'MarkerEdgeColor', 'k');
    hold on
end
h2 = scatter(ThoraxPos(:,1),ThoraxPos(:,2),MarkerSize/2,'r','filled','MarkerFaceAlpha',0.6);
if Outline==1
    set(h2,'MarkerEdgeColor', 'k');
    hold on
end

%Plot current surface particles
h = scatter(SurfPos(:,1),SurfPos(:,2),MarkerSize,'r','filled');
hold on
if Outline==1
    set(h,'MarkerEdgeColor', 'k');     
    hold on
end

PlotFigureText0(FontSize)

if ShowNodeNo==1
    text(xs,ys,[num2str(NodeIndex(:,1))],'FontSize',5)
end
set(gcf,'color','w');
daspect([1 1 1])
xlim([-NWindow/2 NWindow/2])
ylim([-NWindow/2 NWindow/2])
xticks([-NWindow/2:10:NWindow/2])
yticks([-NWindow/2:10:NWindow/2])
grid on
ax = gca;
ax.GridLineStyle = '-';
ax.GridColor = 'k';
box on

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PlotFigureText0(FontSize)

global it dt  

Time = it*dt;

%Top Left
xlabel(['$t$ = ',num2str(Time/60,'%.2f'),' min'],'FontSize',FontSize,...
    'Interpreter','Latex')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PlotFigureText1(phi_t,nRaft,nFree,nAdd,FontSize)

global NWindow it dt parkrate unparkrate R0 WorkLengthScale...
    PlotTime RadNeighbors TrajectoryNoise tau_dec

Time = it*dt;
TimeRange = 60;
TimeDomain = find(PlotTime(:,1)<=(Time-TimeRange));

if isempty(TimeDomain)
    AvgAlpha = mean(parkrate);
    AvgDelta = mean(unparkrate);
else
    AvgAlpha = mean(parkrate(TimeDomain(end):end));
    AvgDelta = mean(unparkrate(TimeDomain(end):end));
end

%Top Left
text(-0.95*NWindow/2,0.9*NWindow/2,...
    ['$t$ = ',num2str(Time/60,'%.0f'),' min'],'FontSize',FontSize,...
    'Interpreter','Latex')
text(-0.95*NWindow/2,0.8*NWindow/2,...
    ['$\langle \alpha(t) \rangle$ = ',num2str(AvgAlpha,'%.4f'),...
     ' $n_{p}$ min$^{-1}$ $\zeta^{-2}$'],...
     'FontSize',FontSize,'Interpreter','latex')
text(-0.95*NWindow/2,0.7*NWindow/2,...
    ['$\langle \delta(t) \rangle$ = ',num2str(AvgDelta,'%.4f'),...
     ' $n_{u}$ min$^{-1}$ $\zeta^{-2}$'],...
     'FontSize',FontSize,'Interpreter','latex')

%Top Right
text(0.95*NWindow/2,0.9*NWindow/2,...
    ['$R_0$ = ',num2str(RadNeighbors),' $\zeta$'],'FontSize',FontSize,...
    'HorizontalAlignment','right','Interpreter','Latex')
text(0.95*NWindow/2,0.8*NWindow/2,...
    ['$\xi$ = $\pm$',num2str(TrajectoryNoise*360/2),'$^{\circ}$'],'FontSize',FontSize,...
    'HorizontalAlignment','right','Interpreter','Latex')

%Bottom Left
text(-0.95*NWindow/2,-0.9*NWindow/2,...
    ['$\phi_{ss}$ = ',num2str(0.238,'%.2f')],...
    'FontSize',FontSize,'Interpreter','Latex')
text(-0.95*NWindow/2,-0.8*NWindow/2,...
    ['$\phi(t)$ = ',num2str(phi_t,'%.2f')],...
    'FontSize',FontSize,'Interpreter','Latex')
text(-0.95*NWindow/2,-0.7*NWindow/2,...
    ['$\phi_{0}$ = ',num2str(1,'%.2f')],...
    'FontSize',FontSize,'Interpreter','Latex')

%Bottom Right
text(0.95*NWindow/2,-0.7*NWindow/2,['$N_{surf}$ = ',num2str(nFree)],...
    'FontSize',FontSize,'HorizontalAlignment','right','Interpreter','Latex')
text(0.95*NWindow/2,-0.8*NWindow/2,['$N_{raft}$ = ',num2str(nRaft)],...
    'FontSize',FontSize,'HorizontalAlignment','right','Interpreter','Latex')
text(0.95*NWindow/2,-0.9*NWindow/2,['$N_{add}$ = ',num2str(nAdd)],...
    'FontSize',FontSize,'HorizontalAlignment','right','Interpreter','Latex')

xlabel(['$R_0$ = ',num2str(R0,'%.2f'),'$\zeta$'...
    ', $L$ = ',num2str(WorkLengthScale,'%.2f'),'$\zeta$',...
    ', $\tau$ = ',num2str(tau_dec/60,'%.0f'),' min'],...
    'FontSize',FontSize,'Interpreter','latex')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PlotFigureText2(phi_t,nRaft,nFree,nAdd,FontSize)

global NWindow it dt parkrate unparkrate R0 WorkLengthScale...
    PlotTime RadNeighbors TrajectoryNoise tau_dec MaxSpeed InputInterval

Time = it*dt;
TimeRange = 60;
TimeDomain = find(PlotTime(:,1)<=(Time-TimeRange));

if isempty(TimeDomain)
    AvgAlpha = mean(parkrate);
    AvgDelta = mean(unparkrate);
else
    AvgAlpha = mean(parkrate(TimeDomain(end):end));
    AvgDelta = mean(unparkrate(TimeDomain(end):end));
end

%Top Left
text(-0.95*NWindow/2,0.9*NWindow/2,...
    ['$t$ = ',num2str(Time/60,'%.0f'),' min'],'FontSize',FontSize,...
    'Interpreter','Latex')
text(-0.95*NWindow/2,0.8*NWindow/2,...
    ['$\langle \alpha(t) \rangle$ = ',num2str(AvgAlpha,'%.4f'),...
     ' $n_{p}$ min$^{-1}$ $\zeta^{-2}$'],...
     'FontSize',FontSize,'Interpreter','latex')
text(-0.95*NWindow/2,0.7*NWindow/2,...
    ['$\langle \delta(t) \rangle$ = ',num2str(AvgDelta,'%.4f'),...
     ' $n_{u}$ min$^{-1}$ $\zeta^{-2}$'],...
     'FontSize',FontSize,'Interpreter','latex')

%Top Right
text(0.95*NWindow/2,0.9*NWindow/2,...
    ['$R_0$ = ',num2str(RadNeighbors),' $\zeta$'],'FontSize',FontSize,...
    'HorizontalAlignment','right','Interpreter','Latex')
text(0.95*NWindow/2,0.8*NWindow/2,...
    ['$\xi$ = $\pm$',num2str(TrajectoryNoise*360/2),'$^{\circ}$'],'FontSize',FontSize,...
    'HorizontalAlignment','right','Interpreter','Latex')

%Bottom Left
text(-0.95*NWindow/2,-0.9*NWindow/2,...
    ['$\phi_{ss}$ = ',num2str(0.238,'%.2f')],...
    'FontSize',FontSize,'Interpreter','Latex')
text(-0.95*NWindow/2,-0.8*NWindow/2,...
    ['$\phi(t)$ = ',num2str(phi_t,'%.2f')],...
    'FontSize',FontSize,'Interpreter','Latex')
text(-0.95*NWindow/2,-0.7*NWindow/2,...
    ['$\phi_{0}$ = ',num2str(1,'%.2f')],...
    'FontSize',FontSize,'Interpreter','Latex')

%Bottom Right
text(0.95*NWindow/2,-0.7*NWindow/2,['$N_{surf}$ = ',num2str(nFree)],...
    'FontSize',FontSize,'HorizontalAlignment','right','Interpreter','Latex')
text(0.95*NWindow/2,-0.8*NWindow/2,['$N_{raft}$ = ',num2str(nRaft)],...
    'FontSize',FontSize,'HorizontalAlignment','right','Interpreter','Latex')
text(0.95*NWindow/2,-0.9*NWindow/2,['$N_{add}$ = ',num2str(nAdd)],...
    'FontSize',FontSize,'HorizontalAlignment','right','Interpreter','Latex')

if it==InputInterval
    xlabel(['$R_0$ = ',num2str(R0,'%.2f'),'$\zeta$'...
        ', $L$ = ',num2str(WorkLengthScale,'%.2f'),'$\zeta$',...
        ', $\dot{\epsilon}_{ij}$ = ',num2str(1/tau_dec*60,'%.2f'),' min$^{-1}$',...
        ', $\dot{R}_{max}$ = NA $\zeta$ min$^{-1}$'],...
        'FontSize',FontSize,'Interpreter','latex')
else
    xlabel(['$R_0$ = ',num2str(R0,'%.2f'),'$\zeta$'...
        ', $L$ = ',num2str(WorkLengthScale,'%.2f'),'$\zeta$',...
        ', $\dot{\epsilon}_{ij}$ = ',num2str(1/tau_dec*60,'%.2f'),' min$^{-1}$',...
        ', $\dot{R}_{max}$ = ',num2str(MaxSpeed*60,'%.2f'),' $\zeta$ min$^{-1}$'],...
        'FontSize',FontSize,'Interpreter','latex')
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CalculateParkingRate

%Count number of new raft nodes at each time 

global RaftNodesPrev RaftNodes dt ct parkrate InputInterval

%Calculate Avg Area Between Frames
CurrentArea = size(RaftNodes,1);
PrevArea = size(RaftNodesPrev,1);

AvgArea = mean([CurrentArea PrevArea]);

%ID New Water Nodes
NewRaftNodes = RaftNodes(~ismember(RaftNodes(:,1),RaftNodesPrev(:,1)),1);

NoParks = size(NewRaftNodes,1);

parkrate(ct,1) = NoParks/AvgArea/(dt*InputInterval)*60;    %parks /zeta^2 /min 

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CalculateUnparkingRate

global RaftNodesPrev SurfNodes RaftNodes  dt ct unparkrate InputInterval

%Calculate Avg Area Between Frames
CurrentArea = size(RaftNodes,1);
PrevArea = size(RaftNodesPrev,1);

AvgArea = mean([CurrentArea PrevArea]);

Check1 = SurfNodes(:,1);
Check2 = RaftNodesPrev(:,1);

%ID New Water Nodes
UnparkedNodes = RaftNodesPrev(~ismember(RaftNodesPrev(:,1),RaftNodes(:,1)),1);

NoParks = size(UnparkedNodes,1);

unparkrate(ct,1) = NoParks/AvgArea/(dt*InputInterval)*60;    %parks /zeta^2 /min 

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function IdentifyEdgeNodes

%Function identifies edge particles in raft by singling out any particles
%with a neighbor that eq. 0 in connect matrix

global NodeIndex EdgeNodes Neighbors MaxDOF

%Check Dist. Between All Nodes & Their Neighbors

%First, crop data to only look at Raft Nodes
NeighborsTemp = Neighbors;
NodeRorW = NodeIndex(ismember(NodeIndex(:,1),NeighborsTemp(:,end)),2);
NeighborsTemp(NodeRorW==0,:) = [];

EdgeNodes = [];

for i=1:MaxDOF
    %Define Neighbors and Nodes
    NeighborSet = NeighborsTemp(:,i);
    NodeSet = NeighborsTemp(:,end);
    
    %Identify the pertinent information of Nodes and Neighb in NodeIndex
    [~,NodeRows] = ismember(NodeSet,NodeIndex(:,1));
    [~,NeighbRows] = ismember(NeighborSet,NodeIndex(:,1));
    
    %Ensure the Node has an actual Neighbor - crop both sets of rows where
        %NeighbRows==0
    NodeRows(NeighbRows==0,:) = [];
    NeighbRows(NeighbRows==0,:) = [];
    
    NodeInfo = NodeIndex(NodeRows,:);
    NeighbInfo = NodeIndex(NeighbRows,:);
    
    %Crop out raft neighbors 
    NodeInfo(NeighbInfo(:,2)==1,:) = [];
    NeighbInfo(NeighbInfo(:,2)==1,:) = [];
    
    %Check dist. between nodes and neighbors
    Disp = NodeInfo(:,4:5)-NeighbInfo(:,4:5);
    Dist = vecnorm(Disp,2,2);
    
    EdgeNodes_temp = NodeInfo(Dist<1.4,1);
    EdgeNodes = [EdgeNodes;EdgeNodes_temp];   
end
EdgeNodes = unique(EdgeNodes);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculat the local pressure on surf particle from neighboring surf particles
function [NetPressure,AlphaPress] = CalculateLocalPressure(NOI)

global NodeIndex PressureRadius

%ID other Raft Nodes
OtherRaftNodesInfo = NodeIndex(NodeIndex(:,2)==1,:);
OtherRaftNodesInfo(OtherRaftNodesInfo(:,1)==NOI,:)=[];
OtherRaftNodesPos = OtherRaftNodesInfo(:,4:5);

%Find Neigboring Raft Nodes
NOI_Pos = NodeIndex(NodeIndex(:,1)==NOI,4:5); 
PosMatrix = ones(size(OtherRaftNodesPos)).*NOI_Pos;

Disp = PosMatrix-OtherRaftNodesPos; %Disp from Other Nodes to Node of Interest (NOI)
Dist = vecnorm(Disp,2,2);
NeighboringRaftNodes = OtherRaftNodesInfo(Dist<=PressureRadius,:);
NeighboringSurfNodes = NeighboringRaftNodes(NeighboringRaftNodes(:,3)==1,:);

%Crop displacement and distance vectors to other nodes to only include
    %neighbors
Disp(Dist>PressureRadius,:) = [];
Dist(Dist>PressureRadius,:) = [];

%Calculate Weighting Function for Each of Neighboring Nodes
WeightType=0;
if WeightType==0
    Weight = (PressureRadius-Dist)/PressureRadius;
elseif WeightType==1
    Weight = exp(-Dist/PressureRadius);
end
Weight(NeighboringRaftNodes(:,3)==0,:) = 0;    %Set weight of empty neighbors to zero

%Calculate displacement unit vectors from neighbors to NOI
D_hat = Disp./Dist;

%Calculate Pressure
Pressure = Weight.*D_hat;
NetPressure = sum(Pressure,1);

AlphaPress = atan2d(NetPressure(2),NetPressure(1));
if AlphaPress<0
    AlphaPress = 360+AlphaPress;
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculate the updated velocity of each particle as predicted by Vicsek
    %model
function [VicsekVelocity,AlphaNew] = CalculateLocalVelocity(NOI)

global SurfNodes RadNeighbors TrajectoryNoise dt

SurfData = SurfNodes(SurfNodes(:,1)==NOI,:); 
    %[CurrentNode XCurrent YCurrent PreviousNode XPrev YPrev];

CurrentPos = [SurfData(:,2) SurfData(:,3)];
PrevPos = [SurfData(:,5) SurfData(:,6)];

OtherSurfNodeInfo = SurfNodes(~ismember(SurfNodes(:,1),NOI),:);

OtherSurfNodes = OtherSurfNodeInfo(:,1);
OtherSurfNodePos = OtherSurfNodeInfo(:,2:3);
OtherSurfNodePrevPos = OtherSurfNodeInfo(:,5:6);

%Find Dist to Other Surf Particles
PosMatrix = ones(size(OtherSurfNodePos)).*CurrentPos;
Dist = vecnorm([PosMatrix-OtherSurfNodePos],2,2);

%ID Particles Within Range
SortMatrix = [Dist OtherSurfNodes OtherSurfNodePos OtherSurfNodePrevPos];

%Calculate the New Velocity Direction Based on Neighbors' Velocities
%Isolate Information for Influencing Particles Including The Node
%of Interest
SurfNeighborsInfo = ...
    SortMatrix(SortMatrix(:,1)<=RadNeighbors,:);
SurfNeighborsInfo(end+1,:) = [0 NOI CurrentPos PrevPos];
%If the Previous Pos. is [0 0], define it as the current position
XZeros = SurfNeighborsInfo(SurfNeighborsInfo(:,5)==0,2);
YZeros = SurfNeighborsInfo(SurfNeighborsInfo(:,6)==0,2);
UpdateNodes = XZeros(ismember(XZeros,YZeros),1);
SurfNeighborsInfo(ismember(SurfNeighborsInfo(:,2),UpdateNodes),5:6) = ...
    SurfNeighborsInfo(ismember(SurfNeighborsInfo(:,2),UpdateNodes),3:4);

%Define Weighting Factor for Influence Based on Distance & Calc. Wtd
%Velocities
DistanceWeightingFactor = (RadNeighbors-SurfNeighborsInfo(:,1))/RadNeighbors;

SurfNeighborsVel = (SurfNeighborsInfo(:,3:4)-SurfNeighborsInfo(:,5:6))/dt;

if ~any(SurfNeighborsVel(:)) %if particle wasn't moving neither was its neighbors
    RandOrPress=1;
    if RandOrPress==0 %Randomly select new direction
        V_0 = 1;    %The set it's effective speed to 1
        Alpha = rand*360;  %And pick a random direction for it to be facing
    else %Pick direction based on pressure
        [NetPressure,PressAngle] = CalculateLocalPressure(NOI);
        if NetPressure(1)==0 && NetPressure(2)==0
            V_0 = 1;    %The set it's effective speed to 1
            Alpha = rand*360;  %And pick a random direction for it to be facing
        else
            Alpha = PressAngle;
            V_0 = 1;
        end
    end
    
    if Alpha==0
        SurfNeighborsVel(end,1:2) = V_0*[1 0];
    elseif 0<Alpha && Alpha<90
        SurfNeighborsVel(end,1:2) = V_0*[cosd(Alpha) sind(Alpha)];
    elseif Alpha==90
        SurfNeighborsVel(end,1:2) = V_0*[0 1];
    elseif 90<Alpha && Alpha<180
        SurfNeighborsVel(end,1:2) = V_0*[-cosd(180-Alpha) sind(180-Alpha)];
    elseif Alpha==180
        SurfNeighborsVel(end,1:2) = V_0*[-1 0];
    elseif 180<Alpha && Alpha<270
        SurfNeighborsVel(end,1:2) = V_0*[-cosd(180+Alpha) -sind(180+Alpha)];
    elseif Alpha==270
        SurfNeighborsVel(end,1:2) = V_0*[0 -1];
    elseif 270<Alpha
        SurfNeighborsVel(end,1:2) = V_0*[cosd(360-Alpha) -sind(360-Alpha)];
    end   
end

WtdVelocities = DistanceWeightingFactor.*SurfNeighborsVel;
WtdVelocity = sum(WtdVelocities,1);

%Find Corresponding Unit Vector
NormWtdVelocity = WtdVelocity/(vecnorm(WtdVelocity,2,2));

%Speed of particle at previous timestep
if any(SurfNeighborsVel(end,:))
    V_0 = vecnorm(SurfNeighborsVel(end,:),2,2);
else
    V_0 = 1;
end

%Updated Velocity of Particle
V = V_0*NormWtdVelocity;
Alpha = atan2(V(2),V(1));  %measured from x-axis with (-)=> Quad. III or IV

V_mag = V_0;

%Pick a Random Rotation amt. Based on TrajectoryNoise
eta = (rand-0.5)*TrajectoryNoise; %Random value within +0.5/-0.5 of noise range
if V(1)~=0 && V(2)~=0 && ~isnan(V(1))
    dAlpha = eta*2*pi;
    dAlphad = dAlpha*180/pi;
end

Alphad = Alpha*180/pi;
%Update Velocity With Noise
AlphaNew = Alphad+dAlphad;
if AlphaNew>=360
    AlphaNew = AlphaNew-360;
elseif AlphaNew<0
    AlphaNew = AlphaNew+360;
end
if AlphaNew==0
    VicsekVelocity = V_mag*[1 0];
elseif 0<AlphaNew && AlphaNew<90
    VicsekVelocity = V_mag*[cosd(AlphaNew) sind(AlphaNew)];
elseif AlphaNew==90
    VicsekVelocity = V_mag*[0 1];
elseif 90<AlphaNew && AlphaNew<180
    VicsekVelocity = V_mag*[-cosd(180-AlphaNew) sind(180-AlphaNew)];
elseif AlphaNew==180
    VicsekVelocity = V_mag*[-1 0];
elseif 180<AlphaNew && AlphaNew<270
    VicsekVelocity = V_mag*[-cosd(180+AlphaNew) -sind(180+AlphaNew)];
elseif AlphaNew==270
    VicsekVelocity = V_mag*[0 -1];
elseif 270<AlphaNew
    VicsekVelocity = V_mag*[cosd(360-AlphaNew) -sind(360-AlphaNew)];
end

DNew = VicsekVelocity*dt; %Relative displacement - use to calc best DOF using a dist. check

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function HealTheDomainOld
%This function ID's nodes with neighbor indices that are zero and allows
%for repairing between any other particle within Raft NewNeighDisThresh

%It functionally redefines the neighbors after nodes have been removed

global NodeIndex Neighbors NewNeighDisThresh GridSize ct

MaxNoNeighb = size(Neighbors,2)-1;

NoNeighbors = size(Neighbors,1);

%Find Neighbors Rows Containing Zeros in one of their Neighbors Slots and
    %Not Belonging to Edge of Domain
NeighborsWithZeros = Neighbors(sum(Neighbors==0,2)>=1,:);

%Exclude the Edge Nodes of the Domain for Expediency
MaxX = max(NodeIndex(:,4));
MinX = min(NodeIndex(:,4));
MaxY = max(NodeIndex(:,5));
MinY = min(NodeIndex(:,5));

RNodes = NodeIndex(NodeIndex(:,4)>(MaxX-1.5*GridSize),1);
LNodes = NodeIndex(NodeIndex(:,4)<(MinX+1.5*GridSize),1);
TNodes = NodeIndex(NodeIndex(:,5)>(MaxY-1.5*GridSize),1);
BNodes = NodeIndex(NodeIndex(:,5)<(MinY+1.5*GridSize),1);
EdgeNodes = unique([RNodes;LNodes;TNodes;BNodes]);

NeighborsWithZeros(ismember(NeighborsWithZeros(:,end),EdgeNodes),:) = [];

%Scramble the Neighbors with Zeros
NoNeighborsWithZeros = size(NeighborsWithZeros,1);
ScrambleIndex = randperm(NoNeighborsWithZeros,NoNeighborsWithZeros)';
NeighborsWithZeros = NeighborsWithZeros(ScrambleIndex,:);
        
if mod(ct,4)~=0 %only check raft particles every other iteration
    OnlyCheckRaft=1;
    disp('Healing the Raft')
else
    OnlyCheckRaft=0;    %Do not heal nodes that aren't connected to raft
    disp('Healing the Domain')
end

if OnlyCheckRaft==1
   RaftNodes = NodeIndex(NodeIndex(:,2)==1,1);
   RaftNeighborsWithZeros = ...
       NeighborsWithZeros(ismember(NeighborsWithZeros(:,end),RaftNodes),:);
end

%Positions of All NeighborsWithZeros
% NeighborsPositions = ...
%     NodeIndex(ismember(NodeIndex(:,1),NeighborsWithZeros(:,end)),4:5);
%Check uniqueness of neighbors
if OnlyCheckRaft==1
    while size(RaftNeighborsWithZeros,1)>0 %While there is at least 1 raft node with zeros for neighbors
        TempNode = RaftNeighborsWithZeros(1,end);
        ExistingNeighb = (Neighbors(Neighbors(:,end)==TempNode,1:MaxNoNeighb));
        TempNeighborsWithZeros = NeighborsWithZeros;
        
        %Remove already-existing neighbors
        TempNeighborsWithZeros(ismember(TempNeighborsWithZeros(:,end),ExistingNeighb'),:)=[];
        
        %Remove the TempNode itself
        TempNeighborsWithZeros(TempNeighborsWithZeros(:,end)==TempNode,:) = [];
        
        if ~isempty(TempNeighborsWithZeros)
            TempNodeInfo = NodeIndex(NodeIndex(:,1)==TempNode,:);
            Pos = TempNodeInfo(1,4:5);
            
            %No. Missing Neighbors
            ZerosColumns = find(RaftNeighborsWithZeros(1,:)==0);
            NoNeeded = size(ZerosColumns,2);
            
            %Positions of other nodes excluding the TempNode
            OtherNodeInfo = NodeIndex(...
                ismember(NodeIndex(:,1),TempNeighborsWithZeros(:,end)),:);
            OtherNodes = OtherNodeInfo(:,1);
            OthersPositions = OtherNodeInfo(:,4:5);
            
            TempNodePosMatrix = ones(size(OthersPositions)).*Pos;
            
            %Calculate distance to each candidate
            Disp = OthersPositions-TempNodePosMatrix;
            Dist = vecnorm(Disp,2,2);
            
            %Crop the needed number of new neighbors from the closest pool
            SortMatrix = [Dist,OtherNodes];
            SortMatrix = sortrows(SortMatrix,1);
            
            if size(SortMatrix,1)>=NoNeeded
                SortMatrixCandidates = SortMatrix(1:NoNeeded,:);
            else
                SortMatrixCandidates = SortMatrix(1:end,:);
                NoNeeded = size(SortMatrix,1);
            end
            SortMatrixCandidates(SortMatrixCandidates(:,1)>NewNeighDisThresh,:) = [];
            NewDist = SortMatrixCandidates(:,1);
            NewNeighbors = SortMatrixCandidates(:,2);
            NoNewNeighbs = size(NewNeighbors,1);
            
            if size(SortMatrixCandidates,1)>=NoNeeded
                Neighbors(Neighbors(:,end)==TempNode,ZerosColumns)...
                    = NewNeighbors';
            else
                Neighbors(Neighbors(:,end)==TempNode,ZerosColumns(1:NoNewNeighbs))...
                    = NewNeighbors';
            end
            
            for i=1:size(NewNeighbors,1)
                %                 if NewDist(i,1)<=NewNeighDisThresh
                %Place new neighbor in the correct position of TempNode
                Neighbors(Neighbors(:,end)==TempNode,ZerosColumns(i)) = NewNeighbors(i);
                
                %Place TempNode in the correct position of NewNeighbor
                NNTempInfo = Neighbors(Neighbors(:,end)==NewNeighbors(i),:);
                ZerosOfNN = find(NNTempInfo==0);
                NoZerosOfNN = size(ZerosOfNN,2);
                Neighbors(Neighbors(:,end)==NewNeighbors(i),ZerosOfNN(1)) = TempNode;
                
                %If this was the only neighbor needed by NewNeighbor(i),
                %remove it from the pool, as well. It's Number of Zeros
                %will be automatically downgraded
                if NoZerosOfNN == 1
                    NeighborsWithZeros(NeighborsWithZeros(:,end)==NewNeighbors(i),:)=[];
                end
                %                 end
            end
        end
        RaftNeighborsWithZeros(RaftNeighborsWithZeros(:,end)==TempNode,:)=[];
        NeighborsWithZeros(NeighborsWithZeros(:,end)==TempNode,:)=[];
        NoRaftNeighborsWithZeros = size(RaftNeighborsWithZeros,1);
    end
else
    while size(NeighborsWithZeros,1)>1 %While there are at least 2 nodes with zeros for neighbors
        %     if ~isempty(NodeIndex(ismember(NodeIndex(:,1),NeighborsWithZeros(1,end)),4:5))
        TempNode = NeighborsWithZeros(1,end);
        ExistingNeighb = (Neighbors(Neighbors(:,end)==TempNode,1:MaxNoNeighb));
        TempNeighborsWithZeros = NeighborsWithZeros;
        
        %Remove already-existing neighbors
        TempNeighborsWithZeros(ismember(TempNeighborsWithZeros(:,end),ExistingNeighb'),:)=[];
        
        %Remove the TempNode itself
        TempNeighborsWithZeros(TempNeighborsWithZeros(:,end)==TempNode,:) = [];
        
        if ~isempty(TempNeighborsWithZeros)
            TempNodeInfo = NodeIndex(NodeIndex(:,1)==TempNode,:);
            Pos = TempNodeInfo(1,4:5);
            
            %No. Missing Neighbors
            ZerosColumns = find(NeighborsWithZeros(1,:)==0);
            NoNeeded = size(ZerosColumns,2);
            
            %Positions of other nodes excluding the TempNode
            OtherNodeInfo = NodeIndex(...
                ismember(NodeIndex(:,1),TempNeighborsWithZeros(:,end)),:);
            OtherNodes = OtherNodeInfo(:,1);
            OthersPositions = OtherNodeInfo(:,4:5);
            
            TempNodePosMatrix = ones(size(OthersPositions)).*Pos;
            
            %Calculate distance to each candidate
            Disp = OthersPositions-TempNodePosMatrix;
            Dist = vecnorm(Disp,2,2);
            
            %Crop the needed number of new neighbors from the closest pool
            SortMatrix = [Dist,OtherNodes];
            SortMatrix = sortrows(SortMatrix,1);
            
            if size(SortMatrix,1)>=NoNeeded
                SortMatrixCandidates = SortMatrix(1:NoNeeded,:);
            else
                SortMatrixCandidates = SortMatrix(1:end,:);
                NoNeeded = size(SortMatrix,1);
            end
            SortMatrixCandidates(SortMatrixCandidates(:,1)>NewNeighDisThresh,:) = [];
            NewDist = SortMatrixCandidates(:,1);
            NewNeighbors = SortMatrixCandidates(:,2);
            NoNewNeighbs = size(NewNeighbors,1);
            
            if size(SortMatrixCandidates,1)>=NoNeeded
                Neighbors(Neighbors(:,end)==TempNode,ZerosColumns)...
                    = NewNeighbors';
            else
                Neighbors(Neighbors(:,end)==TempNode,ZerosColumns(1:NoNewNeighbs))...
                    = NewNeighbors';
            end
            
            for i=1:size(NewNeighbors,1)
                %                 if NewDist(i,1)<=NewNeighDisThresh
                %Place new neighbor in the correct position of TempNode
                Neighbors(Neighbors(:,end)==TempNode,ZerosColumns(i)) = NewNeighbors(i);
                
                %Place TempNode in the correct position of NewNeighbor
                NNTempInfo = Neighbors(Neighbors(:,end)==NewNeighbors(i),:);
                ZerosOfNN = find(NNTempInfo==0);
                NoZerosOfNN = size(ZerosOfNN,2);
                Neighbors(Neighbors(:,end)==NewNeighbors(i),ZerosOfNN(1)) = TempNode;
                
                %If this was the only neighbor needed by NewNeighbor(i),
                %remove it from the pool, as well. It's Number of Zeros
                %will be automatically downgraded
                if NoZerosOfNN == 1
                    NeighborsWithZeros(NeighborsWithZeros(:,end)==NewNeighbors(i),:)=[];
                end
                %                 end
            end
        end
        NeighborsWithZeros(NeighborsWithZeros(:,end)==TempNode,:)=[];
        NoNeighborsWithZeros = size(NeighborsWithZeros,1);
    end
end

%Enforce uniqueness of neighbors
for i=1:size(Neighbors,1)
    TempNeighbors = Neighbors(i,1:(size(Neighbors,2)-1));
    Neighbors(i,1:(size(Neighbors,2)-1)) = zeros(size(TempNeighbors));
    TempNeighbors(TempNeighbors(:)==0)=[];
    UniqueTempNeighbors = unique(TempNeighbors);
    Neighbors(i,1:size(UniqueTempNeighbors,2)) = UniqueTempNeighbors;
end

DeleteOldNodesFromNeighbors;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function NormVel = CalculateNormalizedVelocity(Corners,t0,tf)

global SurfNodes SurfNodesPrev dt InputInterval it

Time = it*dt/60;
DeltaT = InputInterval*dt;
if t0<=Time && Time<=tf && it>InputInterval
    %Isolate particles in ROI
    
    BL = Corners(1,:); BR = Corners(2,:); 
    TR = Corners(3,:); TL = Corners(4,:);
    
    SurfNodesInROI_rows = inpolygon(SurfNodes(:,4),SurfNodes(:,5),...
        Corners(:,1),Corners(:,2)); %1 where SurfNode is inside ROI
    SurfNodesInROI0_rows = inpolygon(SurfNodesPrev(:,4),SurfNodesPrev(:,5),...
        Corners(:,1),Corners(:,2)); 
    
    NewSurfNodes = SurfNodes(SurfNodesInROI_rows,:);    %SurfNodes in ROI at t
    OldSurfNodes = SurfNodesPrev(SurfNodesInROI0_rows,:);   %SurfNodes in ROI at t-dt
    
    ConservedSurfNodeInfo = ...
        OldSurfNodes(ismember(OldSurfNodes(:,1),NewSurfNodes(:,1)),:);
    ConservedSurfNodes = unique(ConservedSurfNodeInfo(:,1),'rows');   
    
    [~,UniqueRows] = unique(NewSurfNodes(:,1),'rows');
    [~,UniqueRowsPrev] = unique(OldSurfNodes(:,1),'rows');
    
    UniqueSurfNodes = NewSurfNodes(UniqueRows,:);
    UniqueSurfNodesPrev = OldSurfNodes(UniqueRowsPrev,:);
    
    X = UniqueSurfNodes(ismember(UniqueSurfNodes(:,1),ConservedSurfNodes),4);
    X0 = ...
        UniqueSurfNodesPrev(ismember(UniqueSurfNodesPrev(:,1),ConservedSurfNodes),4);
    Y = UniqueSurfNodes(ismember(UniqueSurfNodes(:,1),ConservedSurfNodes),5);
    Y0 = ...
        UniqueSurfNodesPrev(ismember(UniqueSurfNodesPrev(:,1),ConservedSurfNodes),5);
    
    Disp = [X,Y]-[X0,Y0];
    Velocities = Disp/DeltaT;
    Speeds = vecnorm(Velocities,2,2);
    
    NoNodesWithHist = size(X,1);
    SumVel = sum(Velocities,1);
    MeanSpeed = mean(Speeds);
    
    NormVel(1) = Time;
    NormVel(2:3) = 1/(NoNodesWithHist*MeanSpeed)*SumVel;
    NormVel(4) = norm(NormVel(2:3));
    if NoNodesWithHist<=1
        NormVel = [Time 0 0 0];
    end
else
    NormVel = [Time 0 0 NaN];
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CalculateOutputs
  
global RaftNodes SurfNodes EdgeNodes ct parkrate unparkrate...
    Nr Ns Ntot phi Area Circumference SurfExcess...
    gamma alpha delta tau

%Geometric Outputs
Nr(ct,1) = size(RaftNodes,1);
Ns(ct,1) = size(SurfNodes,1);
Ntot(ct,1) = Ns(ct,1) + Nr(ct,1);

phi(ct,1) = Ns(ct,1)/Nr(ct,1);

Area(ct,1) = Nr(ct,1);
Circumference(ct,1) = size(EdgeNodes,1);

SurfExcess(ct,1) = Circumference(ct,1)/(2*sqrt(Area(ct,1)*pi));

%Growth Kinetics
if ct>1
    gamma(ct,1) = Area(ct,1)*parkrate(ct,1)/Circumference(ct,1);
    alpha(ct,1) = parkrate(ct,1);
    delta(ct,1) = unparkrate(ct,1);
end

%Calculate Global Decay rate
if ct>1
    tau(ct,1) = CalculateDecayRate;
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tau = CalculateDecayRate

global NodeIndex NodeIndexPrev dt InputInterval

%Identify Raft Nodes 

NodeIndex_temp = NodeIndex(NodeIndex(:,2)==1,:);
PrevNodeIndex_temp = NodeIndexPrev(NodeIndexPrev(:,2)==1,:);

CommonNodes = intersect(NodeIndex_temp(:,1),PrevNodeIndex_temp(:,1));

[~,NewNodeRows] = ismember(CommonNodes,NodeIndex_temp(:,1));
[~,OldNodeRows] = ismember(CommonNodes,PrevNodeIndex_temp(:,1));

NewPos = NodeIndex_temp(NewNodeRows,4:5);
OldPos = PrevNodeIndex_temp(OldNodeRows,4:5);
Disp = NewPos-OldPos;
Dist = vecnorm(Disp,2,2);

DecayFraction = vecnorm(NewPos,2,2)./vecnorm(OldPos,2,2);

AvgDecayFraction = mean(DecayFraction);

tau = -(dt*InputInterval)/log(AvgDecayFraction)/60; %min

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PlotColoredSurfTraffic

global NodeIndex RaftNodes RedExp SurfNodes...
    mmPerNode LengthScale NWindow ParticleSize...
    NewNeighDisThresh InputInterval dt ct

FontSize = 12;
LineWidth = 2;
Outline = 1;

%ID Which nodes to plot
PercentPlotted = 100;

UptoDateNodes = SurfNodesArray{:,ct};
MaxNodeNo = max(UptoDateNodes(:,1));
ShowNodes = (1:ceil(100/PercentPlotted):MaxNodeNo)';

TimeRange = 13; %s
NoConcurrentFrames = ceil(TimeRange/(dt*InputInterval));

%Define Marker Size

BckGrndMarkerSize = 1000*1/NWindow^RedExp;

xs = NodeIndex(:,4); ys = NodeIndex(:,5);

nFree = size(SurfNodes,1);
nRaft = size(RaftNodes,1);
nAdd = nFree+nRaft; %total added plus original raft size
phi_t = nFree/nRaft;

%Plot background lattice
scatter(xs,ys,BckGrndMarkerSize,'k.')
hold on
daspect([1 1 1])
xlim([-NWindow/2 NWindow/2])
ylim([-NWindow/2 NWindow/2])

% Get plot size in pixel units
pos = get(gcf, 'Position'); %// gives x left, y bottom, width, height
FigWidth = pos(3);  %figure width in pixels
MarkerDiam = ParticleSize*FigWidth/NWindow;
MarkerSize = floor((pi/4*MarkerDiam)^2); %in num pixels
if MarkerSize==0
    MarkerSize=1;
end
clf

%Plot a length scale bar
if LengthScale==1
    LengthBarFract = 5;
    LengthBar = NWindow/LengthBarFract;
    LengthBar_mm = round(LengthBar*mmPerNode,-1);
    LengthBar = round(LengthBar_mm/mmPerNode);
    StartBar_x = -LengthBar/2;
    EndBar_x = StartBar_x+LengthBar;
    StartBar_y = -0.9*NWindow/2;
    EndBar_y = StartBar_y;
    
    PlotBarX = [StartBar_x EndBar_x];
    PlotBarY = [StartBar_y EndBar_y];
    
    plot(PlotBarX,PlotBarY,'k','LineWidth',LineWidth)
    hold on
    
    LengthBarText = [num2str(LengthBar_mm,'%.0f'),' mm'];
    text(0,-0.85*NWindow/2,LengthBarText,'FontSize',FontSize,...
        'HorizontalAlignment','center',...
        'Interpreter','Latex')
end

%Plot Surface History
if ct>=NoConcurrentFrames
    StartFrame = ct-NoConcurrentFrames+1;
    for f=StartFrame:ct
        TempFrameNo = f-StartFrame+1;
        Transparency = TempFrameNo/NoConcurrentFrames;
        
        %Plot Surf Particles With Decreasing Transparency
        SurfNodes_temp = SurfNodesArray{:,f};
        SurfNodes_temp = SurfNodes_temp(ismember(SurfNodes_temp(:,1),ShowNodes),:);
        
        SurfPos = [SurfNodes_temp(:,4) SurfNodes_temp(:,5)];

        SurfNodeHist = SurfNodesHistoryArray{:,f};
        SurfNodeHist = SurfNodeHist(ismember(SurfNodeHist(:,1),ShowNodes),:);

        % Plot previous surface particles
        OldSurfPos = SurfNodeHist(:,5:6);
        NewSurfPos = SurfNodeHist(:,2:3);

        Disp = NewSurfPos-OldSurfPos;
        Dist = vecnorm(Disp,2,2);
        
        PlotDist = Dist;
        PlotDist(Dist>=0.85*NewNeighDisThresh,:) = 0.85*NewNeighDisThresh;
        
        UnitVec = Disp./Dist;
        
        AbdomenPos = NewSurfPos-UnitVec*1.*PlotDist;
        ThoraxPos = NewSurfPos-UnitVec*0.5.*PlotDist;
        AbdomenPos(AbdomenPos(:,1)==0,:)=[];
        ThoraxPos(ThoraxPos(:,1)==0,:)=[];
        
%         AbdomenColors = ones(size(AbdomenPos(:,1)))*TempFrameNo;
%         ThoraxColors = ones(size(ThoraxPos(:,1)))*TempFrameNo;
%         Colors = ones(size(SurfPos(:,1)))*TempFrameNo;        
%         
%         h1 = scatter3(AbdomenPos(:,1),AbdomenPos(:,2),AbdomenColors,...
%             MarkerSize,AbdomenColors,'filled');
%         hold on
%         h2 = scatter3(ThoraxPos(:,1),ThoraxPos(:,2),ThoraxColors,...
%             MarkerSize/2,ThoraxColors,'filled');
%         hold on
%         h = scatter3(SurfPos(:,1),SurfPos(:,2),Colors,...
%             MarkerSize,Colors,'filled');
%         hold on
        
        h1 = scatter(AbdomenPos(:,1),AbdomenPos(:,2),...
            MarkerSize,'b','filled','MarkerFaceAlpha',Transparency);
        hold on
        h2 = scatter(ThoraxPos(:,1),ThoraxPos(:,2),...
            MarkerSize/2,'b','filled','MarkerFaceAlpha',Transparency);
        hold on
        h = scatter(SurfPos(:,1),SurfPos(:,2),...
            MarkerSize,'b','filled','MarkerFaceAlpha',Transparency);
        hold on
        
%         colormap(jet)
        
        if Outline==1
            set(h1,'MarkerEdgeColor', 'k');
            set(h2,'MarkerEdgeColor', 'k');
            set(h,'MarkerEdgeColor', 'k');
            hold on
        end
    end    
else %Plot all available SurfNodeInfo
    for f=1:ct
        TempFrameNo = f;
        Transparency = TempFrameNo/NoConcurrentFrames;
        
        %Plot Surf Particles With Decreasing Transparency
        SurfNodes_temp = SurfNodesArray{:,f};
        SurfNodes_temp = SurfNodes_temp(ismember(SurfNodes_temp(:,1),ShowNodes),:);
        
        SurfPos = [SurfNodes_temp(:,4) SurfNodes_temp(:,5)];

        SurfNodeHist = SurfNodesHistoryArray{:,f};
        SurfNodeHist = SurfNodeHist(ismember(SurfNodeHist(:,1),ShowNodes),:);

        % Plot previous surface particles
        OldSurfPos = SurfNodeHist(:,5:6);
        NewSurfPos = SurfNodeHist(:,2:3);

        Disp = NewSurfPos-OldSurfPos;
        Dist = vecnorm(Disp,2,2);
        
        PlotDist = Dist;
        PlotDist(Dist>=0.85*NewNeighDisThresh,:) = 0.85*NewNeighDisThresh;
        
        UnitVec = Disp./Dist;
        
        AbdomenPos = NewSurfPos-UnitVec*1.*PlotDist;
        ThoraxPos = NewSurfPos-UnitVec*0.5.*PlotDist;
        AbdomenPos(AbdomenPos(:,1)==0,:)=[];
        ThoraxPos(ThoraxPos(:,1)==0,:)=[];
        AbdomenColors = ones(size(AbdomenPos(:,1)))*TempFrameNo;
        ThoraxColors = ones(size(ThoraxPos(:,1)))*TempFrameNo;
        Colors = ones(size(SurfPos(:,1)))*TempFrameNo;
        
%         h1 = scatter3(AbdomenPos(:,1),AbdomenPos(:,2),AbdomenColors,...
%             MarkerSize,AbdomenColors,'filled');
%         hold on
%         h2 = scatter3(ThoraxPos(:,1),ThoraxPos(:,2),ThoraxColors,...
%             MarkerSize/2,ThoraxColors,'filled');
%         hold on
%         h = scatter3(SurfPos(:,1),SurfPos(:,2),Colors,...
%             MarkerSize,Colors,'filled');
%         hold on
 
        h1 = scatter(AbdomenPos(:,1),AbdomenPos(:,2),...
            MarkerSize,'b','filled','MarkerFaceAlpha',Transparency);
        hold on
        h2 = scatter(ThoraxPos(:,1),ThoraxPos(:,2),...
            MarkerSize/2,'b','filled','MarkerFaceAlpha',Transparency);
        hold on
        h = scatter(SurfPos(:,1),SurfPos(:,2),...
            MarkerSize,'b','filled','MarkerFaceAlpha',Transparency);
        hold on
        
%         colormap(jet)
        
        if Outline==1
            set(h1,'MarkerEdgeColor', 'k');
            set(h2,'MarkerEdgeColor', 'k');
            set(h,'MarkerEdgeColor', 'k');
            hold on
        end
    end
end

% caxis([0 NoConcurrentFrames])

PlotFigureText1(phi_t,nRaft,nFree,nAdd,FontSize);

set(gcf,'color','w');
daspect([1 1 1])
xlim([-NWindow/2 NWindow/2])
ylim([-NWindow/2 NWindow/2])
xticks([])
yticks([])
box on

end
