function z = zernpol(n,m,r,nflag)
5gV%jQgkC
%ZERNPOL Radial Zernike polynomials of order N and frequency M. S.bB.<
% Z = ZERNPOL(N,M,R) returns the radial Zernike polynomials of 0Bx.jx0?
% order N and frequency M, evaluated at R. N is a vector of ad). X:Qs
% positive integers (including 0), and M is a vector with the tl |Qw";I
% same number of elements as N. Each element k of M must be a Dz4fP;n
% positive integer, with possible values M(k) = 0,2,4,...,N(k) ALqP;/
% for N(k) even, and M(k) = 1,3,5,...,N(k) for N(k) odd. R is \Lxsg!wtJ
% a vector of numbers between 0 and 1. The output Z is a matrix w{J0K;L
% with one column for every (N,M) pair, and one row for every -MrEJ
% element in R. P>/n!1c
% 0p\cDrB?
% Z = ZERNPOL(N,M,R,'norm') returns the normalized Zernike poly- 6mr5`5~w
% nomials. The normalization factor Nnm = sqrt(2*(n+1)) is 1=x4m=wV
% chosen so that the integral of (r * [Znm(r)]^2) from r=0 to /xmUu0H$R
% r=1 is unity. For the non-normalized polynomials, Znm(r=1)=1 v%|^\A"V
% for all [n,m]. XOQj?Q7)U
% 3+gp_7L
% The radial Zernike polynomials are the radial portion of the &h.E
B
% Zernike functions, which are an orthogonal basis on the unit KS($S(Fi
% circle. The series representation of the radial Zernike &u-H/CU%
% polynomials is okx~F9
% =k#SQ/@
% (n-m)/2 +;7Rz_.6f
% __ [bd fp
a
% m \ s n-2s w)o^?9T
% Z(r) = /__ (-1) [(n-s)!/(s!((n-m)/2-s)!((n+m)/2-s)!)] * r jX5lwP
Q|F
% n s=0 @G0k+
% `a}!t=~#w
% The following table shows the first 12 polynomials. l$$N~F N
% b&BSigrvou
% n m Zernike polynomial Normalization f!;4-.p`
% --------------------------------------------- RkVU^N"
% 0 0 1 sqrt(2) &D,gKT~
% 1 1 r 2 "V!y"yQ
% 2 0 2*r^2 - 1 sqrt(6) rWKc,A[
% 2 2 r^2 sqrt(6) zG|}| //}
% 3 1 3*r^3 - 2*r sqrt(8) ;W6P$@'zs
% 3 3 r^3 sqrt(8) f(Q-W6
% 4 0 6*r^4 - 6*r^2 + 1 sqrt(10) 4R5+"h:
% 4 2 4*r^4 - 3*r^2 sqrt(10) M#II,z>q
% 4 4 r^4 sqrt(10) trL:qD+{(
% 5 1 10*r^5 - 12*r^3 + 3*r sqrt(12) GQXN1R
% 5 3 5*r^5 - 4*r^3 sqrt(12) z}v6!u|iZu
% 5 5 r^5 sqrt(12) "Gx(-NH+
% --------------------------------------------- X(F2 5
% ` z<k7ig
% Example: ]J\tosTi
% wjGD[~mB
% % Display three example Zernike radial polynomials $
BV4 i$
% r = 0:0.01:1; _w8iPL5:
% n = [3 2 5]; ]YcM45xg
% m = [1 2 1]; 6;g_}Zx
% z = zernpol(n,m,r); SXn\k;F<
% figure .Ua|KKK C
% plot(r,z) P#5&D*`}h
% grid on sqw^Hwy=!2
% legend('Z_3^1(r)','Z_2^2(r)','Z_5^1(r)','Location','NorthWest') cx?t C#t
% HMT^gmF)
% See also ZERNFUN, ZERNFUN2. 3*9<JHu
aW-'Jg=@H^
% A note on the algorithm. 5}FPqyK"
% ------------------------ 1Wzm51RU
% The radial Zernike polynomials are computed using the series
r=YprVX
% representation shown in the Help section above. For many special a=3?hVpB
% functions, direct evaluation using the series representation can JAM4
R_
% produce poor numerical results (floating point errors), because u!TVvc
% the summation often involves computing small differences between ||TKo967]
% large successive terms in the series. (In such cases, the functions ?k)(~Y&@p
% are often evaluated using alternative methods such as recurrence SB
R=
% relations: see the Legendre functions, for example). For the Zernike \Ub=Wm\
% polynomials, however, this problem does not arise, because the uY+N163i
% polynomials are evaluated over the finite domain r = (0,1), and ^Wk.D-
% because the coefficients for a given polynomial are generally all "|&xUWJ!)
% of similar magnitude. s,UccA@
% ^W-03
% ZERNPOL has been written using a vectorized implementation: multiple [ix45xu7
% Zernike polynomials can be computed (i.e., multiple sets of [N,M] M$j]VZ
% values can be passed as inputs) for a vector of points R. To achieve ajFSbi)l
% this vectorization most efficiently, the algorithm in ZERNPOL S~auwY ,<
% involves pre-determining all the powers p of R that are required to S'"(zc3=
% compute the outputs, and then compiling the {R^p} into a single 7XLz Ewa
% matrix. This avoids any redundant computation of the R^p, and 5yO%| )
% minimizes the sizes of certain intermediate variables. QF 2Eg
% sr(f9Vl
% Paul Fricker 11/13/2006 C"w>U
,<]X0;~oB
|ho|Kl `=
% Check and prepare the inputs: ao>`[-
% ----------------------------- K1c@]]y)
if ( ~any(size(n)==1) ) || ( ~any(size(m)==1) ) <a_Q1 l
error('zernpol:NMvectors','N and M must be vectors.') Y'6GY*dL
end *gHGi(U(U
OEc$ro=m*
if length(n)~=length(m) G
@ib
error('zernpol:NMlength','N and M must be the same length.') p7y8/m\6
end |1vikG8
J A'C\
n = n(:); =#qf0
m = m(:); qH(3Z^ #.|
length_n = length(n); H% c:f
"_Wv,CYmNr
if any(mod(n-m,2)) XBi}hT
error('zernpol:NMmultiplesof2','All N and M must differ by multiples of 2 (including 0).') '{9nQDgT
end z)$X/v
v{7Jzjd
if any(m<0) Cf#[E~2 4
error('zernpol:Mpositive','All M must be positive.') `em}vdY
end J)R;NYl
>gNVL
(
if any(m>n) R<>ptwy
error('zernpol:MlessthanN','Each M must be less than or equal to its corresponding N.') rOC2 S(m
end f,utA3[
fz
H$`X'M
if any( r>1 | r<0 ) *T(z4RVg
error('zernpol:Rlessthan1','All R must be between 0 and 1.') sBozz #
end NijvFT$V1
0'ha!4h3Z
if ~any(size(r)==1) gc|?$aE
error('zernpol:Rvector','R must be a vector.') +m Plid\
end 4\*!]5i
`/en&l
r = r(:); Y_qRW. k
length_r = length(r); WJ)( *1
\bJ,8J1C
if nargin==4 >U/m/H'
isnorm = ischar(nflag) & strcmpi(nflag,'norm'); 58FjzW
if ~isnorm ,A`.u \f(:
error('zernpol:normalization','Unrecognized normalization flag.') Un{hI`3]
end !F3Y7R
else
%W!C
isnorm = false; 4c"x&x|
end ~@8r-[
t#pF.!9=
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ,'~8{,h5
% Compute the Zernike Polynomials C(( 7
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ROZOX$XM
8*O]
% Determine the required powers of r: 2u0C~s
% ----------------------------------- i|zs
Li/
rpowers = []; |TCHPKN
for j = 1:length(n) QH:PClW![
rpowers = [rpowers m(j):2:n(j)]; -*;-T9
end Rlvb@aXgy
rpowers = unique(rpowers); o&tETJ5Bhe
b(<#n6a}\
% Pre-compute the values of r raised to the required powers, H=2sT +Sp
% and compile them in a matrix: iwJeV J
% ----------------------------- f|eUpf%)
if rpowers(1)==0 2%0J/]n\A"
rpowern = arrayfun(@(p)r.^p,rpowers(2:end),'UniformOutput',false); o[C,fh,$
rpowern = cat(2,rpowern{:}); #:E}Eby/6I
rpowern = [ones(length_r,1) rpowern]; ~";GH20
else G$b*N4yR
rpowern = arrayfun(@(p)r.^p,rpowers,'UniformOutput',false); Gn}G$uk61
rpowern = cat(2,rpowern{:}); ^HpUbZpat)
end {9(#X]'
pwq a/Yi
% Compute the values of the polynomials: G&P[n8Z$
% -------------------------------------- <5oG[1j
z = zeros(length_r,length_n); ')ZM#
:G
for j = 1:length_n tqdw
y.
s = 0:(n(j)-m(j))/2; )n8(U%q$
pows = n(j):-2:m(j); }u"iA^'Ot
for k = length(s):-1:1 FjUf|
p = (1-2*mod(s(k),2))* ... (8bo"{zI
prod(2:(n(j)-s(k)))/ ... C'4gve 7!
prod(2:s(k))/ ... Y",
:u@R
prod(2:((n(j)-m(j))/2-s(k)))/ ... ["N{6d&Q
prod(2:((n(j)+m(j))/2-s(k))); u} y)'eH
idx = (pows(k)==rpowers); eBw6k09C+
z(:,j) = z(:,j) + p*rpowern(:,idx); xWNB/{F
end 9 F"2$;
J!l/!Z>!cF
if isnorm h_O6Z2J1
z(:,j) = z(:,j)*sqrt(2*(n(j)+1)); +k@$C,A
end `/WX!4eR,
end NWK+.{s>m
'`.bmiM
% EOF zernpol