function [design, ddesign, dddesign] = dradial(n, rmax, radiusdata, angledata)
% DRADIAL, returns first two radial derivative of an orthogonal polynomial 
%	design matrix up to zernike order n, by creating the Jacobi and 
%	Fourier polynomials, taking the Jacobi's derivatives, and then
%	evaluating the polynomial/frequency result.
% 	
% Inputs:
%	n 	   - maximum order of Zernike polynomial design to create
%	rmax	   - maximum radius to fit to the data
%	radiusdata - column vector of radii data points
%	angledata  - column vector of angle data points
% Output:
%	design	   - design matrix of evaluated Zernike polynomials.
%	ddesign	   - design matrix of evaluated Zernike polynomials at
%			their first radial derivative.
%	dddesign   - design matrix of evaluated Zernike polynomials at 
%		     	their second radial derivative.
%
% 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.
%
% Copyright (c)  N. Locantore 1998

% last edited 5/4/98

if max(radiusdata > 1);
	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	
djdesign = [];				% initialize jacobi contribution	
ddjdesign = [];				% initialize jacobi contribution	
zdesign = [];				% initialize fourier contribution

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

     	djacobival  = polyder(jacobi(i,j))/rmax;	
     	ddjacobival = polyder(polyder(jacobi(i,j)));
     	ddjacobival = ddjacobival / (rmax^2);	
     	jdesignadd  = polyval(jacobi(i,j), radiusdata);
     	djdesignadd  = polyval(djacobival, radiusdata);
     	ddjdesignadd  = polyval(ddjacobival, 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 
			djdesign = [djdesign, djdesignadd];  %  Jacobi design 1st d. 
			ddjdesign = [ddjdesign, ddjdesignadd];  %  Jacobi design 2nd d. 
			
     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
			djdesign = [djdesign, djdesignadd, djdesignadd];  % Jacobi design 1std
			ddjdesign = [ddjdesign, ddjdesignadd, ddjdesignadd]; % Jacobi des 2ndd
     end

  end

end

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

ddesign = djdesign .* zdesign;	
dddesign = ddjdesign .* zdesign;	

