function [CoopAreas,dCoopAreas,DefectAreas,dDefectAreas,CNL,dCNL,DNL,dDNL,basevector] = NowakLiaoEnsemble(L,InitializationRatio,Vt,Vr,T,R,P,S,basevector,Ntime,Nrandoms,FilePrefix)
%   =============================
%   David Liao
%   2008 February 25
%   Ensemble calculation on Nowak-Liao model
%   =============================
%
%   [CoopAreas,dCoopAreas,DefectAreas,dDefectAreas,CNL,dCNL,DNL,dDNL,baseve
%   ctor] =
%   NowakLiaoEnsemble(L,InitializationRatio,Vt,Vr,T,R,P,S,basevector,Ntime,
%   Nrandoms,FilePrefix)

CoopAreas = []; dCoopAreas = []; DefectAreas = []; dDefectAreas = [];
CNL = []; dCNL = []; DNL = []; dDNL = [];

Nparams = length(basevector);

N = Nparams*Nrandoms;

tic;

for j = 1:Nparams
    
    baseNow = basevector(j);
    
    CoopAreasSingleParam = [];
    DefectAreasSingleParam = [];
    CNLSingleParam = [];
    DNLSingleParam = [];
    
    for i = 1:Nrandoms
        
        RandNow = rand(L);
        InitialLatticeNow = (RandNow > InitializationRatio) - (RandNow <= InitializationRatio);
        [NowLattice,NowDetCNL,NowDetDNL,NowDetClusters] = DetNowakLiaoLoop(InitialLatticeNow,Vt,Vr,T,R,P,S,baseNow,Ntime,[FilePrefix num2str(i,'%06g') 'Rand' num2str(j,'%06g') 'Param']);
        
        CoopAreasSingleParam = [CoopAreasSingleParam NowDetClusters(:,2)];
        DefectAreasSingleParam = [DefectAreasSingleParam NowDetClusters(:,1)];
        
        CNLSingleParam = [CNLSingleParam NowDetCNL];
        DNLSingleParam = [DNLSingleParam NowDetDNL];
        
        Nnow = (j-1)*Nrandoms + i;
        
        fprintf(['===== ELAPSED OPERATIONS: ' ,num2str(Nnow), ' TOTAL OPERATIONS: ' ,num2str(N), ' =====\n']);

    end
    
    CoopAreas = [CoopAreas mean(CoopAreasSingleParam,2)];
    dCoopAreas = [dCoopAreas std(CoopAreasSingleParam,0,2)/sqrt(Nrandoms)];
    
    DefectAreas = [DefectAreas mean(DefectAreasSingleParam,2)];
    dDefectAreas = [dDefectAreas std(DefectAreasSingleParam,0,2)/sqrt(Nrandoms)];
    
    CNL = [CNL mean(CNLSingleParam,2)];
    dCNL = [dCNL std(CNLSingleParam,0,2)/sqrt(Nrandoms)];
    
    DNL = [DNL mean(DNLSingleParam,2)];
    dDNL = [dDNL std(DNLSingleParam,0,2)/sqrt(Nrandoms)];
    
end