%This code was written to deal with "Zernike polynomials" code graciously %donated by Paul Fricker via file exchange. %Here you will find a practical example of a function decomposition by %Zernike basis. %The function is F below, feel free to modify %Unlike Paul's example found in 'zernfun2.m' here the domain is the true %unit circle, without NaN's filling it up to the unit square. %Alex Chtchetinin, December 10, 2007 N=100;%dimension of theta-R grid %r=linspace(1/N/2,1-1/N/2,N);t=linspace(1/N/2,2*pi*(1-1/N/2),N); r=linspace(0,1,N);t=linspace(0,2*pi,N); [t,r]=meshgrid(t,r); [x,y] = pol2cart(t,r); F=x.^3-y.^2;%Function to be decomposed figure,surface(x,y,F);shading interp title('Function to be decomposed');axis square p = 0:8;%indexes of Zernike functions to be used in the decomposition n=length(p); z = zernfun2(p,r(:),t(:),'norm'); zz=reshape(z,N,N,length(p)); figure; for k = 1:n z = zz(:,:,k); subplot(ceil(n^.5),ceil(n^.5),k) surface(x,y,z), shading interp set(gca,'XTick',[],'YTick',[]) axis square title(['Z_' num2str(p(k))]); end for i=1:n zzz=zz(:,:,i).*F.*r; G = @(x,y)interp2(t,r,zzz,x,y); Q(i) = dblquad(G,0,2*pi,0,1); end disp(Q) Z=0*zz(:,:,1); figure; for k = 1:n Z=Z+Q(k)*zz(:,:,k); subplot(ceil(n^.5),ceil(n^.5),k) surface(x,y,Z), shading interp set(gca,'XTick',[],'YTick',[]) axis square if k==1 title(['Z_' num2str(p(k))]); else title(['Z_0 to Z_' num2str(p(k))]); end end figure,surface(x,y,Z-F);shading interp; title('Residual error');axis square