function polyout = jacobi(n, m)
% JACOBI, returns the Jacobi polynomial for order n and 
%	  Fourier frequency m
% Inputs:
%	n - order of Jacobi polynomial
%	m - frequency of Fourier wave
% Output:
%	polyout - Jacobi polynomial in vector form, where the elements of the
%		  vector correspond to the coefficients for each degree
%		  coeff's are for [(n)th, (n-1)st, ... , 1st, 0th]
%
% Note:  To make completely general, need to error check that
%	 n >= m and that they are both > 0.
%		
% Note:  The construction of the Jacobi polynomials use equation #26 of
% 	 Appendix VII of "Principles of Optics", by Born and Wolf (1980). 
% 
% Note:  Function code by N. Locantore, 12/5/97.

rcoef = factoria((n-m)/2);	% creates the outer polynomial coefficient

deg1  = (n+m)/2;		% degree of first argument
poly1 = [1, zeros(1,deg1)];	% creates polynomial of degree deg1

deg2  = (n-m)/2;		% degree of second argument
polyh = [1, -1];		% holding the polynomial (x-1)
poly2 = [1];			% initializes a polynomial of degree 0

				% create second argument's polynomial 
for i = 1:deg2;
   poly2 = conv(poly2, polyh);
end

poly3 = conv(poly1, poly2);	% mulitplies both polynomials together

			% take corresponding % of derivatives of poly3
for j = 1:deg2;
   poly3 = polyder(poly3);
end

% now reconvert your polynomials to the correct degrees
% Note:  it is more convenient to do the algebra in the r^2 realm and 
% 	 double the resulting orders of coefficients

dummy = zeros(1, length(poly3));	% creates dummy vector of zeros
poly3 = [poly3; dummy];	        	% creates matrix for special binding
poly3 = poly3(:);			% makes vector from matrix
poly3 = poly3(1:length(poly3)-1);

outdegree = length(poly3) - m;		% correct degree of output Jacobi

polyout = poly3(1:outdegree);		% final polynomial
polyout = polyout / rcoef;		% multiplied by its constant
polyout = polyout';			% for easier viewing
 



