function []=figureS3()
close all
figure()
hold on
%% Loading Data Lengths
% PTS length, fitting estimate
data_name=sprintf('cavityPTSthres.dat'); % already converted by conversion_factor=1/1.00663684;
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
Ndata=round(size(numbers,1)/2);
Ts_PTS=numbers(2*(1:Ndata)-1);
xi_PTS=numbers(2*(1:Ndata));

% Positional order length
data_name=sprintf('length_st.dat');
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
Nst=round(size(numbers,1)/3);
numbers=reshape(numbers,3,Nst);
T_st=numbers(1,:);
xi_sts=numbers(2,:);
xi_st_errors=numbers(3,:);

% Bond-orientational order length
data_name=sprintf('length6.dat');
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
N6=round(size(numbers,1)/3);
numbers=reshape(numbers,3,N6);
T6=numbers(1,:);
xi6s=numbers(2,:);
xi6_errors=numbers(3,:);

%% Plotting
plot(1./Ts_PTS,xi_PTS,'Marker','s')
errorbar(1./T_st,xi_sts,xi_st_errors,'Marker','d')
errorbar(1./T6,xi6s,xi6_errors,'Marker','.')
xlim([0,40])
ylim([0,8])
xlabel('$1/T$','Interpreter','latex');
ylabel('$\xi$','Interpreter','latex')
eval(sprintf('print -dpdf figureS3.pdf'))
end