function [CurrentLattice,Ctots,Dtots,ClusterData] = MixDetNowakLiaoLoop(InitialLattice,Vt,Vr,T,R,P,S,base,N,FilePrefix)
%   ================================================
%   David Liao
%   2008 January 27
%   Ising Trio Loop
%   ================================================
%
%   Requires NowakLiao.m and IsingTrioFluorescence.m
%   [CurrentLattice,Ctots,Dtots,ClusterData] =
%   NowakLiaoLoop(InitialLattice,Vt,Vr,T,R,P,S,base,N,FilePrefix)

ClusterData=zeros(N,2);

L = size(InitialLattice);
NeighborFilter = zeros(L);
NeighborFilter(2,1,1) = 1; NeighborFilter(1,2,1) = 1; NeighborFilter(1,1,2) = 1;
NeighborFilter(L(1),1,1) = 1; NeighborFilter(1,L(2),1) = 1; NeighborFilter(1,1,L(3)) = 1;
NeighborFilterfft = fftn(NeighborFilter);

NeighborFilter(1,1,1) = 1;
AutoNeighborFilterfft = fftn(NeighborFilter);

ZerothColorImage = IsingTrioFluorescence(InitialLattice);
imwrite(ZerothColorImage,[FilePrefix '000000.png']);
CurrentLattice = InitialLattice;

Ctots = zeros(N,1);
Dtots = zeros(N,1);

tic;

for i = 1:N
    [CurrentLattice,CurrentCtot,CurrentDtot] = MixDetNowakLiao(CurrentLattice,NeighborFilterfft,AutoNeighborFilterfft,Vt,Vr,T,R,P,S,base);
    CurrentColorImage = IsingTrioFluorescence(CurrentLattice);
    
    [a_r,b_r]=CorrelationRadius2(CurrentColorImage(:,:,1),0.8,5);
    ClusterData(i,1)=a_r;
      
    [a_g,b_g]=CorrelationRadius2(CurrentColorImage(:,:,2),0.8,5);
    ClusterData(i,2)=a_g;    
    
    imwrite(CurrentColorImage,[FilePrefix num2str(i,'%06g') '.png']);
    
    Ctots(i) = CurrentCtot;
    Dtots(i) = CurrentDtot;
    
    ElapsedTime = toc;
    
    RemainingTime = (N - i)*ElapsedTime/i;
    fprintf(['Remaining time: ' ,num2str(RemainingTime), ' s.  Elapsed time: ' ,num2str(ElapsedTime), ' s.  \n']);
end