function ps7_3_2(n,k) %Attempt to find a reduced model with stable part of 
                      %order k for a Riemann sum approximation of order n.

C=ones(1,100);       %Create a state-space description of the Ga(s)
B=[];
a=[];
for m=0:99
B=[B; 1/(200+m)];
a=[a -1*(1+m/100)];
end
A=diag(a);
D=0;
Ga=pck(A,B,C,D);


[sysout,sig]=sysbal(Ga,1e-20);  %Balance the system.  The 1e-20 can be set to
                                %any appropriate tolerance param that gives 
                                %enough poles for hankmr to work with!

if (k>length(sig))
     k=length(sig)-1;           %can't choose k greater than what sysbal will
end                             %give us!

sh=hankmr(sysout,sig,k,'d');     %model reduce!
     [Ar,Br,Cr,Dr]=unpck(sh);
     Ghat=tf(ss(Ar,Br,Cr,Dr));
     
w=logspace(-4,4,200);            %numerically compute the infinity norm of 
                                 %G-Ghat.

Ghattemp=freqresp(Ghat,w) ;      %Frequency response of Ghat.
for m=1:200
Ghatfreq(m)=Ghattemp(:,:,m);     %dealing with strange MATLAB notation...
end

o=ones(1,200);                  %Freq response of G minus unstable term.       
Gfreq=o./(j*w-o).*log(2*o.*(j*w+2*o)./(3*o)./(j*w+o));

Gdif=Gfreq+Ghatfreq;            %the plus is correct---remember the minus
                                %sign from the partial fraction expansion!

semilogx(w,abs(Gdif));          %frequency response of difference
grid
title('|G(jw)-Ghat(jw)|')        %plot the frequency response to check the 
xlabel('w')                      %infinity norm of the difference.

s=tf('s')
Ghattf=log(1.5)/(s-1)-Ghat       %the final transfer function