clear all;
close all;

Na = 300; Nuc = Na/2;
b1 = -2.0; b2 = -3.0; a = 1.70;
Fx = 0.02;

H = zeros(Na); V = zeros(Na);
for ai = 1:Na
    if ((ai<Na)&(mod(ai,2)==0)) H(ai,ai+1) = b1; end;
    if ((ai<Na)&(mod(ai,2)==1)) H(ai,ai+1) = b2; end;
    H(ai,ai) = a*Fx/4*(ai-(Na+1)/2);
end
H = H + H';

[V,E] = eig(H); E=real(diag(E)); [E,A] = sort(E); V = V(:,A);
x = a/2*([1:Na]-(Na+1)/2);

hwmin = 1; hwmax = 5; Nw = 1000; hg = 0.05;
hw = linspace(hwmin,hwmax,Nw);
chi = zeros(1,Nw);
Nf = Nuc;

Xmat = zeros(Na); Vx = zeros(Na);
for n = 1:Na
    Vx(:,n) = x(:).*V(:,n);
    for m = 1:Na
        Xmat(m,n) = V(:,m)'*Vx(:,n);
    end
end

Nk = 1000; k1 = 0; k2 = 2*pi/a; k = linspace(k1,k2,Nk);
ecv = 2*sqrt(b1^2+b2^2+2*b1*b2*cos(k*a));
pcv = 1i*a*(b1^2-b2^2)./ecv;

for n = 1:Nf
    for m = Nf+1:Na
        Xnm = Xmat(n,m);
        if (E(m)-E(n)>2*Fx*a) chi(:) = chi(:) + 2*Xnm^2*(E(m)-E(n))./((E(m)-E(n))^2-(hw(:)+1i*hg).^2);
        end;
    end
end
if (Fx<1e-3) chi = chi/(a*Nuc); end;
plot(hw,imag(chi))
hold on;


if (Fx<1e-3)
    for iw = 1:Nw
        chi(iw) = trapz(k(:),abs(pcv(:)).^2./ecv(:)./(ecv(:).^2-(hw(iw)+1i*hg).^2))/pi;
    end
    plot(hw,imag(chi),'g')
end

if (Fx>=1e-3)
    phi = zeros(Nk,1);
    eint = cumtrapz(k,ecv);
    eav = eint(Nk)/(k2-k1);
    i1 = round((hwmin-eav)*(k2-k1)/(2*pi*Fx)); i2 = round((hwmax-eav)*(k2-k1)/(2*pi*Fx));
    for iw = 1:i2-i1+1
        hwp = 2*pi*Fx*(iw+i1-1)/(k2-k1) + eav;
        hw0(iw) = hwp;
        phi = exp(1i*(hwp*(k-k1)-eint)/Fx)/(k2-k1);
        P(iw) = trapz(k,pcv.*phi);
    end
    chi = zeros(1,Nw);
    for iw = 1:i2-i1+1
        chi(:) = chi(:) + abs(P(iw))^2./(hw0(iw)^2-(hw(:)+1i*hg).^2)*2/(a*Fx);
    end
    plot(hw,imag(chi),'g')
end