%% Generate an undecorated version of the bottom row of Figure 1

function []=Figure1bottom()
cavity_size=6.25;
%% Set figure parameters
L=20; % linear size of the figure
Lh=L/2;
Ngrid=101; % number of grid data point
% Coloring
Ncolor=1000;
color_map=hot(Ncolor);
color_map=1-color_map;
color_map(1:300,1)=1; color_map(1:300,2)=1; color_map(1:300,3)=1;

%% Positional overlap profile
% Read in the data
data_name=sprintf('averagepos.dat');
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
numbers=(reshape(numbers,3,Ngrid^2));
xs=numbers(1,1:(Ngrid^2)).'; xs=reshape(xs,Ngrid,Ngrid); % x-coordinate
ys=numbers(2,1:(Ngrid^2)).'; ys=reshape(ys,Ngrid,Ngrid); % y-coordinate
qs=numbers(3,1:(Ngrid^2)).'; qs=reshape(qs,Ngrid,Ngrid); % Q-values

% Focus on values between 0 and 1
I=logical(qs>1); qs(I)=1;
I=logical(qs<0); qs(I)=0;

% Smoothen the surface with creating the finer grid and then interpolating
x_grid=xs(1,2)-xs(1,1); x_finegrid=x_grid/10; x_refine=(xs(1,1):x_finegrid:xs(1,end));
y_grid=ys(2,1)-ys(1,1); y_finegrid=y_grid/10; y_refine=(ys(1,1):y_finegrid:ys(end,1));
[X,Y]=meshgrid(x_refine,y_refine);
Q=interp2(xs,ys,qs,X,Y);

% Trivially extend the profile to the figure range, [-Lh,Lh]*[-Lh,Lh],
x_extend=((-Lh):x_finegrid:Lh); % x-extended grid
y_extend=((-Lh):y_finegrid:Lh); % y-extended grid
[X,Y]=meshgrid(x_extend,y_extend);
[~,x_begin]=min((x_extend-xs(1,1)).^2); % x-coodinate value where we should begin embeding what we have obtained above
[~,x_end]=min((x_extend-xs(1,end)).^2); % x-coodinate value where we should end embeding what we have obtained above
[~,y_begin]=min((y_extend-ys(1,1)).^2); % y-coodinate value where we should begin embeding what we have obtained above
[~,y_end]=min((y_extend-ys(end,1)).^2); % y-coodinate value where we should end embeding what we have obtained above
Q_ex=ones(size(X));
Q_ex(x_begin:1:x_end,y_begin:1:y_end)=Q; % embedding what we have obtained above
Q=Q_ex;

% Set the color outside the cavity to black
rs=sqrt(X.^2+Y.^2);
I=logical(rs>cavity_size);
Q(I)=1;

% Adjust for coloring purposes
I=logical(Q<=0.30); Q(I)=0;

% Generate figure and print it out
figure(1);
colormap(color_map)
surf(X,Y,Q,Q,'EdgeColor','none','FaceColor','interp')
view([0,90])
xlim([-Lh,Lh])
ylim([-Lh,Lh])
axis square off
set(gcf, 'PaperPositionMode', 'auto');
print -depsc2 Figure1bottomfirst.eps
pause(1)

%% l=6 overlap profile
% Read in the data
data_name=sprintf('average6.dat');
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
numbers=(reshape(numbers,3,Ngrid^2));
xs=numbers(1,1:(Ngrid^2)).'; xs=reshape(xs,Ngrid,Ngrid); % x-coordinate
ys=numbers(2,1:(Ngrid^2)).'; ys=reshape(ys,Ngrid,Ngrid); % y-coordinate
qs=numbers(3,1:(Ngrid^2)).'; qs=reshape(qs,Ngrid,Ngrid); % Q-values

% Focus on values between 0 and 1
I=logical(qs>1); qs(I)=1;
I=logical(qs<0); qs(I)=0;

% Smoothen the surface with creating the finer grid and then interpolating
x_grid=xs(1,2)-xs(1,1); x_finegrid=x_grid/10; x_refine=(xs(1,1):x_finegrid:xs(1,end));
y_grid=ys(2,1)-ys(1,1); y_finegrid=y_grid/10; y_refine=(ys(1,1):y_finegrid:ys(end,1));
[X,Y]=meshgrid(x_refine,y_refine);
Q=interp2(xs,ys,qs,X,Y);

% Trivially extend the profile to the figure range, [-Lh,Lh]*[-Lh,Lh],
x_extend=((-Lh):x_finegrid:Lh); % x-extended grid
y_extend=((-Lh):y_finegrid:Lh); % y-extended grid
[X,Y]=meshgrid(x_extend,y_extend);
[~,x_begin]=min((x_extend-xs(1,1)).^2); % x-coodinate value where we should begin embeding what we have obtained above
[~,x_end]=min((x_extend-xs(1,end)).^2); % x-coodinate value where we should end embeding what we have obtained above
[~,y_begin]=min((y_extend-ys(1,1)).^2); % y-coodinate value where we should begin embeding what we have obtained above
[~,y_end]=min((y_extend-ys(end,1)).^2); % y-coodinate value where we should end embeding what we have obtained above
Q_ex=ones(size(X));
Q_ex(x_begin:1:x_end,y_begin:1:y_end)=Q; % embedding what we have obtained above
Q=Q_ex;

% Set the color outside the cavity to black
rs=sqrt(X.^2+Y.^2);
I=logical(rs>cavity_size);
Q(I)=1;

% Adjust for coloring purposes
I=logical(Q<=0.30); Q(I)=0;

% Generate figure and print it out
figure(6);
colormap(color_map)
surf(X,Y,Q,Q,'EdgeColor','none','FaceColor','interp')
view([0,90])
xlim([-Lh,Lh])
ylim([-Lh,Lh])
axis square off
set(gcf, 'PaperPositionMode', 'auto');
print -depsc2 Figure1bottomsecond.eps
pause(1)

%% l=8 overlap
% Read in the data
data_name=sprintf('average8.dat');
fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
numbers=(reshape(numbers,3,Ngrid^2));
xs=numbers(1,1:(Ngrid^2)).'; xs=reshape(xs,Ngrid,Ngrid); % x-coordinate
ys=numbers(2,1:(Ngrid^2)).'; ys=reshape(ys,Ngrid,Ngrid); % y-coordinate
qs=numbers(3,1:(Ngrid^2)).'; qs=reshape(qs,Ngrid,Ngrid); % Q-values

% Focus on values between 0 and 1
I=logical(qs>1); qs(I)=1;
I=logical(qs<0); qs(I)=0;

% Smoothen the surface with creating the finer grid and then interpolating
x_grid=xs(1,2)-xs(1,1); x_finegrid=x_grid/10; x_refine=(xs(1,1):x_finegrid:xs(1,end));
y_grid=ys(2,1)-ys(1,1); y_finegrid=y_grid/10; y_refine=(ys(1,1):y_finegrid:ys(end,1));
[X,Y]=meshgrid(x_refine,y_refine);
Q=interp2(xs,ys,qs,X,Y);

% Trivially extend the profile to the figure range, [-Lh,Lh]*[-Lh,Lh],
x_extend=((-Lh):x_finegrid:Lh); % x-extended grid
y_extend=((-Lh):y_finegrid:Lh); % y-extended grid
[X,Y]=meshgrid(x_extend,y_extend);
[~,x_begin]=min((x_extend-xs(1,1)).^2); % x-coodinate value where we should begin embeding what we have obtained above
[~,x_end]=min((x_extend-xs(1,end)).^2); % x-coodinate value where we should end embeding what we have obtained above
[~,y_begin]=min((y_extend-ys(1,1)).^2); % y-coodinate value where we should begin embeding what we have obtained above
[~,y_end]=min((y_extend-ys(end,1)).^2); % y-coodinate value where we should end embeding what we have obtained above
Q_ex=ones(size(X));
Q_ex(x_begin:1:x_end,y_begin:1:y_end)=Q; % embedding what we have obtained above
Q=Q_ex;

% Set the color outside the cavity to black
rs=sqrt(X.^2+Y.^2);
I=logical(rs>cavity_size);
Q(I)=1;

% Adjust for coloring purposes
I=logical(Q<=0.30); Q(I)=0;

% Generate figure and print it out
figure(8);
colormap(color_map)
surf(X,Y,Q,Q,'EdgeColor','none','FaceColor','interp')
view([0,90])
xlim([-Lh,Lh])
ylim([-Lh,Lh])
axis square off
set(gcf, 'PaperPositionMode', 'auto');
print -depsc2 Figure1bottomthird.eps
pause(1)


end