%--------------------------------------------------------------------------
%  Gegenbauer Reconstruction of f(x) which is not periodic in x = [-1,1]
%  (Total number of Fourier coefficients = 2N+1 (2Nmax+1)).
%  The convergence constraint [lambda = m = N/4] is used. 
%
%  - Jae-Hun Jung
%
% D. Gottlieb, C.-W. Shu, A. Solomonoff, H. Vandeven,
%        On the Gibbs phenomenon I: Recovering exponential accuracy from the Fourier partial
%        sum of a nonperiodic analytic function,
%        J. Comput. Appl. Math. 43 (1992) 81-92.
% D. Gottlieb, C.-W Shu,
%       On the Gibbs phenomenon and its resolution, SIAM Rev. 39 (1997) 644-668.
%--------------------------------------------------------------------------
 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 = Nmax/4;
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;
% Normalization Factor
   for iy = 0:NL-1
       h(iy+1) = L(iy+1,end)/(lambda+iy);
   end
   h = h*sqrt(pi)*gamma(lambda+0.5)/gamma(lambda);
   
% Transformation Matrix from Gegenbauer space to Fourier space
for iy = 1:NL
  for ix = -Nmax:Nmax
   Tq(iy,ix+Nmax+1) = sum((1-(xi).^2).^(lambda-0.5).*exp(i*ix*pi*xi).*L(iy,:).*omega')/h(iy);
  end
end
%------- Gegenbauer Reconstruction (Projection)
  Ntrunc = Nmax/4;
  C = Tq(1:Ntrunc,:)*a;
  fm = (L(1:Ntrunc,:))'*C;

  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)')
  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')