function ps61_6242_2004(n)
% function ps61_6242_2004(n)
%
% solves problem 6.1(b) from 6.242/Fall2004

if nargin<1, n=4; end  % default argument
tol=1e-13;              % tolerance
N=10000;               % number of samples for H-Infinity norm calculation

a=[(1./(1:floor(n/2)))';(2:(ceil(n/2)+1))'];         % the vector of 1/k (k=1,...,n)
%a=1./(1:n)';
g=(a.^(1/3))./(1+a);   % vector of G(1/k)
W=1./(repmat(a,1,n)+repmat(a',n,1));      % W
[U,D]=eig(W);          % W=UDU' where U is orthogonal, D diagonal
d=diag(D);  
d=max(d,tol);          % regularization
ug=U'*g;               % solving for c
c=U*(ug./d); 
H2Ghat=sum((abs(ug).^2)./d); % square of H2 norm of Ghat
%c=W\g; H2Ghat=g'*c;



w=linspace(0,10,N)';  % column vector of frequency samples
Ghat=(1./(repmat(a',N,1)+repmat(j*w,1,n)))*c;  % samples of Ghat
G=((j*w).^(1/3))./(1+j*w);                    % samples of G

disp(['numeric H2 error: ',num2str((10/N)*sum(abs(G-Ghat).^2))])
disp(['H2 error:         ',num2str(1-H2Ghat)])
disp(['H-Infinity error: ',num2str(max(abs(G-Ghat)))])
close(gcf)
subplot(2,1,1); plot(w,real(G),w,real(Ghat)); grid
subplot(2,1,2); plot(w,imag(G),w,imag(Ghat)); grid