%% Generate an undecorated version of Figure 3c

function []=Figure3c()
ylimhigh=0.06; % y-axis range
angular_index=8; % angular_momentum

figure();
hold on
for temperature_index=1:5
    % temperature_index=1 corresponds to temperature T1=1.000;
    % temperature_index=2 corresponds to temperature T2=0.800;
    % temperature_index=3 corresponds to temperature T3=0.600;
    % temperature_index=4 corresponds to temperature T4=0.510;
    % temperature_index=5 corresponds to temperature T5=0.450;
    
    %% Read in the data
    data_name=sprintf('./l=%dPTS_T%d.dat',angular_index,temperature_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).'; % positional bond-orientational correlation functions
    gPTS_error=numbers(3,1:Nr).'; % error bars
    
    %% Stretched exponential fit, A*exp(-(r/xiPTS)^3)
    gamma=3;
    x=R; y=gPTS;
    Init=[ylimhigh, 1]; weight=diag(ones(size(x))); Outliers=false(size(x));
    [fit_values,confidence_intervals]=fixed_streched_exponential_decay_fit(x,y,weight,Init,Outliers,gamma);
    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).^gamma);    

    %% 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([1,4])
ylim([0,ylimhigh])
axis square
set(gcf, 'PaperPositionMode', 'auto');
print -depsc2 Figure3c.eps
    
end