clear
[s1, s2, s3, s0] = Utility.PauliMats{:};

mu1 = 1;
% mu2 = 3;
t = 0.9;
m = 1.2;
d = 0.05;
V = 0.3;

ham = @(k,mu2) kron([1,0;0,0],kron(s3,k*s1-mu1*s0)+kron(s1,d*s0))...
          +kron([0,0;0,1],-(t*k^2+mu2)*kron(s3,s0)+m*kron(s0,s3))...
          +V*kron(s1,kron(s3,s0));
      
kmin = -5;
kmax = 5;
kpts = 1001;
ks = linspace(kmin,kmax,kpts);

%%
% mu2 = 3;
% ens = zeros(8,kpts);
% for ik = 1:kpts
%     ens(:,ik) = Utility.eigen(ham(ks(ik),mu2));
% end
% figure
% plot(ks,ens)

%%
mu2min = -2;
mu2max = 3;
mu2pts = 101;
mu2s = linspace(mu2min,mu2max,mu2pts);

ps = zeros(8,8,mu2pts);
gs = zeros(mu2pts,1);
for imu = 1:mu2pts
    tic
    p = 1;
    g = 100;
for ik = 1:kpts
    [evecs, evals] = Utility.eigen(ham(ks(ik),mu2s(imu)));
    p = p*evecs(:,1:4)*evecs(:,1:4)';
    g = min(g,evals(5)-evals(4));
end
    ps(:,:,imu) = p;
    gs(imu) = g;
    toc
end

%%
p0 = ps(:,:,end);
wf0 = evecs(:,1:4);
ees = zeros(4,mu2pts-1);
for imu = 1:mu2pts-1
    ees(:,imu) = eig((p0*wf0)'*ps(:,:,imu)*wf0);
end

%%
figure
plot(mu2s,gs*(1/(2*d)),'-kx')
ns = sum(angle(ees),1)/pi;
hold on
sn = sign((mu1^2+d^2)*(mu2s.^2-m^2)-2*mu1*mu2s*V^2+V^4);
plot(mu2s,angle(sn)/pi,'ro')
% plot(mu2s(1:end-1),ns,'ro')
hold off
xlabel('\mu_d')
ylabel('Gap/2\Delta / \nu')
legend('Gap','Topo. Number')