clear all;
close all;

N=432;

Pd={'0.01','0.02','0.03'};
cM={'sr','xb','og','vm','*c'};
%cML={'s:r','x:b','o:g','v:m','*:c'};
cL={'-r','-b','-g','-m','-c'};

crystal = 'HS1';

topDir = './jammingHS1/';

allSpectraFile = 'allSpectraHessianGeometric.dat';

nBinsLog = 10;
nBinsLin = 100;
for i =1:numel(Pd)

    dataDir = [ topDir 'N' num2str(N) '_s' Pd{i} filesep ];
    allSpectraFilename = [ dataDir allSpectraFile ];

    S=[];
    Bs=[];
    Hs=[];
    Hsum=[];
    Xs=[];
    
    S = load(allSpectraFilename);

    S = S(S>0);
    S = sqrt(S);
    S = S(S>1.0e-7);

    Smin = min(S(:));
    Sminexp = floor(log10(Smin))-1;
    Smax = max(S(:));
    Smaxexp = round(log10(Smax))+1;
    So = 5.0e-1;
    Soexp = log10(So);
    Blog = logspace(Sminexp,Soexp,nBinsLog+1);
    Blin = linspace(So,Smax,nBinsLin+1); 
    Blog(end) = [];
    Bs = [ Blog Blin];
    
    Hs = histcounts(S, Bs, 'Normalization','pdf');
    Bs(end) = [];

    idx = [];
    idx = (Hs<2e-1);
    idx(10:end)=0;
    idx;
    Hs(idx) = [];
    Bs(idx) = [];
    
    loglog( Bs, Hs, cL{i}, 'markers', 4, 'DisplayName', ['Pd=' Pd{i}]);
    hold on;

end

xlim([1e-4,5e+0])
ylim([1e-4,2e+0])
%axis equal
xlabel('\omega');
ylabel('PDF(\omega)');
legend('-DynamicLegend', 'Location', 'Best');
title([ crystal '  All Spectra Hessian Geometric' '   N = ' num2str(N) ]); 




axes('Position',[.42 .37 .30 .30]);
box on;

meanORmedian = 0; %0 -> mean , 1 -> median 

allIPRFile = 'allIPRGeometric.dat';

%nBinsX = 30;
nBinsLog = 10;
nBinsLin = 100;
for i =1:numel(Pd)

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

    S=[];
    X=[];
    Y=[];
    Xout=[];
    xBins=[];
    xCounts=[];
    idx=[];
    xCenter=[];
    yMean=[];
    yMedian=[];
    
    S = load(allIPRFilename);
    X = S(:,1);
    Y = S(:,2);

    Xout = (X>1e-5);
    X = X(Xout);
    Y = Y(Xout);

    %First, trim out any zeros
    zeros = (X == 0) | (Y == 0);
    X = X(~zeros);
    Y = Y(~zeros);

    Xmin = min(X(:));
    %Xminexp = floor(log10(Xmin))-1;
    Xminexp = log10(Xmin) - 0.01;
    Xmax = max(X(:));
    %Xmaxexp = round(log10(Xmax))+1;
    Xmaxexp = log10(Xmax) + 0.01;

    Xo = 5.0e-1;
    Xoexp = log10(Xo);
    xBinslog = logspace(Xminexp,Xoexp,nBinsLog+1);
    xBinslin = linspace(Xo,Xmax,nBinsLin+1); 
    xBinslog(end) = [];
    xBins = [ xBinslog xBinslin];

    [xCounts,~,idx] = histcounts(X, xBins);
    xCountSum = sum(xCounts(:));
    xBinWid = diff(xBins);
    xCounts = xCounts./xBinWid;
    xBins(end) = [];
    xCounts = xCounts./xCountSum;
    xCenter = xBinWid/2 + xBins;

    if meanORmedian == 0
        yMean = accumarray(idx(:),Y,[],@mean);
        loglog( xCenter, yMean, cL{i});
            
    else
        yMedian = accumarray(idx(:),Y,[],@median);
        loglog( xCenter, yMedian, cL{i}, 'DisplayName', ['Pd=' Pd{i}]);
    end
    
    hold on;

end

xlim([1e-4,5e+0]);
ylim([1e-3,1e+0]);
xticks([1e-4,1e-3,1e-2,1e-1,1e+0]);
yticks([1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1e+0]);
%axis equal
xlabel('\omega');
ylabel('Y(\omega)');
%legend('-DynamicLegend', 'Location', 'Best');