function [zMatrix,nVec,mVec,jVec] = zStd(rho,theta,maxTerm)

% Evaluates Zernike polynomials over the points specified by col vectors
% rho and theta.
% rho = column vector of normalized radius
% theta = column vector of angles (rad)
% maxTerm = maximum number of standard Zernike polynomials
% Standard Zernike coefficients

points = length(rho);

for j = 1:maxTerm
    i = 0;  % index for polynomials
    n = 0; m = 0;

    while i<=j
        flag = 0;
        for m = 0:n
            if mod(n-m, 2)==0
                i = i+1;
                if i==j
                    flag = 1;
                    break   % out of for
                end
                if m~=0
                    i = i+1;
                    if i==j
                        flag = 1; 
                        break    % out of for
                    end
                end
            end
        end
        if flag == 1
%             [j m n]
            break   % out of while
        end
        n=n+1;
    end

    radialPoly = zeros(points,1);
    for s = 0:(n-m)/2
        coef = (-1)^s*factorial(n-s) / ...
            (factorial(s)*factorial((n+m)/2-s)*factorial((n-m)/2-s));
        power = n-2*s;
        radialPoly = radialPoly + coef * rho.^power;
    end
    if m==0
        zMatrix(:,j) = sqrt(n+1)*radialPoly;
        nVec(j,1) = n;
        mVec(j,1) = m; 
        jVec(j,1) = j;
    elseif mod(j,2)==0
        zMatrix(:,j) = sqrt(2)*sqrt(n+1)*radialPoly .* cos(m*theta);
        nVec(j,1) = n;
        mVec(j,1) = m;
        jVec(j,1) = j;
    else
        zMatrix(:,j) = sqrt(2)*sqrt(n+1)*radialPoly .* sin(m*theta);
        nVec(j,1) = n;
        mVec(j,1) = m;
        jVec(j,1) = j;
    end     
end