function []=figureS15b()
close all
figure()
hold on
%% Data parameters
conversion_factor=1;
sample_indices=[0,1,2,3,4];
phis=[0.700, 0.740,0.800,0.820,0.840]; % area fraction
Zs=[11.01665,14.5729,24.50337,30.35702,39.50089]; % reduced pressure
Ns=size(sample_indices,2);
xipeakPTS=zeros(size(sample_indices));

%% 
for sample_index_dummy=1:Ns
    sample_index=sample_indices(sample_index_dummy);
    %% Color
    colorB=(sample_index_dummy-1)/(Ns-1);
    colorR=1-colorB;
    color=[colorR,0,colorB];
    %% Loading data
    % PTS Susceptibility    
    data_name=sprintf('./chiTs%d.dat',sample_index);
    fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
    NCP=round(size(numbers,1)/3);
    cavity_sizes=numbers(3*(1:NCP)-2);
    chiTs=numbers(3*(1:NCP)-1);
    chiTs_errors=numbers(3*(1:NCP));
    bx=conversion_factor*cavity_sizes;
    by=chiTs;
    be=chiTs_errors;
    
    %% Fitting (not plotted)
    bx_fit=bx;
    by_fit=by;
    Npick=5;
    nearmax_x=zeros(1,Npick);
    nearmax_y=zeros(1,Npick);
    for i=1:Npick
        [~,where]=max(by_fit);
        nearmax_x(i)=bx_fit(where);
        nearmax_y(i)=by_fit(where);
        deletion=ones(size(bx_fit));
        deletion(where)=0;
        deletion=logical(deletion);
        bx_fit=bx_fit(deletion);
        by_fit=by_fit(deletion);
    end
    p = polyfit(nearmax_x,nearmax_y,2);
    a=p(1);
    b=p(2);
    %     c=p(3);
    x_peak=-b/(2*a);
    %     y_peak=c-(b^2/(4*a));
    xipeakPTS(sample_index_dummy)=x_peak;
    
    %% Plotting
    errorbar(bx,by,be,'Color',color,'Marker','.','MarkerSize',20)
    xlim([0,4])
    ylim([0,0.065])
    xlabel('$R$','Interpreter','latex');
    ylabel('$\chi_{\mathrm{PTS}}$','Interpreter','latex')
end
print(gcf, '-dpdf', 'figureS15b.pdf');
end