function z = zernpol(n,m,r,nflag) d\]Yk]r
%ZERNPOL Radial Zernike polynomials of order N and frequency M. wSEWwU[
% Z = ZERNPOL(N,M,R) returns the radial Zernike polynomials of %<0eA`F4
% order N and frequency M, evaluated at R. N is a vector of 7fHc[,
% positive integers (including 0), and M is a vector with the "%qzj93>
% same number of elements as N. Each element k of M must be a 5T@aCC@$h
% positive integer, with possible values M(k) = 0,2,4,...,N(k) "4o=,$E=
% for N(k) even, and M(k) = 1,3,5,...,N(k) for N(k) odd. R is A1"SLFY
% a vector of numbers between 0 and 1. The output Z is a matrix JH8}Ru%Z
% with one column for every (N,M) pair, and one row for every `=UWqb(K_
% element in R. a5YIUVCv
% _oG&OJ@
% Z = ZERNPOL(N,M,R,'norm') returns the normalized Zernike poly- FAsFjRS
% nomials. The normalization factor Nnm = sqrt(2*(n+1)) is W,XTF
% chosen so that the integral of (r * [Znm(r)]^2) from r=0 to bD^b
% r=1 is unity. For the non-normalized polynomials, Znm(r=1)=1 wE]K~y!`
% for all [n,m]. ,m_WR7!$E
% Hnk:K9u.B:
% The radial Zernike polynomials are the radial portion of the X5LBEOG
% Zernike functions, which are an orthogonal basis on the unit lf(`SYQnOY
% circle. The series representation of the radial Zernike 6eFp8bANN#
% polynomials is (o5j'2:.
% qpIC{'A.
% (n-m)/2 }e2VY
% __ M{g%cR0
% m \ s n-2s `d7n?|pD
% Z(r) = /__ (-1) [(n-s)!/(s!((n-m)/2-s)!((n+m)/2-s)!)] * r aJ Z"D8C
% n s=0 V!v:]E
% b~=0[Rv
% The following table shows the first 12 polynomials. d.UQW
yLG
% NIdZ
% n m Zernike polynomial Normalization WOzf]3Xcj
% --------------------------------------------- [AMAa]^
% 0 0 1 sqrt(2) >NYW{(j
% 1 1 r 2 U_x )#,4
% 2 0 2*r^2 - 1 sqrt(6) BTgG4F/)
% 2 2 r^2 sqrt(6) I[)% , jd
% 3 1 3*r^3 - 2*r sqrt(8) Wbr+KX8)
% 3 3 r^3 sqrt(8) W-NDBP:
% 4 0 6*r^4 - 6*r^2 + 1 sqrt(10) Q></`QWpoB
% 4 2 4*r^4 - 3*r^2 sqrt(10) o1M$.*
% 4 4 r^4 sqrt(10) J'$>Gk]
% 5 1 10*r^5 - 12*r^3 + 3*r sqrt(12) A,sr[Pa@
% 5 3 5*r^5 - 4*r^3 sqrt(12) ry9T U
% 5 5 r^5 sqrt(12) /xtq_*I1S
% --------------------------------------------- .X# `k
% hn#1%p6t
% Example: y;_% W
% n]E?3UGD@W
% % Display three example Zernike radial polynomials ,]bB9tid
% r = 0:0.01:1; zR2B-
&]H
% n = [3 2 5]; aDO!
% m = [1 2 1]; T5a*z}L5
% z = zernpol(n,m,r); '>r7V
% figure i3rH'B-I.
% plot(r,z) Fu!RhsW5j
% grid on UQ{L{H
% legend('Z_3^1(r)','Z_2^2(r)','Z_5^1(r)','Location','NorthWest') *98$dQR$
% ^0/j0]O
% See also ZERNFUN, ZERNFUN2. 9eH$XYy
0u\GO;
% A note on the algorithm. feQ **wI
% ------------------------ JvY}-}?c
% The radial Zernike polynomials are computed using the series dqN5]Sb2B
% representation shown in the Help section above. For many special ~l)-wNqR4r
% functions, direct evaluation using the series representation can FCnm1x#
% produce poor numerical results (floating point errors), because Lhqz\ o
% the summation often involves computing small differences between @Y1s$,=xB
% large successive terms in the series. (In such cases, the functions C=eF.FB;'
% are often evaluated using alternative methods such as recurrence r-:Uz\gM
% relations: see the Legendre functions, for example). For the Zernike P9T}S
% polynomials, however, this problem does not arise, because the r zt Ru
% polynomials are evaluated over the finite domain r = (0,1), and l<A|d{" ]
% because the coefficients for a given polynomial are generally all oH_;4QU4y
% of similar magnitude. |UX(+;n
% ax (c#
% ZERNPOL has been written using a vectorized implementation: multiple 2 B
% Zernike polynomials can be computed (i.e., multiple sets of [N,M] *s:(jDlv
% values can be passed as inputs) for a vector of points R. To achieve 6_FE 4RR[
% this vectorization most efficiently, the algorithm in ZERNPOL YJ\Xj56gv
% involves pre-determining all the powers p of R that are required to /pT=0=
% compute the outputs, and then compiling the {R^p} into a single Ymg|4%O@
% matrix. This avoids any redundant computation of the R^p, and 0N$v"uX@
% minimizes the sizes of certain intermediate variables. 1,% R;7J=g
% Y\No4w ^|d
% Paul Fricker 11/13/2006 b45-:mi!
NxsBX:XDn
>U[j]V]
% Check and prepare the inputs: ~vyf4TF<#
% ----------------------------- d8c=L8~jt
if ( ~any(size(n)==1) ) || ( ~any(size(m)==1) ) 5yP\I+Fm
error('zernpol:NMvectors','N and M must be vectors.') B!?%O
end * -8&[D0
*2,VyY
if length(n)~=length(m) AJF#Aw `o
error('zernpol:NMlength','N and M must be the same length.') m)'=G%y
end LGw$v[wb
Y2~nBb
n = n(:); OG>}M$Ora
m = m(:); OWg(#pZk
length_n = length(n); <nT
+$
cWe"%I
if any(mod(n-m,2)) 5Ou`z5S\k
error('zernpol:NMmultiplesof2','All N and M must differ by multiples of 2 (including 0).') %5"9</a&G
end YwjKAyLU
Na]:_K5Dp
if any(m<0) ) QU
error('zernpol:Mpositive','All M must be positive.') <+?7H\b
end GkQpELO:
J;8IY=
if any(m>n) 0"28'
error('zernpol:MlessthanN','Each M must be less than or equal to its corresponding N.') j~[z2tV
end jK& h~)
o <pf#tifv
if any( r>1 | r<0 ) :&V h?
error('zernpol:Rlessthan1','All R must be between 0 and 1.') "hH.#5j
end /rnu<Q#iH
j i7[nY
if ~any(size(r)==1) !Fd~~v
error('zernpol:Rvector','R must be a vector.') Z8K?
end wsI`fO^A8
61k"p2?+
r = r(:); Ku5\]
length_r = length(r); \Je0CD=e`
|B.Y6L6l
if nargin==4 ) l:[^$=,
isnorm = ischar(nflag) & strcmpi(nflag,'norm'); =5~jx
if ~isnorm `ZyI!"
error('zernpol:normalization','Unrecognized normalization flag.') (MxQ+D\
end ,St#Vla
else @tNz Q8
isnorm = false; "n(hfz0y%
end S2sQOM@
hK L4cpK4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !Qu"BF
% Compute the Zernike Polynomials ib#KpEk
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R3gdLa.
r*2+xDoEi
% Determine the required powers of r: L*xhGoC=
% ----------------------------------- x[~b2o
rpowers = []; [+T.at
for j = 1:length(n) vWzm@
rpowers = [rpowers m(j):2:n(j)]; tkeoNuAM
end %[Wh [zZy
rpowers = unique(rpowers); CkOz
&^ 1$^=
% Pre-compute the values of r raised to the required powers, N} G[7Rp8l
% and compile them in a matrix: AG`L64B
% ----------------------------- y\4L{GlBM
if rpowers(1)==0 46_xyz3+
rpowern = arrayfun(@(p)r.^p,rpowers(2:end),'UniformOutput',false); LDgrR[
rpowern = cat(2,rpowern{:}); X(
\AB
rpowern = [ones(length_r,1) rpowern]; LM~[@_j
else ]!?;@$wx
rpowern = arrayfun(@(p)r.^p,rpowers,'UniformOutput',false); 4;e5H_}Oo
rpowern = cat(2,rpowern{:}); md)c0Bg8~
end ]7W&JKmA&
!;q&NHco
% Compute the values of the polynomials: {f3YsM;]C
% -------------------------------------- at#ja_ hd
z = zeros(length_r,length_n); 0L2 F[TN
for j = 1:length_n )WNzWUfn=z
s = 0:(n(j)-m(j))/2; HSjlD{R
pows = n(j):-2:m(j); LO9=xGj.
for k = length(s):-1:1 ?GKb7Oj
p = (1-2*mod(s(k),2))* ... 7Wf/$vRab
prod(2:(n(j)-s(k)))/ ... 6M @[B|Q(
prod(2:s(k))/ ... [M]
prod(2:((n(j)-m(j))/2-s(k)))/ ... {f/~1G[M
prod(2:((n(j)+m(j))/2-s(k))); I667Gz$j5
idx = (pows(k)==rpowers); H%qsjB^
z(:,j) = z(:,j) + p*rpowern(:,idx); F~R;n_IJ
end q6McG HT
`uv2H$
if isnorm R1adWBD>
z(:,j) = z(:,j)*sqrt(2*(n(j)+1)); Q|S.R1L^
end B3pCy~*5
end ZQ4p(6a
Z3"%`*Tmq-
% EOF zernpol