function design = zernike(n, rmax, radiusdata, angledata)
% ZERNIKE, returns an orthogonal polynomial design matrix up to 
%	zernike order n, by creating the Jacobi and Fourier polynomials
% 	and evaluating them inside this function.
% Inputs:
%	n 	   - maximum order of Zernike polynomial design to create
%	rmax	   - maximum radius desired (DUMMY)
%	radiusdata - column vector of radii data points
%	angledata  - column vector of angle data points
% Output:
%	design     - a design matrix of evaluated Zernike polynomials.
%
% Note:  Need to make error checking so that the data vectors are of same
%	 length.  Also need to create a way to make a uniform theta vector
% 	 from inside this function.
%
% Written by:  N. Locantore (12/5/97), modified (4/6/98)

if max(radiusdata > 1.05);
	disp ('#############################################')
	disp ('##  Radius has not been normalized to the  ##')
	disp ('##  unit disk. Results may be invalid.     ##')
	disp ('#############################################')
	end

design  = [];				% initialize design matrix
jdesign = [];				% initialize jacobi contribution	
zdesign = [];				% initialize fourier contribution

				%  loop to get design matrix
for i = 0:n;	
						%  for each order up to n
  k = 2*((i/2) - floor(i/2));	
  for j = k:2:i;				%  get design column (n,m)

     jdesignadd = polyval(jacobi(i,j), radiusdata);
     zdesigncos = cos(j*angledata);
     zdesignsin = sin(j*angledata);     

% add design columns

     if j ==0 ;			% Jacobi .* cosine only		
	zcoef = sqrt(i+1);		  %  Zernike coefficient		
	zaddcos = zcoef * zdesigncos;	  %  Fourier cosine
	zdesign = [zdesign, zaddcos];	  %  Fourier  design column (cosine)
	jdesign = [jdesign, jdesignadd];  %  Jacobi design column 
			
     else
	zcoef = sqrt(2*(i+1));		  %  Zernike coefficient
	zaddcos = zcoef * zdesigncos;	  %  Fourier cosine
	zaddsin = zcoef * zdesignsin;	  %  Fourier sine
	zdesign = [zdesign, zaddcos, zaddsin];	      % Fourier design columns 
	jdesign = [jdesign, jdesignadd, jdesignadd];  % Jacobi design columns

     end

  end

end

design = jdesign .* zdesign;	% Zernike design is element by element product
				% of the Jacobi design .* Fourier design

