%
%  w09leg.m -- djm, 05 july 2005
%

clear
N = 64;  dx = 2*pi/N;

%  matlab fft works well on 0 < x < 2*pi

x = (0:dx:2*pi-dx);  z = 0.6;
k0 = (-N/2:N/2-1);
k1 = fftshift(k0);

%  single fourier mode (note normalizations of fft & ifft !!!)
j = 1;  a = 10^(-log10(1e5*eps)/N);
y0 = exp(i*j*x);
y0 = 1./sqrt(1 - 2*z*y0/a + y0.^2/(a^2));

yf = fft(y0)/N;
y1 = ifft(yf)*N;

figure(1);  clf
subplot(2,1,1);  hold on
plot(x,real(y1),'y',x,imag(y1),'y')
plot(x,real(y0),'r--',x,real(y0),'k.',x,imag(y0),'b--',x,imag(y0),'k.')
pm = max(abs(y0))*1.05;
axis([0 2*pi -pm pm])
title('\bf real (r) & imag (b) parts of y(x)')
subplot(2,1,2);  hold on
P = (a.^(0:N-1)).*yf;
plot(0:N-1,real(P),'r',0:N-1,real(P),'k.')
axis([0 N -0.5 1.05])
title(['\bf legendre polynomials, P_j(' num2str(z) ')'])

M = N-1;

Q = zeros(N,1);
for j = 0:M
    a = legendre(j,z);
    Q(j+1) = a(1);
end

plot(0:N-1,Q,'mo')

figure(2);  clf
plot(0:N-1,log(abs(real(P)-Q')),'.')

%  recursion
R = zeros(size(Q));  S = zeros(size(Q));
R(1) = 1;  R(2) = z;  S(1) = 1;  S(2) = cos(z);
for j=1:M-1
    R(j+2) = ( (2*j+1)*z*R(j+1) - j*R(j) ) / (j+1);
    S(j+2) = (  2*cos(z)*S(j+1) -   S(j) );
end

figure(3);  clf
subplot(2,1,1);  hold on
plot(0:N-1,Q,'mo',0:N-1,R,'k+')
axis([0 N -0.5 1.05])
title(['\bf legendre polynomials, P_j(' num2str(z) ')'])
subplot(2,1,2);  hold on
plot(0:N-1,cos((0:N-1)*z),'mo',0:N-1,S,'k+')
axis([0 N -0.5 1.05])
title(['\bf cos(j ' num2str(z) ') via 3-term recursion identity'])