%
%  legendre kernel -- djm, 04 july 2005
%

N = 100;

x = -1:0.002:1.0;  y = interp1(x,x,0.2,'nearest');

%  find zero denominator
denom = x-y;  pt = find(denom==0);  denom(pt) = 1;
denom = 1./denom;  denom(pt) = 0;

PNpy = legendre(N+1,y);  PNpy = PNpy(1);
PN0y = legendre(N  ,y);  PN0y = PN0y(1);
PNmy = legendre(N-1,y);  PNmy = PNmy(1);

PNpx = legendre(N+1,x);  PNpx = PNpx(1,:);
PN0x = legendre(N  ,x);  PN0x = PN0x(1,:);
PNmx = legendre(N-1,x);  PNmx = PNmx(1,:);

%PN0xp = N*(x.*PN0x-PNmx)./(x.^2 -1);

kernel = (N+1)/2 * (PN0y*PNpx - PNpy*PN0x) .* denom;

%  l'Hopital rule
kernel(pt) = (N+1)/2 * (PN0y*(N+1)*(y*PNpy-PN0y)-PNpy*N*(y*PN0y-PNmy))/(y^2 - 1);
%kernel(pt) = 0;

figure(1);  clf;  hold on
plot(x,kernel,'k.')
plot(x,kernel)
plot(x(pt),kernel(pt),'ro')
title('\bf fourier-legendre kernel')
ylabel(['\bf y = ' num2str(y)])
xlabel('\bf x-axis')

%figure(2);  clf;  hold on
%plot(x,PN0x,'r',x,PN0xp,'b')