function L = legendre(n,x,varargin)
% Based on the first recurrence relation found in 
% http://en.wikipedia.org/wiki/Associated_Legendre_function
x = x(:).';
L = zeros(n+1,length(x));
if (nargin == 2)
  normalization = "unnorm";
else
  normalization = varargin{1};
end
for m = 1:n+1
  if (strcmp(normalization,"norm"))
    scale = sqrt((n+0.5)/prod(n-m+2:n+m-1));
  elseif (strcmp(normalization,"sch"))
    scale = sqrt(2/prod(n-m+2:n+m-1));
  else
    scale = (-1).^(m-1);
  end
  LP = ones(n+1,length(x))*scale;
  LP(m,:) = prod(1:2:2*(m-1)-1)*LP(m,:).*sqrt(1-x.^2).^(m-1);
  if (m == n+1)
    L(m,:) = LP(m,:);
    break
  end  
  LP(m+1,:) = (2*m-1).*x.*LP(m,:);
  for l = m+1:n
    LP(l+1,:) = ((2*(l-1)+1).*x.*LP(l,:)-(l-1+m-1).*LP(l-1,:))./(l-m+1);
  end
  L(m,:) = LP(n+1,:);
end
if (strcmp(normalization,"sch"))
  L(1,:) = L(1,:)/sqrt(2/prod(n+1:n));
end