%Script to create set up an visualize matrix
%to solve for coefficients in rational approximation 

n = 3;
m = 5;
NN = n + 2*m + 1;
N = m + n;
f=@(x) exp(-x).*sin(pi*x);

a = zeros(NN,1);
for k = 1:NN
   a(k) = quad( @(theta) f(cos(theta)).*cos((k-1)*theta), 0, pi, 1e-12) / pi;
end
a(2:end) = 2*a(2:end);

a = a(:);
A(2:end,1) = a(2:N+1);
a(2:end)   = a(2:end)/2;
A(1,1:m+1) = a(1:m+1)';
B = toeplitz(a(1:N),a(1:m));
C = hankel(a(3:N+2),a(N+2:N+m+1));
A(2:end,2:m+1)   = B+C;
A(1:n+1,m+2:N+2) = -eye(n+1);
x = -A(:,2:end)\A(:,1);
q = [1;x(1:m)];
p = x(m+1:m+n+1);

figure(1)
spy(A)
title('Sparsity Pattern of A')
colormap jet
figure(2)
subplot(2,1,1)
imagesc(B);
title('Toeplitz Matrix')
colorbar
subplot(2,1,2)
%use log of entries to see them better
imagesc(log(abs(C)))
title('Hankel Matrix')
colormap hsv
colorbar