function [V_pq] = Zernike_polynomial(rho,theta, p,q)
%
% Note I could do these calculations in a vectorized manner but for time 
% constraints choose to just to the simplest thing and assume the input 
% arguments are scalars.
% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-04
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

assert( p>=0, 'p must be a nonnegative integer' );
assert( mod(p-abs(q),2)==0, 'must have the difference p - |q| even' );
assert( abs(q)<=p, 'incorrect sizes for p and q' );

% Evaluate the radial part 
sUpperLimt = 0.5 * ( p - abs(q) );
sArray     = 0:sUpperLimt;

num_fact_1 = (-1).^sArray;
num_fact_2 = factorial( p - sArray );
num_fact_3 = rho.^( p - 2 * sArray );

den_fact_1 = factorial( sArray );
den_fact_2 = factorial( 0.5*(p + abs(q)) - sArray );
den_fact_3 = factorial( 0.5*(p - abs(q)) - sArray );

num = num_fact_1 .* num_fact_2 .* num_fact_3;
den = den_fact_1 .* den_fact_2 .* den_fact_3;

R_pq = sum( num ./ den );

V_pq = R_pq * exp( i * q * theta );