%
%  This is a test driver for MOR via 
%  IRKA.   It is applied
%  to a CD-Player Model of order n = 120
%  It should produce a reduced model of
%  order k = 20.  Bode plots comparing the
%  output of the reduced vs full system are
%  produced.   Also, eigenvalue distribution
%  of the reduced model is plotted.
%

%clear;

%
%  This will input matrices A, B, C for CD-Player
%
   load CD.mat;
%  C = B';
%
%  Set k = requested order of the reduced model
%
%   k = 20;
    k = input('input k value ')
% 
%  A positive real number s1 is input and
%  
%    shift = s1*i  = s1*sqrt(-1)   is constructed 
%                                  and displayed
%

   n = size(A,1); I = eye(n);
   tol = 1.0e-5;

%
%  Run the Arnoldi process for k-steps  to
%  construct reduced basis V
%

   [Ak,Bk,Ck] = IrkaN(A,B,C,k, tol);
%
%  Construct the k-th order reduced model  
%


 sys = ss(A,B,C,0);
 sys_r = ss(Ak,Bk,Ck,0);

%
%  Display Bode plots from Matlab sigma function
%
   figure(1)
   clf
   sigma(sys)
   hold
   sigma(sys_r,'r')
   hold
   legend('full system','reduced system')
   title('CD Player')
   
   n = size(A,1);
   k =  size(Ak,1);
%
%  Display Bode plots from our own routine 
%
   
   figure(2)
   clf
   e1 = 0; e2 = 7;
   textstring = [ 'Freq-Response CD-Player :  n = ', num2str(n),  ' k = ', num2str(k)];
   Mysigma_log(A,B,C,D,e1,e2,textstring,'b');
   hold
   Mysigma_log(Ak,Bk,Ck,D,e1,e2,textstring,'r-.');
   legend('Original', 'Reduced')
   
   Dimension_Full_System = size(A,1) 
   Dimension_Reduced_System =  size(Ak,1)
%
%  Display eigenvalue plots for original and reduced systems
%

   t = eig(A); s = eig(Ak);
   figure(3)
   plot(real(t), imag(t), '*' ,'linewidth',2.5);
   hold on
   
   plot(real(s), imag(s), 'ro','linewidth',2.5);
   
   hold off