function [] = phiJCompareSSHSPaper()

currentFile = mfilename( 'fullpath' );
[pathstr,~,~] = fileparts( currentFile );
GMfilename = [pathstr,filesep,'..',filesep,'Data',filesep,'softSpheres',filesep,'GM',filesep,'allPhiJ.dat'];

% GMfilename = '/Users/pkmorse/Documents/Data/Duke/softSpheres/GM/allPhiJ.dat';
% IQfolder = '/Users/pkmorse/Documents/Data/Duke/IQFinSizeSweep';
data = readtable(GMfilename);
dimInput = data{:,1};
nInput = data{:,2};
seedInput = data{:,3};
phiJInput = data{:,4};

for dim = 3:6
    figure(dim*1000)
    clf; hold all
    figure(dim*100)
    clf; hold all
    if dim == 3
        nP = [64 128 256 512 1024 2048 4096 8192 16384];
        fitRange = 2:numel(nP)-1;
        fitRangeGM = 2:numel(nP)-1;
        JonM = [0.634306591259341, 0.638132148277319, 0.640220021317927, 0.641964239626639, 0.643124028148778, 0.644182727124316, 0.645099187297707, 0.645759619073016, NaN];
        JonS = 1e-3*[0.219631842474421, 0.175486210451159, 0.158688655186312, 0.146616905418173, 0.136749591665652, 0.128890586645558, 0.124434457985934, 0.120091970657432, NaN];
        JonInf = 0.6487;
    elseif dim == 4
        nP = [128 256 512 1024 2048 4096 8192 16384 32768];
        fitRange = 3:numel(nP)-2;
        fitRangeGM = 3:numel(nP)-2;
        JonM = [0.440099186341750, 0.443627633773675, 0.446187133759240, 0.447648927417493, 0.448680155821864, 0.449672505114451, 0.450368243669639, NaN, NaN];
        JonS = 1e-4*[1.62354993699543, 1.10766159769682, 0.910447250798696, 0.869299113712696, 0.840926081286386, 0.811263294352993, 0.810657297763417, NaN, NaN];
        JonM(1) = NaN;
        JonS(1) = NaN;
        JonInf = 0.4564;
    elseif dim == 5
        nP = [256 512 1024 2048 4096 8192 16384 32768];
        fitRange = 1:numel(nP)-3;
        fitRangeGM = 1:numel(nP)-3;
        JonM = [0.297827512843182, 0.298987578158429, 0.300186811239760, 0.301439933944474, 0.302208074420703, NaN, NaN, NaN];
        JonS = 1e-3*[0.148266230764804, 0.136631676852739, 0.132573308438855, 0.127111833578664, 0.124135587768526, NaN, NaN, NaN];
        JonInf = 0.3083;
    else %dim = 6
        nP = [256 512 1024 2048 4096 8192 16384 32768];
        fitRange = 2:numel(nP)-3;
        fitRangeGM = 2:numel(nP)-2;
        JonM = [NaN, 0.193306207445564, 0.194096770840507, 0.194970184113463, 0.195443397939522,NaN, NaN, NaN];
        JonS = 1e-4*[NaN, 0.904829836495813, 0.924562802635254, 0.829082747740295, 0.899852041804372,NaN, NaN, NaN];
        JonInf = 0.2008;
    end
    
    JGMM = NaN(numel(nP),1);
    JGMS = NaN(numel(nP),1);
    JIQM = NaN(numel(nP),1);
    JIQS = NaN(numel(nP),1);
    for i = 1:numel(nP)
        try
            JList = phiJInput(dimInput==dim&nInput==nP(i));
            JGMM(i) = mean(JList);
            JGMS(i) = std(JList)/sqrt(numel(JList));
%             thisIQfolder = [IQfolder,filesep,num2str(dim),filesep,num2str(nP(i)),filesep];
%             numAvgs = numel(dir(thisIQfolder)) - 3;
%             phiDirList = dir([thisIQfolder,filesep,num2str(1),filesep,'0.*']);
%             numPhi = numel(phiDirList);
%             phiList = zeros(numPhi,1);
%             for j = 1:numPhi
%                 phiList(j) = str2double(phiDirList(j).name);
%             end
%             isJammedFilename = [thisIQfolder,filesep,num2str(j),filesep,'JList.dat'];
%             if exist(isJammedFilename,'file')
%                 isJammed = load(isJammedFilename);
%             else
%                 isJammed = zeros(numPhi,1);
%                 for j = 0:numAvgs-1
%                     for k = 1:numPhi
%                         isJammed(k) = isJammed(k) + load([thisIQfolder,filesep,num2str(j),filesep,phiDirList(k).name,filesep,'isJammed.dat']);
%                     end
%                 end
%                 isJammed = isJammed/numAvgs;
%                 saveascii(isJammed,isJammedFilename)
%             end
%             [res,err] = erfPercFit(phiList,isJammed);
%             JIQM(i) = res.phiLR;
%             JIQS(i) = res.width;
        catch
        end
    end
    figure(dim*1000)
%     errorbar(nP.^(-1/dim),JIQM,JIQS,'MarkerSize',10,'Linewidth',1.5,'Color','b')
    errorbar(nP(fitRangeGM).^(-1/dim),JGMM(fitRangeGM),JGMS(fitRangeGM),'.','MarkerSize',10,'Linewidth',1.5,'Color','k')
    errorbar(nP(fitRange).^(-1/dim),JonM(fitRange),JonS(fitRange),'.','MarkerSize',10,'Linewidth',1.5,'Color','r')
    x50 = linspace(min(nP(fitRange))^(-1/dim),max(nP(fitRange))^(-1/dim),100);
    [res,err] = fit(nP(fitRange)'.^(-1/dim),JonM(fitRange)','poly1');
    plot(x50',res.p1*x50'+res.p2,'-','Linewidth',1.5,'Color','r')
    x50 = linspace(min(nP(fitRangeGM))^(-1/dim),max(nP(fitRangeGM))^(-1/dim),100);
    [res,err] = fit(nP(fitRangeGM)'.^(-1/dim),JGMM(fitRangeGM),'poly1');
    plot(x50',res.p1*x50'+res.p2,'-','Linewidth',1.5,'Color','k')
    figure(dim*100)
    errorbar(nP.^(-1/dim),JonInf - JonM,JonS,'MarkerSize',10,'Linewidth',1.5,'Color','r')
    figure(dim*1000)
    xlabel('N^{-1/d}')
    ylabel('\phi_{J0}')
    plotSpecs()
    set(gca,'Linewidth',2)
    set(gca,'FontSize',20)
    if dim == 3
        legend('SS algorithm','HS algorithm','location','northeast')
        axis([0 0.36 0.636 0.647])
    elseif dim == 4
        axis([0 0.36 0.444 0.452])
    elseif dim == 5
        axis([0 0.36 0.297 0.3045])
    else %dim == 6
        axis([0 0.36 0.193 0.198])
    end
       
end