clc;fontsize=22;

N=40; %used for paper for now

N=100;

figure(4);
fontsize=22;
clf(gcf)


ar_alpha=-[0.1,0.5,1];
letters='abcdef'
for n_alpha=1:numel(ar_alpha)
phi=-pi/2;
alpha=ar_alpha(n_alpha);

gamma=0.3 ;

sigma=0.05;
Nk=500;
export=1;

 path_head='';
 fig_name='FigS9.pdf';
%path_head='';
%%%%%%%%%%%%%%%%%%%%%%%%%%%
H0=zeros(N);
for j=1:N
for k=1:N
	if j==k+1
		H0(j,k)=exp(-1i*alpha);
	end
    	if j==k-1
		H0(j,k)=exp(1i*alpha);
	end
end
end


k=linspace(-pi,pi,Nk);

fk=1i*gamma*(1/pi*(atan((k-phi)/sigma)-atan((k+phi)/sigma)))
%figure(3);
%plot(k,imag(fk))


dk=k(2)-k(1);
U=zeros(N);
for j=1:N
    for j1=1:N
        for n=1:Nk

            U(j,j1)=U(j,j1)+exp(1i*(j-j1)*k(n))*dk/2/pi*fk(n);
        end
    end
end

% figure(11);clf(gcf);
% plot(k,imag(fk))


H=H0 + U;

[V,D]=eig(H);
D=diag(D);
[tmp,ind]=sort(real(D));
D=D(ind);
V=V(:,ind);

[tmp,ind_alpha]=min(abs(k+alpha))
ll=1:numel(k);
ll=circshift(ll',-ind_alpha+round(numel(k)/2));
range_m=ll(1:round(numel(k)/2));
range_p=ll(round(numel(k)/2)+2:numel(k));

k_p=[];
k_m=[];
ipr=[];
for j=1:N
    psi=V(:,j);
    psik=k*0;
    for nk=1:numel(k)
        psik(nk)=sum(psi.*exp(-1i*k(nk)*[1:N]'));
    end
    psik_p=psik(range_p);
    psik_m=psik(range_m);
    [tmp,ind_kp]=max(abs(psik_p));
    [tmp,ind_km]=max(abs(psik_m));

    k_p(j)=k(range_p(ind_kp));
    k_m(j)=k(range_m(ind_km));

    %if j==30
    %    figure(6);clf(gcf);
    %    subplot(1,2,1)
    %    plot(1:N,abs(V(:,j)).^2)
    %    subplot(1,2,2)
    %    plot(k(range_p),abs(psik_p));hold on;
    %    plot(k(range_m),abs(psik_m))
    %end
    %

    ipr(j)=sum(abs(psi.^4))/sum(abs(psi.^2))^2;
end
sp=[];sm=[];
for j=1:N
    sm(j)=sum(abs(V(1:N/2,j).^2));
    sp(j)=sum(abs(V(N/2+1:N,j).^2));
end

%figure(3);
%plot(k,abs(psik));



Ek=exp(-1i*alpha)*exp(-1i*k)+  exp(1i*alpha)*exp(1i*k);
[tmp,ind_p]=min(abs(k-phi));
[tmp,ind_m]=min(abs(k+phi));
Ep=Ek(ind_p)
Em=Ek(ind_m)

Ek1=Ek+fk;


subplot(numel(ar_phi),2,1+2*(n_alpha-1))
plot(k,real(Ek),'-b','linewidth',3);hold on;

plot(k,Ep+Ek*0,'--k','linewidth',1);hold on;
plot(k,Em+Ek*0,'--k','linewidth',1);hold on;

xlim([-pi,pi])

set(gca,'fontsize',fontsize)
set(gca,'linewidth',1)
set(gca,'ticklength',[0.03,0.03]/2)
xlabel('$K$','fontsize',fontsize,'interpreter','latex')
ylabel('Re $\varepsilon_{\rm pair}(K)$','fontsize',fontsize,'interpreter','latex')


plot(k_p,real(D),'vr','markersize',10,'linewidth',1);
plot(k_m,real(D),'ok','markersize',10,'linewidth',1);

set(gca,'xtick',[-pi,0,pi]);
set(gca,'xticklabels',{'-\pi','0','\pi'});

yl=ylim;
fill([-phi,phi,phi,-phi],[yl(2),yl(2),yl(1),yl(1)],'b','facealpha',0.4, 'EdgeColor', 'none')

title(['(',letters(2*n_alpha-1),') $\Phi=',sprintf('%2.1f',alpha ),'$'],'fontsize',fontsize,'interpreter','latex')
set(gca,'titlehorizontalalignment','left')

subplot(numel(ar_alpha),2,2+2*(n_alpha-1))

[X,Y]=meshgrid(1:N,real(D));
pcolor(X,Y,abs(V.^2)');
shading faceted;
caxis([0,0.04])
xlabel('$k$','fontsize',fontsize,'interpreter','latex')
set(gca,'ydir','normal')
cb=colorbar('ticklength',0.03);colormap(hot);
%ylabel(cb,'$|\psi_n^2|$','fontsize',fontsize,'interpreter','latex')
colorTitleHandle = get(cb,'Title');
titleString = '$|\psi_n^2|$';
set(colorTitleHandle ,'String',titleString,'fontsize',fontsize,'interpreter','latex');
hold on;
set(gca,'linewidth',1)
set(cb,'linewidth',1)
set(gca,'ticklength',[0.03,0.03]/2)
set(gca,'tickdir','out')

set(gca,'fontsize',fontsize)
set(gca,'linewidth',1)
set(gca,'ticklength',[0.03,0.03]/2)
xlabel('$n$','fontsize',fontsize,'interpreter','latex')
%ylabel('Re $\varepsilon(K)$','fontsize',fontsize,'interpreter','latex')

xl=xlim;
plot(xl,Ep+xl*0,'--w','linewidth',1);hold on;
plot(xl,Em+xl*0,'--w','linewidth',1);hold on;

title(['(',letters(2*n_alpha),')'],'fontsize',fontsize,'interpreter','latex')
set(gca,'titlehorizontalalignment','left')
end

set(gca,'yticklabels',[])
%annotation('textbox',[0.06 0.9 0.3,0.05],'string','{(a)}','fontsize',fontsize,'Interpreter','Latex','LineStyle','none');
%annotation('textbox',[0.52  0.9 0.3,0.05],'string','{(b)}','fontsize',fontsize,'Interpreter','Latex','LineStyle','none');



set(gcf,'Position',[  1052         105         742         872])

if export
set(gcf, 'Renderer', 'painters'); % Use vector graphics renderer
%print('-depsc','Fig3.eps');system('epstopdf Fig3.eps');
set(gcf,'PaperUnits','centimeters','PaperPositionMode','manual');
exportgraphics(gcf, [path_head,fig_name], 'ContentType', 'vector');
%print('-dpdf','Fig3.pdf')

end
