clear all:
close all;

nParticles=432;
Pd={'0.01','0.02','0.03'};

cM={'sr','xb','og','vm','*c'};

localizedFlag=0; % 0=extended, 1=localized

dim1=500;
dim2=400;
fig=figure('position',[0,0,dim1,dim2]);

lattice = 'HS1';

topDir = './jammingHS1/';
allExtendedForcesFile = 'allModeExtendedForces.dat';
allLocalizedForcesFile = 'allModeLocalizedForces.dat';

for i =1:numel(Pd)

    dataDir = [ topDir 'N' num2str(nParticles) '_s' strjoin(Pd(i)) filesep ];

    allExtendedForcesFilename = [ dataDir allExtendedForcesFile ];
    allLocalizedForcesFilename = [ dataDir allLocalizedForcesFile ];

    allExtendedForces = load(allExtendedForcesFilename);
    allLocalizedForces = load(allLocalizedForcesFilename);

    fext = allExtendedForces(:);
    fext = fext(fext > 1e-16 );
    floc = allLocalizedForces(:);
    floc = floc(floc > 1e-16 );
    fall = [fext; floc];

    favg = mean(fall(:));
    floc = floc/favg;
    fext = fext/favg;
    
    Mext = numel(fext);
    Gext = sort(fext);
    xext = (1:Mext)/Mext;

    floc = floc(floc > 1e-6 );
    Mloc = numel(floc);
    Gloc = sort(floc);
    xloc = (1:Mloc)/Mloc;

    if localizedFlag == 0
        loglog( Gext, xext, strjoin(cM(i)), 'DisplayName', ['s=' Pd{i}]);
    elseif localizedFlag == 1
        loglog( Gloc, xloc, strjoin(cM(i)), 'DisplayName', ['s=' Pd{i}]);
    end
        
    hold on;
    %loglog( G, x, 'sb');

end

if localizedFlag == 0
    a1 = 1.42311;
    b1 = 5e+0;
    G = Gext;
elseif localizedFlag == 1
    a1 = 1.17462;
    b1 = 8e+1;
    G = Gloc;
end

PL1 = b1*(G.^a1);
loglog( G, PL1, 'k-', 'DisplayName', ['slope =' num2str(a1)],'LineWidth',2);

%axis equal
if localizedFlag == 0
    xlim([1e-6,1e+2])
    ylim([1e-6,2e+0])
elseif localizedFlag == 1
    xlim([1e-6,1e+2])
    ylim([1e-6,2e+0])
end
xlabel('f/<f>')
ylabel('CDF(f/<f>)')
legend('-DynamicLegend', 'Location', 'Best');

if localizedFlag == 0
    title([ lattice '  all extended forces' '   N = ' num2str(nParticles) ]); 
elseif localizedFlag == 1
    title([ lattice '  all localized forces' '   N = ' num2str(nParticles) ]);
end