function S = ZernikeAnnularEval( Z, x, y, e )

[m1, n1] = size( x );
[m2, n2] = size( y );
if( m1 ~= m2 || n1 ~= n2 )
    error( 'x & y are not the same size' );
end

if( length( Z ) > 22 ) 
    disp('ZernikeAnnularEval() is not implemented with >22 terms');
    return;
elseif length(Z)<22
    Z(22)=0;
end

r2 = x.*x + y.*y;
r = sqrt( r2 );
r3 = r2 .* r;
r4 = r2 .* r2;
r5 = r3 .* r2;
r6 = r3 .* r3;

t = atan2( y, x );
s = sin( t );
c = cos( t );
s2 = sin( 2*t );
c2 = cos( 2*t );
s3 = sin( 3*t );
c3 = cos( 3*t );
s4 = sin( 4*t );
c4 = cos( 4*t );
s5 = sin( 5*t );
c5 = cos( 5*t );

e2=e*e;
e4=e2*e2;
e6=e4*e2;
e8=e6*e2;
e10=e8*e2;
e12=e10*e2;

S =     Z(1)  * (1 + 0*x); % 0*x to set NaNs properly

den=sqrt(1+e2);
S = S + Z(2)  * 2 * r .* c/den;
S = S + Z(3)  * 2 * r .* s/den;

den=1-e.^2;
S = S + Z(4)  * sqrt( 3 )  * (2*r2 - 1-e2)/den;

den=sqrt(1+e2+e4);
S = S + Z(5)  * sqrt( 6 )  * r2 .* s2/den;
S = S + Z(6)  * sqrt( 6 )  * r2 .* c2/den;

den=sqrt( (1-e2)^2*(1+e2)*(1+4*e2+e4) );
S = S + Z(7)  * sqrt( 8 )  * (3*r3 - 2*r -2*e4*r + e2*r.*(3*r2-2) ) .* s/den;
S = S + Z(8)  * sqrt( 8 )  * (3*r3 - 2*r -2*e4*r + e2*r.*(3*r2-2) ) .* c/den;

den=sqrt(1+e2+e4+e6);
S = S + Z(9)  * sqrt( 8 )  * r3.* s3/den;
S = S + Z(10) * sqrt( 8 )  * r3.* c3/den;

den=(1-e2)^2;
S = S + Z(11) * sqrt( 5 )  * (6*r4 - 6*r2 + 1+e4+e2*(4-6*r2) )/den;

den=(1-e2)^3*(1+e2+e4);
num=sqrt((1-e2)^4*(1+e2+e4)/(1+4*e2+10*e4+4*e6+e8));
S = S + Z(12) * sqrt( 10 ) * (4*r4 - 3*r2-3*e6*r2-e2*r2.*(3-4*r2)-e4*r2.*(3-4*r2)) .* c2*num/den;
S = S + Z(13) * sqrt( 10 ) * (4*r4 - 3*r2-3*e6*r2-e2*r2.*(3-4*r2)-e4*r2.*(3-4*r2)) .* s2*num/den;

den=sqrt(1+e2+e4+e6+e8);
S = S + Z(14) * sqrt( 10 ) * r4 .* c4/den;
S = S + Z(15) * sqrt( 10 ) * r4 .* s4/den;

den=(1-e2)^3*(1+4*e2+e4);
numE=sqrt((1-e2)^2*(1+4*e2+e4)/(1+9*e2+9*e4+e6) );
numQ=10*r5-12*r3 +3*r+3*e8*r-12*e6*r.*(r2-1)+2*e4*r.*(15-24*r2+5*r4)+4*e2*r.*(3-12*r2+10*r4);
S = S + Z(16) * sqrt( 12 ) * numE*numQ.* c/den;
S = S + Z(17) * sqrt( 12 ) * numE*numQ.* s/den;

den=(1-e2)^4*(1+e2)*(1+e4);
numE=sqrt((1-e2)^6*(1+e2)*(1+e4)/(1+4*e2+10*e4+20*e6+10*e8+4*e10+e12));
numQ=r3.*(5*r2-4-4*e8-e2*(4-5*r2)-e4*(4-5*r2)-e6*(4-5*r2));
S = S + Z(18) * sqrt( 12 ) * numE*numQ.* c3/den;
S = S + Z(19) * sqrt( 12 ) * numE*numQ.* s3/den;

den=sqrt(1+e2+e4+e6+e8+e10);
S = S + Z(20) * sqrt( 12 ) * r5 .* c5/den;
S = S + Z(21) * sqrt( 12 ) * r5 .* s5/den;

den=(1-e2)^3;
S = S + Z(22) * sqrt(  7 ) * (20*r6 - 30*r4 + 12*r2 - 1-e6+3*e4*(-3+4*r2)-3*e2*(3-12*r2+10*r4))/den;


