clear; 

% load matrix data, and exact values of transfer function
load brake834M; 
load brake834K; 
% set the (1,2) and (2,1) blocks of the matrix K to zero, so that
% K has the following structure: 
%
%   K = [ K11  0  K13 ]
%       [ 0   K22 K23 ]
%       [ K31 K32 K33 ]

K12 = K(1:618,619:726) 
K21 = K(619:726,1:618) 

K11 = K(1:618,1:618);
K22 = K(619:726,619:726);
K33 = K(727:834,727:834);

K13 = K(1:618,727:834);
K23 = K(619:726,727:834); 

K31 = K13';
K32 = K23';  

Kr = [ K11 K12 K13
       K21 K22 K23 
       K31 K32 K33 ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C = M;
G = K;
N = size(C,1);
l = zeros(N,1);
l(1)=1;
b = l;


%%%%%%%

% expansion point 
s0 = 0;

% order of the reduced model 
n = 30;

% Reduced-order modeling 

[Cn,Gn,bn,ln] = romlin(N,C,G,b,l,n,s0);

% Bode plot 

npts = 501;
freqlow = 0;
freqhig = 10000;
freq = linspace(freqlow,freqhig,npts);
s = 2*pi*sqrt(-1)*freq;
freqpts = max(size(s));

% evaluate the transfer function of the reduced-order model 

for k = 1:freqpts,
    Hn(k) = ln'*((Gn+s(k)^2*Cn)\bn);
end

% ``exact'' transfer function 
%for k = 1:freqpts,
%    HN(k) = l'*((G+s(k)^2*C)\b);
%end
% precalculated
load brake834freqexact; 

figure(1)
subplot(2,1,1) 
semilogy(abs(freq),abs(HN),'r-');
xlabel('Frequency (Hz)')
ylabel('Magnitude')
hold on;
semilogy(abs(freq),abs(Hn),'b--');
hold off;

subplot(2,1,2) 
plot(abs(freq),angle(HN)*180/pi,'r-');
xlabel('Frequency (Hz)')
ylabel('phase(degree)')
hold on;
plot(freq,angle(Hn)*180/pi,'b--');
hold off;