clear; 

% load data
load ex306v4;   % matrix data, and exact values of transfer function
C = CM;
G = GM;
l = LM(:,1);
b = LM(:,2);
[n,n]=size(C); 


% order of the reduced model 
kmax = 60;

% expansion point 
s0 = 2*pi*1e+9;

% reduced-order modeling 

[Ck,Gk,bk,lk]=romlin(n,C,G,b,l,kmax,s0); 

% Bode plot 

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

% ``exact'' transfer function

tic;
for k = 1:freqpts,
    HN(k) = l'*((G+s(k)*C)\b);
end
Fullmodel_cputime=toc
% precalculated
% load ex306v4HN;

% transfer function of reduced-order model 

tic; 
for k = 1:freqpts,
    Hn(k) = lk'*((Gk+s(k)*Ck)\bk);
end
Reducedordermodel_cputime = toc

% This is not exactly bode polt, but commonly as presented in literature
figure(1)     
plot(abs(freq),abs(HN),'r-');
xlabel('Frequency (Hz)')
ylabel('Voltage (Volts)')
hold on;
plot(abs(freq),abs(Hn),'b--');
hold off;