% Written by Ping Zhou
% September, 2007
% This function generates the anunular zernikes, which are orthonormal in
% an annular aperture.
% This annular zernikes have the same numbering as Durango software

function [zMatrix,nVec,elVec] = zAnnular(maxDegree,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);
i = 1;  % index for polynomials
for n = 0:maxDegree % degree of radial polynomial
    for el = n:-2:-n 
        m = abs(el); % angular order
        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 el>0
            zMatrix(:,i) = sqrt(2)*sqrt(n+1)*radialPoly .* cos(m*theta);
            nVec(i,1) = n;
            elVec(i,1) = el;
        elseif el==0
            zMatrix(:,i) = sqrt(n+1)*radialPoly;
            nVec(i,1) = n;
            elVec(i,1) = el;
        else
            zMatrix(:,i) = sqrt(2)*sqrt(n+1)*radialPoly .* sin(m*theta);
            nVec(i,1) = n;
            elVec(i,1) = el;
        end
        i = i+1;
    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