clear; 

% load matrix data, and exact values of transfer function
load brake834M; 
load brake834K; 
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 = 40;

% 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); 
semilogy(abs(freq),abs(HN-Hn)./abs(HN),'b-');
xlabel('Frequency (Hz)')
ylabel('Relative error')