function y = jacpoly(n,a,b,z)
% JACPOLY Calculate Jacobi polynomial
% Y = JACPOLY(N,A,B,Z) 
%   Calculates Jacobi polynomials of order N, type (A,B) at the points Z. 
%   Z must be a column vector.
%
% See also: djacpoly

% Peter Rucz
% 2010-2015

m = 0:n;
nz = length(z);
temp= repmat(binomial(n,m).*factorial(a+b+n+m)./factorial(a+m),nz,1) ...
    .*repmat((z-1)/2,1,n+1).^repmat(m,nz,1);
y = repmat(gamma(a+n+1)/(factorial(n)*gamma(a+b+n+1)),nz,1).*sum(temp,2);

end % of function jacpoly