clear; 

% read in data
awebookeg; 
N = size(C,1); 

% expansion point 
s0 = 1*pi*1e+9;
% order of reduced model 
n = 5;    % n=6 gives perfect plots 

% Reduced-order modeling 

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

% Bode plot 

npts = 1001;
freqlow = 10^0;
freqhig = 2*10^(10);
linst = (freqhig-freqlow)/(npts-1);
freq = freqlow:linst:freqhig;
s = 2*pi*sqrt(-1)*freq;
freqpts = max(size(s));

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

% Transfer function - exact 
for k = 1:freqpts,
    HN(k) = l'*((G+s(k)*C)\b);
end

% Transfer function - approximation 
for k = 1:freqpts,
    Hn(k) = ln'*((Gn + s(k)*Cn)\bn);
end

figure(1)
subplot(2,1,1) 
plot(abs(freq),20*log10(abs(HN)),'r-');
xlabel('Frequency (Hz)')
ylabel('Gain (dB)')
hold on;
plot(abs(freq),20*log10(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;