function [Znm] = zernike(n,m,x,y,d)

// Ouput variables initialisation (not found in input variables)
Znm=[];

// Display mode
mode(0);

// Display warning for floating point exception
ieee(1);


//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// 
// Class:     Psych 221/EE 362
// Function:  zernike.m
// Author:    Patrick Maeda
// Purpose:   Compute Zernike Polynomial
// Date:      02.28.03
// 
// Matlab 6.1:  03.03.03
// 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// 
// zernike.m  is a function that computes the values of a Zernike Polynomial
//            over a circular pupil of diameter d
// 
// Output:
// Znm is the Zernike polynomial term of order n and frequency m
// 
// Input:
// n = highest power or order of the radial polynomial term, [a positive integer]
// m = azimuthal frequency of the sinusoidal component, [a signed integer, |m| <= n]
// x = 1-D array of pupil x-coordinate values, [length(x) must equal length(y)]
// y = 1-D array of pupil y-coordinate values, [length(y) must equal length(x)]
// d = pupil diameter
// 
// The Zernike Polynomial definitions used are taken from:
// Thibos, L., Applegate, R.A., Schweigerling, J.T., Webb, R., VSIA Standards Taskforce Members,
// """"Standards for Reporting the Optical Aberrations of Eyes""""
// OSA Trends in Optics and Photonics Vol. 35, Vision Science and its Applications,
// Lakshminarayanan,V. (ed) (Optical Society of America, Washington, DC 2000), pp: 232-244. 
// 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// initialize circular pupil function
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Imax = max(size(mtlb_double(x)));
Jmax = max(size(mtlb_double(y)));
a = d/2;

for I = 1:Imax //initialize circular pupil
 for J = 1:Jmax
   p(I,J) = mtlb_logic(sqrt(mtlb_a(mtlb_e(x,I)^2,mtlb_e(y,J)^2)),"<=",a);
 end;
end;

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Compute Normalization Constant
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nnm = sqrt((2*mtlb_a(n,1))/mtlb_a(1,mtlb_logic(m,"==",0)));

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Compute Zernike polynomial
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if mtlb_logic(n,"==",0) then
  Znm = p;
else
  Znm = zeros(Imax,Jmax);
  for I = 1:Imax
    for J = 1:Jmax
      r = sqrt(mtlb_a(mtlb_e(x,I)^2,mtlb_e(y,J)^2));
      if mtlb_logic(mtlb_e(x,I),">=",0) & mtlb_logic(mtlb_e(y,J),">=",0) | mtlb_logic(mtlb_e(x,I),">=",0) & mtlb_logic(mtlb_e(y,J),"<",0) then
        theta = atan(mtlb_e(y,J)/mtlb_a(mtlb_e(x,I),1.000000000D-20));
      else
        theta = mtlb_a(%pi,atan(mtlb_e(y,J)/mtlb_a(mtlb_e(x,I),1.000000000D-20)));
      end;
      for s = mtlb_imp(0,0.5*mtlb_s(n,abs(mtlb_double(m))))
      
        Znm(I,J) = mtlb_a(Znm(I,J),((((-1)^s)*factorial(mtlb_s(n,s)))*((r/a)^mtlb_s(n,2*s)))/((factorial(s)*factorial(mtlb_s(0.5*mtlb_a(n,abs(mtlb_double(m))),s)))*factorial(mtlb_s(0.5*mtlb_s(n,abs(mtlb_double(m))),s))));
      end;
      Znm(I,J) = ((p(I,J)*Nnm)*Znm(I,J))*mtlb_s(mtlb_logic(m,">=",0)*cos(m*theta),mtlb_logic(m,"<",0)*sin(m*theta));
    end;
  end;
end;

endfunction