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 = mtlb_double(d)/2;

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

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

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

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

if mtlb_logic(mtlb_double(n),"==",0) then
  Znm = p;
else
  Znm = zeros(Imax,Jmax);
  for I = 1:Imax
    for J = 1:Jmax
      r = sqrt(mtlb_a(mtlb_double(mtlb_e(x,I))^2,mtlb_double(mtlb_e(y,J))^2));
      if bool2s(bool2s(mtlb_logic(mtlb_double(mtlb_e(x,I)),">=",0))&bool2s(mtlb_logic(mtlb_double(mtlb_e(y,J)),">=",0)))|bool2s(bool2s(mtlb_logic(mtlb_double(mtlb_e(x,I)),">=",0))&bool2s(mtlb_logic(mtlb_double(mtlb_e(y,J)),"<",0))) then
        theta = atan(mtlb_double(mtlb_e(y,J))/mtlb_a(mtlb_double(mtlb_e(x,I)),1.000000000D-20));
      else
        theta = mtlb_a(%pi,atan(mtlb_double(mtlb_e(y,J))/mtlb_a(mtlb_double(mtlb_e(x,I)),1.000000000D-20)));
      end;
      for s = mtlb_imp(0,0.5*mtlb_s(mtlb_double(n),abs(mtlb_double(m))))
      
        // !! L.74: Matlab function factorial not yet converted, original calling sequence used
        // !! L.74: Matlab function factorial not yet converted, original calling sequence used
        // !! L.74: Matlab function factorial not yet converted, original calling sequence used
        // !! L.74: Matlab function factorial not yet converted, original calling sequence used
        Znm(I,J) = mtlb_a(Znm(I,J),((((-1)^s)*mtlb_double(factorial(mtlb_s(mtlb_double(n),s))))*((r/a)^mtlb_s(mtlb_double(n),2*s)))/((mtlb_double(factorial(s))*mtlb_double(factorial(mtlb_s(0.5*mtlb_a(mtlb_double(n),abs(mtlb_double(m))),s))))*mtlb_double(factorial(mtlb_s(0.5*mtlb_s(mtlb_double(n),abs(mtlb_double(m))),s)))));
      end;
      Znm(I,J) = ((p(I,J)*Nnm)*Znm(I,J))*mtlb_s(mtlb_logic(mtlb_double(m),">=",0)*cos(mtlb_double(m)*theta),mtlb_logic(mtlb_double(m),"<",0)*sin(mtlb_double(m)*theta));
    end;
  end;
end;
endfunction