%
%  w11tch.m -- djm, 11 july 2005
%

N = 64;

%  define tchebyshev points
dy = pi/N;  y = 0:dy:pi;  x = cos(y);

f = cos(8*x.^2 - x);
fp = -(16*x - 1).*sin(8*x.^2 - x);
%f = x.^5 - x.^2 + 1;
%fp = 5*x.^4 - 2*x;
%  t5
%f = (16*x.^5 - 20*x.^3 + 5*x);
%  t6
%f = (32*x.^6 - 48*x.^4 + 18*x.^2 - 1);
f2 = [f fliplr(f(2:N))];  y2 = 0:dy:2*pi-dy;

figure(1);  clf;
subplot(2,1,1);  hold on
plot(x,f,'b',x,f,'k.',x,fp,'r')
axis([-1 1 -16 17])
title('\bf tchebychev discretization for cos(8x^2 - x)')
xlabel('\bf x-axis')
ylabel('\bf f(x)')

subplot(2,1,2);  hold on
plot(y2,f2,'b',y2,f2,'k.')
axis([0 2*pi -1.1 1.1])
title('\bf equi-spaced grid (even extension)')
xlabel('\bf y-axis (y = cos^{-1}x)')
ylabel('\bf f(cos y)')

%  fft for cosine coefficients
cj = real(fft(f2))/N;  tj = cj(1:N);

%  differentiation
kk = [0:N-1 -N:-1];
fp = N*real(ifft((i*cj.*kk)));  fp = [fp(N+1:end) fp(1)];
fp(2:end-1) = fp(2:end-1)./sin(y(2:end-1));

fp(end) = sum((0:N-1).^2 .* tj);
fp(1) = -sum((-1).^(0:N-1) .* (0:N-1).^2 .* tj);

subplot(2,1,1)
plot(x,fliplr(fp),'kx')