%% Generate an undecorated version of Figure 3a

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

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

%% Load the cycle-expansion data
% Nn=2; k=3
data_name=sprintf('cycle_lambdas%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); % betas
y2_I=numbers(2*(1:Ndata)-0); % lambdas

% Nn=2; k=7
data_name=sprintf('cycle_lambdas%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); % betas
y2_II=numbers(2*(1:Ndata)-0); % lambdas

% Nn=2; k=11
data_name=sprintf('cycle_lambdas%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); % betas
y2_III=numbers(2*(1:Ndata)-0); % lambdas

% 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_lambdas%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); % betas
y3_I=numbers(2*(1:Ndata)-0); % lambdas

% Nn=3; k=7
data_name=sprintf('cycle_lambdas%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); % betas
y3_II=numbers(2*(1:Ndata)-0); % lambdas

% Nn=3; k=11
data_name=sprintf('cycle_lambdas%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); % betas
y3_III=numbers(2*(1:Ndata)-0); % lambdas


y2_MC=-y2_MC./x2_MC;
y3_MC=-y3_MC./x3_MC;
x2_MC=1./x2_MC;
x3_MC=1./x3_MC;
y2_I=-y2_I./x2_I;
y3_I=-y3_I./x3_I;
x2_I=1./x2_I;
x3_I=1./x3_I;
y2_II=-y2_II./x2_II;
y3_II=-y3_II./x3_II;
x2_II=1./x2_II;
x3_II=1./x3_II;
y2_III=-y2_III./x2_III;
y3_III=-y3_III./x3_III;
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([-1.2,-0.9])
xlabel('$\tilde{T}$','Interpreter','latex')
ylabel('$\tilde{f}$','Interpreter','latex')

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