% Written by Ping Zhou
% September, 2007
% This function generates the anunular zernikes with Zernike Standard
% Polynomial numbering.

function [zMatrix,nVec,mVec] = zAnnularStd(maxTerm,rho,theta,e)

% rho = column vector of normalized radius
% theta = column vector of angles (rad)
% maxDegree = maximum degree of radial polynomials. All polynomials through
% maxDegree are evaluated.
% e = obscuration ratio
% Reference: 
% 1. Mahajan, V.N., Zernike annular polynomials for imaging systems
% with annular pupils, J. Opt. Soc. Am., 71, 75, 1981
% 2. Brett Patterson, Zernike polynomials for circular and anuualar domains

points = length(rho);
for a = 1:maxTerm
    i = 0;  % index for polynomials
    n = 0; m = 0;
    a
    while i<=a
        flag = 0;
        for m = 0:n
            if mod(n-m, 2)==0
                i = i+1;
                if i==a
                    flag = 1;
                    break   % out of for
                end
                if m~=0
                    i = i+1;
                    if i==a
                        flag = 1;
                        break    % out of for
                    end
                end
            end
        end
        if flag == 1
%                         [a m n]
            break   % out of while
        end
        n=n+1;
    end

    radialPoly = zeros(points,1);

    if m==0 && rem(n,2)==0
        radialPoly = R(0,n,sqrt((rho.^2-e^2)/(1-e^2)));
    elseif m==n
        epsilon = 0;
        for p = 0:n
            epsilon = epsilon + e^(2*p);
        end
        radialPoly = rho.^n/sqrt(epsilon);
    else
        j = (n-m)/2;
        u = rho.^2;
        radialPoly = sqrt((1-e^2)/2/(2*j+m+1)/h(m,j,e)).*rho.^m.*Q(m,j,u,e);
    end

    radialPoly(rho<e) = NaN;

    if m==0
        zMatrix(:,a) = sqrt(n+1)*radialPoly;
        nVec(a,1) = n;
        mVec(a,1) = m;
 
    elseif mod(a,2)==0
        zMatrix(:,a) = sqrt(2)*sqrt(n+1)*radialPoly .* cos(m*theta);
        nVec(a,1) = n;
        mVec(a,1) = m;
     
    else
        zMatrix(:,a) = sqrt(2)*sqrt(n+1)*radialPoly .* sin(m*theta);
        nVec(a,1) = n;
        mVec(a,1) = m;
    
    end
end

%%%%%%%%%%%%% Function of Q %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Output = Q(m,j,u,e)
% return a vector if rho is a vector
points = length(u);
Output = zeros(points,1);
if m==0
    Output = R(0,2*j,sqrt((u-e^2)/(1-e^2)));
else
    Coef = 2*(2*j+2*m-1)*h(m-1,j,e)/(j+m)/(1-e^2)./Q(m-1,j,0,e);
    for i = 0:j
        Output = Output + Q(m-1,i,0,e).*Q(m-1,i,u,e)./h(m-1,i,e);
    end
    Output = Output*Coef;
end

%%%%%%%%%%%% Function of R %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function radialPoly = R(m,n,rho)
% return a vector if rho is a vector
points = length(rho);
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


%%%%%%%%%%%% Function of h %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function output = h(m,j,e)
% always return a number   

if m==0
    output = (1-e^2)/2/(2*j+1);
else
    output = -2*(2*j+2*m-1)*Q(m-1,j+1,0,e)/(j+m)/(1-e^2)/Q(m-1,j,0,e)*h(m-1,j,e);
end