function varargout=jacobi(l,m,x)
% P=JACOBI(l,m,x)
%
% Calculates the Jacobi polynomial P_l^{m,0}(x) by summation.
%
% INPUT:
%
% l      The degree                           [scalar]
% m      The first upper index, usually alpha [scalar]
% x      The abscissas                        [vector]
%
% OUTPUT:
%
% P      The Jacobi polynomial at x          [vector]
%
% EXAMPLE:
%
% jacobi('demo1') % Checks the equivalence with Legendre polynomials
% jacobi('demo2') % Checks the equivalence with hypergeometric functions
%
% Last modified by fjsimons-at-alum.mit.edu, 03/26/2009
% Last modified by dongwang-at-princeton.edu, 06/02/2008

% Should modify this to use the three-term recursion and return the whole
% set; this will be faster

% Supply default values for the arguments
defval('l',5)

if ~isstr(l)
  % Supply default values for the arguments
  defval('l',5)
  
  if l>64
    warning('High-degree Jacobi polynomials may be inaccurate')
  end
  
  defval('m',0)
  defval('x',linspace(0,1,20))
  % Make sure the argument is a row vector
  x=x(:)';
  
  % Initialize the matrix
  P=zeros(1,length(x));
  
  % Loop through the iteration
  for n=0:l
    % This according to De Villiers (2003) eq. (44) which appears to have
    % a typo - or, thus, ultimately, Szego - and see the demos.
    P=P+(x-1).^n.*(x+1).^(l-n)/...
      factorial(n+m)/factorial(l-n)^2/factorial(n);
    % This according to Abramowitz and Stegun p. 775, same thing
    % P=P+(x-1).^(l-n).*(x+1).^n/...
    %  factorial(n)/factorial(l+m-n)/factorial(l-n)/factorial(n);
  end
  % Apply the constants
  P=P*factorial(l+m)*factorial(l)/2^l;
  % Provide output
  vars={P};
  varargout=vars(1:nargout);
elseif strcmp(l,'demo1')
  % Calculate the Legendre functions
  x=linspace(0,1,20); m=0;
  for l=0:15
    difer(jacobi(l,m,x)-rindeks(legendre(l,x),m+1))
  end
elseif strcmp(l,'demo2')
  % Calculate the hypergeometric functions
  x=linspace(0,1,20); m=0;
  for l=0:15
    % This according to Mathworld eq. (25) which used to have a typo
    difer(jacobi(l,m,x)-factorial(l+m)/factorial(l)/factorial(m)*...
	  hypergeom([-l,l+m+1],m+1,1/2*(1-x)))
    % This according to Mathworld eq. (27) which was typo-free
    %    difer(jacobi(l,m,x)-factorial(l+m)/factorial(l)/factorial(m)*...
    %	  ((x+1)/2).^l.*hypergeom([-l,-l],m+1,(x-1)./(x+1)))
  end
end