load CD.mat;
test = input('1 for match inf or 2 for match 0: ');
bc   = input('1 for original B and C; 2 for random B and C: ');
kmax = 120;
 
if (test == 1),
   [V,W,T,f,g] = NSLanczos(A,kmax,B,C');
end

if (test == 2),
   [V,W,T,f,g] = NSLanczos(inv(A),kmax,B,C');
end
 
if (bc == 2),
    B = randn(120,1);
    C = randn(1,120);
 end

 figure(1), clf
 e1 = 0; e2 = 7;
 Mysigma_log(A,B,C,D,e1,e2,'','b');
 n = size(A,1);
 axis([10^e1 10^e2 1e-4 1e6])
 hold on

 figure(2), clf
 ew = eig(A);
 plot(real(ew),imag(ew),'b.','markersize',10)
% axis equal 
 hold on

 for k=1:kmax
   Ak = W(:,1:k)'*A*V(:,1:k);
   Bk = W(:,1:k)' * B;
   Ck = C*V(:,1:k);
   figure(1)
   textstring = [ 'Frequency Response, CD Player:  n = ', num2str(n),  ' k = ', num2str(k)];
   newfig = Mysigma_log(Ak,Bk,Ck,D,e1,e2,textstring,'r-');
   if k>1, delete(fig); end; fig = newfig;
   axis([10^e1 10^e2 1e-4 1e6])
   legend('Original', 'Reduced')
   figure(2)
   if k>1, delete(ewfig); end
   ew_r = eig(Ak); 
   ewfig = plot(real(ew_r),imag(ew_r),'ro','markersize',15,'linewidth',2);
   title('Eigenvalues of Full and Reduced Models','Fontsize',20)
   drawnow
 end