%% Generate an undecorated version of Figure 3b

function []=figure3b()
%% Load the Monte Carlo data
data_name=sprintf('chis%d.dat',2);
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
Ndata=round(size(numbers,1)/2);
x2_MC=numbers(2*(1:Ndata)-1); % betas
y2_MC=numbers(2*(1:Ndata)-0); % chis

data_name=sprintf('chis%d.dat',3);
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
Ndata=round(size(numbers,1)/2);
x3_MC=numbers(2*(1:Ndata)-1); % betas
y3_MC=numbers(2*(1:Ndata)-0); % chis

%% Load the cycle-expansion data
% Nn=2; k=3
data_name=sprintf('cycle_chis%d_%d.dat',2,3);
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
Ndata=round(size(numbers,1)/2);
x2_I=numbers(2*(1:Ndata)-1);
y2_I=numbers(2*(1:Ndata)-0);

% Nn=2; k=7
data_name=sprintf('cycle_chis%d_%d.dat',2,7);
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
Ndata=round(size(numbers,1)/2);
x2_II=numbers(2*(1:Ndata)-1);
y2_II=numbers(2*(1:Ndata)-0);

% Nn=2; k=11
data_name=sprintf('cycle_chis%d_%d.dat',2,11);
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
Ndata=round(size(numbers,1)/2);
x2_III=numbers(2*(1:Ndata)-1);
y2_III=numbers(2*(1:Ndata)-0);

% Restriction to xi<10^5
x2_I=x2_I(1:(Ndata-5)); y2_I=y2_I(1:(Ndata-5));
x2_II=x2_II(1:(Ndata-5)); y2_II=y2_II(1:(Ndata-5));
x2_III=x2_III(1:(Ndata-5)); y2_III=y2_III(1:(Ndata-5));

% Nn=3; k=3
data_name=sprintf('cycle_chis%d_%d.dat',3,3);
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
Ndata=round(size(numbers,1)/2);
x3_I=numbers(2*(1:Ndata)-1);
y3_I=numbers(2*(1:Ndata)-0);

% Nn=3; k=7
data_name=sprintf('cycle_chis%d_%d.dat',3,7);
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
Ndata=round(size(numbers,1)/2);
x3_II=numbers(2*(1:Ndata)-1);
y3_II=numbers(2*(1:Ndata)-0);

% Nn=3; k=11
data_name=sprintf('cycle_chis%d_%d.dat',3,11);
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
Ndata=round(size(numbers,1)/2);
x3_III=numbers(2*(1:Ndata)-1);
y3_III=numbers(2*(1:Ndata)-0);

x2_MC=1./x2_MC;
x3_MC=1./x3_MC;
x2_I=1./x2_I;
x3_I=1./x3_I;
x2_II=1./x2_II;
x3_II=1./x3_II;
x2_III=1./x2_III;
x3_III=1./x3_III;

%% Plot
figure();
hold on
plot(x2_MC,y2_MC);
plot(x2_I,y2_I);
plot(x2_II,y2_II);
plot(x2_III,y2_III);
plot(x3_MC,y3_MC);
plot(x3_I,y3_I);
plot(x3_II,y3_II);
plot(x3_III,y3_III);
xlim([0,1])
ylim([0,10])
xlabel('$\tilde{T}$','Interpreter','latex')
ylabel('$\chi_{\rm L}$','Interpreter','latex')

set(gcf, 'PaperPositionMode', 'auto');
print -depsc2 Figure3b.eps
end