%% Generate an undecorated version of Figure 2c

function []=Figure2c()
angular_index=8; % angular momentum
ylimhigh=0.06; % y-axis range

figure();
hold on
for packing_index=1:5
    % packing_index=1 corresponds to packing fraction phi1=0.600;
    % packing_index=2 corresponds to packing fraction phi2=0.650;
    % packing_index=3 corresponds to packing fraction phi3=0.680;
    % packing_index=4 corresponds to packing fraction phi4=0.690;
    % packing_index=5 corresponds to packing fraction phi5=0.695;
    
    %% Read in the data
    data_name=sprintf('./l=%dPTS_phi%d.dat',angular_index,packing_index); fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
    Nr=round(size(numbers,1)/3);% number of radial data points
    numbers=reshape(numbers,3,Nr);
    R=numbers(1,1:Nr).'; % cavity radii where bond-orientational point-to-set correlation functions are measured
    gPTS=numbers(2,1:Nr).'; % bond-orientational point-to-set correlation functions
    gPTS_error=numbers(3,1:Nr).'; % error bars
    
    %% Exponential fit, A*exp(-(r/xiPTS))
    x=R; y=gPTS;
    Init=[ylimhigh, 1]; weight=diag(ones(size(x))); Outliers=false(size(x));
    [fit_values,confidence_intervals]=exponential_decay_fit(x,y,weight,Init,Outliers);
    A=fit_values(1);
    xiPTS=fit_values(2);
    xiPTS_error=1.96*(confidence_intervals(2,2)-confidence_intervals(1,2))/2;  % 95-percent confidence interval
    R_fit=0:0.01:20;
    gPTS_fit=A*exp(-(R_fit/xiPTS));
    
    %% Plot them
    plot(R_fit,gPTS_fit)
    errorbar(R,gPTS,gPTS_error,'Linestyle','none','Marker','.')

end
xlabel('$R$','Interpreter','latex')
ylabel('$g^{\rm PTS}_{8}$','Interpreter','latex');
xlim([0,20])
ylim([0,ylimhigh])
axis square
set(gcf, 'PaperPositionMode', 'auto');
print -depsc2 Figure2c.eps
    
end