%--------------------------------------------------------------------------
%  Inverse Polynomial Reconstruction of f(x) which is not periodic in x = [-1,1]
%  (Total number of Fourier coefficients = 2N+1 (2Nmax+1))
%
%  - Jae-Hun Jung
%
% B. D. Shizgal and J.-H. Jung, "Towards the resolution of the Gibbs phenomenon", 
%    Journal of Computational and Applied Mathematics, Vol. 161(1), pp. 41-65.  
% J.-H. Jung and B. D. Shizgal, "Generalization of the Inverse Polynomial Reconstruction Method 
%    in the Resolution of the Gibbs Phenomena", Journal of Computational and Applied Mathematics, 
%    Vol. 172(1), pp. 131-151.  
% J.-H. Jung and B .D. Shizgal, "On the numerical convergence with the inverse polynomial 
%    reconstruction method for the resolution of the Gibbs phenomenon", 
%    Journal of Computational Physics, Vol. 224, pp. 477-488 
%--------------------------------------------------------------------------
 clear all, i = sqrt(-1);

  Nmax = 12; N = Nmax; NL = 2*N+1;

  % Chebyshev Quadrature
  Nq=500;
  xi=-cos(pi*[0:Nq]/Nq);
  omega=xi'*0+pi/Nq;omega(1)=omega(1)/2;omega(Nq+1)=omega(1);
  omega=omega.*sqrt(1-xi'.^2); 
  
% Function to be reconstructed
  fq=xi';  

% Fourier coefficients
  a=zeros(2*Nmax+1,1); 
  for ix=-Nmax:Nmax
     a(ix+Nmax+1) = sum(fq.*omega.*exp(-i*ix*pi*xi'))/2;
  end
  
  % Fourier approximation
    fN = exp(i*(xi')*[-N:N]*pi)*a;

% Gegenbauer polynomials
%lambda = input('lambda ');
lambda = 0.5; % Legendre polynomials

L(1,:)=1.0+ xi*0.0;              %L0
L(2,:)=2.*lambda*xi;             %L1

for k=1:2*N+1-2;
    index=k+1;
    L(index+1,:)= (2*(index+lambda-1)/index)*xi.*L(index,:) - L(index-1,:)*(index+2*lambda-2)/index;
end;

% Transformation Matrix 
for ix = -Nmax:Nmax
     for iy = 1:NL
        Tq(ix+Nmax+1,iy)=sum(L(iy,:)'.*exp(-i*ix*pi*xi').*omega)/2;
     end
end

% Inverse Reconstruction 
  C = Tq\a;
  fm = (L(1:2*N+1,:))'*C;

  
  figure
  subplot(1,2,1), imagesc(real(fN)')
  subplot(1,2,2), imagesc(real(fm)')
  pause
  
  figure
  plot(xi,real(fm),'b',xi, real(fN),'r',xi,fq,'g')
  axis([-1.1, 1.1 -1.1 1.1]), xlabel('x'), ylabel('f(x), f_m(x), f_N(x)')
  
  pause
  figure
  plot(xi,log10(abs(fq-real(fm))),'b', xi, log10(abs(fN-fq)),'r')
  axis([-1.1 1.1 -18 2]), xlabel('x'), ylabel('Pointwise Errors')