function [] = onsetPlotsLogLSQPaper()
currentFile = mfilename( 'fullpath' );
[pathstr,~,~] = fileparts( currentFile );
filename = [pathstr,filesep,'..',filesep,'Data',filesep,'phiJ.dat'];

clearFigs = 1;
recalculate = 1;
phiOnInf = zeros(6,1);
phiJInf = zeros(6,1);
phiOnInfErr = zeros(6,1);
phiJInfErr = zeros(6,1);
phid = [NaN
    NaN
    0.57700000000000
   0.403600000000000
   0.265
   0.172300000000000];
if clearFigs == 1
    figure(333333)
    clf; hold all
    figure(333334)
    clf; hold all
    figure(1113)
    clf; hold all
    figure(52)
    clf; hold all
    for dim = [3 4 5 6]
        figure(dim * 10);
        clf; hold all;
    end
end

if recalculate
    data = readtable(filename);
    dimInput = data{:,1};
    nInput = data{:,2};
    phiEqInput = data{:,3};
    seedInput = data{:,4};
    zeroFInput = data{:,5};
    iterationFactorInput = data{:,6};
    alphaInput = data{:,7};
    phiJInput = data{:,8};
end
nu = 1;

loList = [0 0 0.4 0.2 0.1 0.05 0 0];
hiJList = [0 0 1 0.47 0.314 1 1 1];
colorList = [1 0 0; 0.75 0 0.25; 0.5 0 0.5; 0.25 0 0.75; 0 0 0; 0 0 1];
JList = [0 0 0.64811 0.4556 0.3092 0.2008 0 0.08];
JTrial = [0 0 0.64 0.4567 0.302 0.196677 0 0.08];
onList = [NaN NaN 0.52 0.335 0.205 0.12 NaN 0.039];
onTrial = [NaN NaN 0.55 0.38 0.24 0.13 NaN 0.039];
JMin = [0 0 0.6 0 0 0 0 0];
aGuess = [0 0 30 30 30 30 0 30];
dGuess = [0 0 0.02 0.02 0.02 0.02 0 0.02];

for dim = [3 4 5 6]
    figure(dim * 10)
    if dim == 3
        nP = [64 128 256 512 1024 2048 4096 8192]';
        fitRange = 1:numel(nP);
    elseif dim == 4
        nP = [128 256 512 1024 2048 4096 8192]';
        fitRange = 2:numel(nP);
    elseif dim == 5
        nP = [256 512 1024 2048 4096]';
        fitRange = 1:numel(nP);
    elseif dim == 6
        nP = [512 1024 2048 4096]';
        fitRange = 1:numel(nP);
    else
        nParticles = [512 1024 2048];
        nP = [512 1024 2048 8192 16384]';
    end
    lo = loList(dim);
    dimColor = [0 0 0; 0 0 0; 1 0 0; 0.66 0 0.33; 0.33 0 0.66; 0 0 1];
    
    xcell = cell(1,numel(nP));
    ycell = cell(1,numel(nP));
    wcell = cell(1,numel(nP));
    widthcell = cell(1,numel(nP));
    for j = 1:numel(nP)
        try
            if recalculate
                phi = phiEqInput(dimInput==dim&nInput==nP(j)&phiJInput>JMin(dim));
                phiJ = phiJInput(dimInput==dim&nInput==nP(j)&phiJInput>JMin(dim));
                thisHi = 1;
                % Avoiding the high density low N packings which
                % spontaneously crystalize
                if dim == 4
                    if nP(j) < 250
                        thisHi = 0.36;
                    elseif nP(j) < 500
                        thisHi = 0.4;
                    end
                end
                worthwhile = phi>=lo & phi < thisHi & phiJ < hiJList(dim);
                phiJ = phiJ(worthwhile);
                phi = phi(worthwhile);
                phiEq = unique(phi);
                phiM = NaN(size(phiEq)); phiS = NaN(size(phiEq)); phiW = NaN(size(phiEq));
                for i = 1:numel(phiEq)
                    phiM(i) = mean(phiJ(phi==phiEq(i)));
                    phiS(i) = std(phiJ(phi==phiEq(i)))/sqrt(sum(phi==phiEq(i)));
                    phiW(i) = std(phiJ(phi==phiEq(i)));
                end
                [~,~] = mkdir(['phiJHolder',filesep,num2str(dim),filesep,num2str(nP(j))]);
                saveascii([phiEq phiM phiS],['phiJHolder',filesep,num2str(dim),filesep,num2str(nP(j)),filesep,'phiData.dat']);
            else
                phiData = load(['phiJHolder',filesep,num2str(dim),filesep,num2str(nP(j)),filesep,'phiData.dat']);
                phiEq = phiData(:,1);
                phiM = phiData(:,2);
                phiS = phiData(:,3);
            end
            
            xcell{j} = phiEq;
            ycell{j} = phiM;
            wcell{j} = 1./sqrt(phiS);
%             wcell{j}(isinf(wcell{j})) = 0;
            widthcell{j} = phiW;
            figure(dim * 10)
            h = errorbar(phiEq,phiM,phiS,'.','Color',[(j-1)/(numel(nP)-1) 0 1 - (j-1)/(numel(nP)-1)],'Linewidth',1.5);
        catch
        end
    end
    y = linspace(0,1.53866667*dim*2^(-dim),1000);
    
    beta0 = zeros(3*numel(nP)+1,1);
    beta0(1) = 0.5;
    if dim == 6
        beta0(1) = 0.2;
    end
    beta0(2:(numel(nP)+1)) = 0.02;
    beta0(numel(nP)+2:2*numel(nP)+1) = onTrial(dim);
    beta0(2*numel(nP)+2:end) = JTrial(dim);
    mdl_cell1param = cell(1,numel(nP));
%     mdl_cell1paramH = cell(1,numel(nP));
    for k = 1:numel(nP)
        mdl_cell1param{k} = str2func(['@(beta,x) beta(1) * beta(',num2str(1+k),') * log(1+exp((x-beta(',num2str(numel(nP)+1+k),'))/beta(',num2str(1+k),'))) + beta(',num2str(2*numel(nP)+1+k),')']);
    end
    
    try
        y = linspace(0,phid(dim),1000);
        [beta1param,r,J,Sigma,~,~,~] = nlinmultifit(xcell(1:numel(nP)), ycell(1:numel(nP)), mdl_cell1param(1:numel(nP)), beta0,0,'weights', wcell(1:numel(nP)));
        ci1param = nlparci(beta1param,r,'Jacobian',J);
        bVals = beta1param(2:numel(nP)+1);
        onVals = beta1param(numel(nP)+2:2*numel(nP)+1); 
        JVals = beta1param(2*numel(nP)+2:end);
        bValsErr = ci1param(2:numel(nP)+1,2) - bVals;
        onValsErr = ci1param(numel(nP)+2:2*numel(nP)+1,2) - onVals;
        JValsErr = ci1param(2*numel(nP)+2:end,2) - JVals;

        [res,err] = fit(nP(fitRange).^(-1/dim),bVals(fitRange),'poly1');
        bInf = res.p2;
        [res,err] = fit(nP(fitRange).^(-1/dim),onVals(fitRange),'poly1');
        onInf = res.p2;
        [res,err] = fit(nP(fitRange).^(-1/dim),JVals(fitRange),'poly1');
        JInf = res.p2;
        
        mdl_cell1paramInf = str2func(['@(beta,x) beta(1) * ',num2str(bInf/JInf),' * log(1+exp((x-',num2str(onInf),')/',num2str(bInf),'))']);
        Sigma = Sigma(1,1);
        [ypred1,delta1] = nlpredci(mdl_cell1paramInf,y,beta1param(1),r,'jacobian',J(:,1));
        
        figure(dim*100000)
        clf; hold all
        al = 0.5;
        Legend = cell(numel(nP),1);
        for i = 1:numel(nP) 
            plot(xcell{i},widthcell{i}*nP(i).^al,'.-','Color',[(i-1)/(numel(nP)-1) 0 1-(i-1)/(numel(nP)-1)],'MarkerSize',10,'Linewidth',1.5)
            Legend{i} = ['N=',num2str(nP(i))];
        end
        xlabel('\phi_{eq}')
        ylabel('$N^{1/2}\sigma_{\phi_J}$','interpreter','latex')
        plotSpecs()
        set(gca,'Linewidth',1.5)
        set(gca,'FontSize',20)
        if numel(nP) >= 5
            legend(Legend,'location','northeast','numcolumns',2,'fontsize',16)
        else
            legend(Legend,'location','northeast','numcolumns',1,'fontsize',16)
        end
        
        offsetErr = 1./[1 2 3 4 5 6].^2/2;
        
        figure(1113)
        a = beta1param(1);
        A = ci1param(1,2)-a;
        [res,err] = fit(nP(fitRange).^(-1/dim),bVals(fitRange),'poly1');
        b = res.p2;
        bslope = res.p1;
        conf = confint(res);
        B = conf(2,2) - b;
        [res,err] = fit(nP(fitRange).^(-1/dim),onVals(fitRange),'poly1');
        c = res.p2;
        conf = confint(res);
        C = conf(2,2) - c;
        [res,err] = fit(nP(fitRange).^(-1/dim),JVals(fitRange),'poly1');
        d = res.p2;
        conf = confint(res);
        D = conf(2,2) - d;
        
        phiJInf(dim) = d;
        phiJInfErr(dim) = D;
        phiOnInf(dim) = c;
        phiOnInfErr(dim) = C;
        delta1 = sqrt(a^2*C^2*exp(2*(y-c)/b)./(1+exp((y-c)/b)) + A^2*b^2*log(1+exp((y-c)/b)) + B^2*(a*log(1+ exp((y-c)/b)) - a*exp((y-c)/b).*(y-c)/(b*(1+exp((y-c)/b)))).^2);
%         shadedErrorBar((y)/onInf,ypred1,delta1,'lineprops',{'-','color',dimColor(dim,:),'Linewidth',2})
%         h2 = plot([phid(dim)/onInf phid(dim)/onInf],[0 1],'Linestyle','--','Color',dimColor(dim,:),'Linewidth',2);
%         set(get(get(h2,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');

        figure(dim*10)
        aVal = beta1param(1);
        for j = 1:numel(nP)
            xset = linspace(min(xcell{j}),max(xcell{j}),100); 
            h = plot(xset,aVal*bVals(j)*log(1+exp((xset - onVals(j))/bVals(j))) + JVals(j),'-','Color',[(j-1)/(numel(nP)-1) 0 1 - (j-1)/(numel(nP)-1)],'Linewidth',1.5); 
            set(get(get(h,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
        end
        shadedErrorBar(y,ypred1 + JInf,delta1,'lineprops',{'-k','Linewidth',2})

        x50 = linspace(min(nP.^(-1/dim)),max(nP.^(-1/dim)),100);
        
        figure(333333)
        offsetVal = 1./[1 2 3 4 5 6].^2;
        onsetConversion = [NaN, NaN, 0.9, 0.89, 0.9, 0.87]; % Actually calculating the crossover, convert to onset
        [res4,err] = fit(nP(fitRange).^(-1/dim),onVals(fitRange),'poly1');
        [res,err] = fit(nP(fitRange).^(-1/dim),(res4.p2 - onVals(fitRange))*offsetVal(dim)*onsetConversion(dim),'poly1');
        errorbar(nP.^(-1/dim),(res4.p2 - onVals)*offsetVal(dim)*onsetConversion(dim),onValsErr*onsetConversion(dim)*offsetErr(dim),'.','Color',dimColor(dim,:),'Linewidth',1.5)
        h50 = plot(x50,x50*res.p1,'-','Color',dimColor(dim,:),'Linewidth',1.5);
        set(get(get(h50,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
        
        
        figure(333334)
        offsetVal = 1./[1 2 3 4 5 6].^2;
        [res,err] = fit(nP(fitRange).^(-1/dim),(JInf - JVals(fitRange))*offsetVal(dim),'poly1');
        errorbar(nP.^(-1/dim),(JInf - JVals)*offsetVal(dim),JValsErr*offsetVal(dim),'.','Color',dimColor(dim,:),'Linewidth',1.5)
        h50 = plot(x50,x50*res.p1,'-','Color',dimColor(dim,:),'Linewidth',1.5);
        set(get(get(h50,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
        
%         boffset = [1 1 1 0.7 0.5 0.5];
        boffset = 1./[1 2 3 4 5 6].^2;
        figure(52)
        errorbar(nP.^(-1/dim),(bVals - b)*boffset(dim),bValsErr*boffset(dim),'.','Color',dimColor(dim,:),'Linewidth',1.5)
        h52 = plot(x50,x50*bslope*boffset(dim),'-','Color',dimColor(dim,:),'Linewidth',1.5);
        set(get(get(h52,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
    catch
    end
       
    
    
    
    try
        nP = nP(fitRange);
        y = linspace(0,phid(dim),1000);
        [beta1param,r,J,Sigma,~,~,~] = nlinmultifit(xcell(1:numel(nP)), ycell(1:numel(nP)), mdl_cell1param(1:numel(nP)), beta0,0,'weights', wcell(1:numel(nP)));
        ci1param = nlparci(beta1param,r,'Jacobian',J);
        bVals = beta1param(2:numel(nP)+1);
        onVals = beta1param(numel(nP)+2:2*numel(nP)+1); 
        JVals = beta1param(2*numel(nP)+2:end);
        bValsErr = ci1param(2:numel(nP)+1,2) - bVals;
        onValsErr = ci1param(numel(nP)+2:2*numel(nP)+1,2) - onVals;
        JValsErr = ci1param(2*numel(nP)+2:end,2) - JVals;

        [res,err] = fit(nP.^(-1/dim),bVals(fitRange),'poly1');
        bInf = res.p2;
        [res,err] = fit(nP(fitRange).^(-1/dim),onVals(fitRange),'poly1');
        onInf = res.p2;
        [res,err] = fit(nP(fitRange).^(-1/dim),JVals(fitRange),'poly1');
        JInf = res.p2;
        
        mdl_cell1paramInf = str2func(['@(beta,x) beta(1) * ',num2str(bInf/JInf),' * log(1+exp((x-',num2str(onInf),')/',num2str(bInf),'))']);
        Sigma = Sigma(1,1);
        [ypred1,delta1] = nlpredci(mdl_cell1paramInf,y,beta1param(1),r,'jacobian',J(:,1));
        
        
        figure(1113)
        a = beta1param(1);
        A = ci1param(1,2)-a;
        [res,err] = fit(nP.^(-1/dim),bVals(fitRange),'poly1');
        b = res.p2;
        bslope = res.p1;
        conf = confint(res);
        B = conf(2,2) - b;
        [res,err] = fit(nP(fitRange).^(-1/dim),onVals(fitRange),'poly1');
        c = res.p2;
        conf = confint(res);
        C = conf(2,2) - c;
        [res,err] = fit(nP(fitRange).^(-1/dim),JVals(fitRange),'poly1');
        d = res.p2;
        conf = confint(res);
        D = conf(2,2) - d;
        
        phiJInf(dim) = d;
        phiJInfErr(dim) = D;
        phiOnInf(dim) = c;
        phiOnInfErr(dim) = C;
        delta1 = sqrt(a^2*C^2*exp(2*(y-c)/b)./(1+exp((y-c)/b)) + A^2*b^2*log(1+exp((y-c)/b)) + B^2*(a*log(1+ exp((y-c)/b)) - a*exp((y-c)/b).*(y-c)/(b*(1+exp((y-c)/b)))).^2);
        shadedErrorBar((y)/onInf,ypred1,delta1,'lineprops',{'-','color',dimColor(dim,:),'Linewidth',2})
        h2 = plot([phid(dim)/onInf phid(dim)/onInf],[0 1],'Linestyle','--','Color',dimColor(dim,:),'Linewidth',2);
        set(get(get(h2,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
    catch
    end
    
    
    
    
    figure(dim * 10)
    plotSpecs()
    set(gca,'Linewidth',1.5)
    xlabel('\phi_{eq}')
    ylabel('\phi_{J}')
    set(gca,'FontSize',20)
    if dim == 3
        axis([0.4 0.6 0.63 0.665])
        yticks([0.63 0.64 0.65 0.66])
        xticks([0.4 0.5 0.6])
    elseif dim == 4
        axis([0.2 0.39 0.438 0.46])
    elseif dim == 5
        axis([0.1 0.25 0.296 0.312])
    elseif dim == 6
        axis([0.05 0.15 0.19 0.2])
    else
        axis([0.01 0.05 0.071 0.076])
    end
end

figure(333333)
xlabel('N^{-1/d}')
ylabel('\phi_{on} - \phi_{on}(N)')
% legend('d=3','d=4','d=5','d=6','d=8','location','southwest','numColumns',2)
set(gca,'xscale','log');
set(gca,'yscale','log');
plotSpecs()
pbaspect([2 0.618 1])
xlabel('N^{-1/d}')
ylabel('[\phi_{on} - \phi_{on}(N)]/d^2')
set(gca,'xscale','log');
set(gca,'yscale','log');
plotSpecs()
pbaspect([2 0.618 1])
axis([0.045 0.41 0.0001 0.0044])
set(gca,'FontSize',20)
set(gca,'Linewidth',1.5)
yticks([1e-4 1e-3])
xticks([0.1 0.2 0.3])

figure(333334)
xlabel('N^{-1/d}')
ylabel('[\phi_{J0} - \phi_{J0}(N)]/d^2')
legend('d=3','d=4','d=5','d=6','d=8','location','northwest','numColumns',2)
set(gca,'xscale','log');
set(gca,'yscale','log');
plotSpecs()
pbaspect([2 0.618 1])
axis([0.045 0.41 0.0001 0.003])
set(gca,'FontSize',20)
set(gca,'Linewidth',1.5)
yticks([1e-4 1e-3])
xticks([0.1 0.2 0.3])

figure(1113)
axis([0.6 1.15 -0.001 0.055])
xlabel('\phi_{eq}/\phi_{co}')
ylabel('(\phi_{J} - \phi_{J0})/\phi_{J0}')
legend('d=3','d=4','d=5','d=6','location','southwest','numColumns',2)
plotSpecs()
set(gca,'FontSize',20)
set(gca,'Linewidth',2)



figure(30)
plot([phid(3) phid(3)],[0 1],'--k','Linewidth',1.5)
legend('N=64','N=128','N=256','N=512','N=1024','N=2048','N=4096','N=8192','N=\infty','\phi_d','Numcolumns',2,'location','northwest','fontsize',12)

figure(40)
plot([phid(4) phid(4)],[0 1],'--k','Linewidth',1.5)
axis([0.2 0.42 0.438 0.47])
legend('N=256','N=512','N=1024','N=2048','N=4096','N=8192','N=\infty','\phi_d','Numcolumns',2,'location','northwest','fontsize',12)

figure(50)
plot([phid(5) phid(5)],[0 1],'--k','Linewidth',1.5)
axis([0.1 0.28 0.296 0.316])
legend('N=256','N=512','N=1024','N=2048','N=4096','N=\infty','\phi_d','Numcolumns',2,'location','northwest','fontsize',12)

figure(60)
plot([phid(6) phid(6)],[0 1],'--k','Linewidth',1.5)
axis([0.05 0.18 0.19 0.21])
legend('N=512','N=1024','N=2048','N=4096','N=\infty','\phi_d','Numcolumns',2,'location','northwest','fontsize',12)

figure(52)
xlabel('N^{-1/d}')
ylabel('[b(N) - b(\infty)]/d^2')
legend('d=3','d=4','d=5','d=6','location','southeast','numColumns',2)
set(gca,'xscale','log');
set(gca,'yscale','log');
plotSpecs()
set(gca,'FontSize',20)
set(gca,'Linewidth',1.5)
axis([0.045 0.41 0.00002 0.003])
xticks([0.1 0.2 0.3])

end