%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% BOUNDARY GR + CHERN INSULATOR
% M. Alvarado & A. Levy Yeyati

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all;
ui = sqrt(-1);

global N Mass

Mass = .5;

n_w = 101;          % Frequency points
n_Q = 101;           % Momenta points

Emax = 1.5;
Emin = 1.5;

N = 2;

%%%%%%%%%%%%%%%

% Periodic Boundary Bloch Phase

dQ = 2*pi/n_Q;
momenta_Q = linspace(-pi, pi, n_Q);

%%%%%%%%%%%%%%%

wmax = Emax;
wmin = -Emin;
freqs = linspace(wmin, wmax, n_w);
eta = 2*(wmax-wmin)/n_w;

%%%%%%%%%%%%%%%

Z = zeros(N);
LDOS_0 = zeros(n_Q, n_w);
LDOS_L = zeros(n_Q, n_w);
LDOS_R = zeros(n_Q, n_w);

%%%%%%%%%%%%%%%

%
% set up the calculation
%

tStart = tic;

disp(' ');
disp('  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
disp(' ');
disp('  	LDOS Calculation');
disp(' ');
disp('  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
disp(' ');

idx_Q = 0;
for k = momenta_Q
    idx_Q = idx_Q+1;
    disp(['      Progress Bar (%) = ', num2str(idx_Q*100/n_Q,'%1.2f')]);
    
    [c, rankH3, M] = fun_FFA(k);

    iw = 0;
    for w=freqs
        iw = iw + 1;
        W = w - ui*eta;
        
        P = fun_polynomial(W, c, rankH3);
        z = roots(P);

        Green_0_0 = fun_residues(z, P, M, W, 0, 0);  
        Green_1_0 = fun_residues(z, P, M, W, 1, 0);
        Green_0_1 = fun_residues(z, P, M, W, 0, 1);
      
        g_L = Green_0_0 - Green_0_1/Green_0_0*Green_1_0;
        g_R = Green_0_0 - Green_1_0/Green_0_0*Green_0_1;
        
        LDOS_0(idx_Q, iw) = trace(imag(Green_0_0))/pi;
        LDOS_L(idx_Q, iw) = trace(imag(g_L))/pi;
        LDOS_R(idx_Q, iw) = trace(imag(g_R))/pi;
    end
end


disp(' ');
disp('------------------------------------------------');
disp(' ');

tEnd = toc(tStart);
fprintf('Elapsed time is :  %d min and %2.3f sec\n',floor(tEnd/60),rem(tEnd,60));

disp(' ');
disp('------------------------------------------------');
disp(' ');


%%


figure(7)

ax_a = gca;

Q = linspace(-1,1,n_Q);

pcolor(Q , freqs, LDOS_R');
% pcolor(Q , freqs, LDOS_L');
% pcolor(Q , freqs, LDOS_0');

% colormap(parula)
shading interp


y_cota = get(gca,'ylim');
x_cota = get(gca,'xlim');

hold on;
plot([0,0],y_cota, 'Color','w', 'Linestyle','--','LineWidth', 1.)
plot(0.5*[1,1],y_cota, 'Color','w', 'Linestyle','--','LineWidth', 1.)
plot(-0.5*[1,1],y_cota, 'Color','w', 'Linestyle','--','LineWidth', 1.)
plot(x_cota,[0,0], 'Color','w', 'Linestyle','--','LineWidth', 1.)
hold off

ax_a = gca;
c = colorbar;

LDOS_max = 6;
caxis(ax_a,[0,LDOS_max]);
title(c,'$\rho_R(k_\parallel)$','Interpreter','latex');
set(c, 'TickLabelInterpreter', 'latex');
set(c, 'fontsize', 12);


box on;
xticks([-1,0,1])
xlabel('$k_\parallel/\pi$','FontSize', 20,'Interpreter','latex');

ylabel('$\omega/t$','FontSize', 20,'Interpreter','latex');
yticks([-1,0,1])
set(gca, 'fontsize', 17)
set(gca, 'TickLabelInterpreter', 'latex')

axis([x_cota, y_cota])

