function dddesign = ddzernike(n, rmax, radiusdata, angledata)
% DDZERNIKE, returns second 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:
%	dddesign   - a design matrix of evaluated Zernike polynomials at 
%		     their first 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.
%
% Written by:  N. Locantore (1/5/98), edited (1/19/98)

% from the /matlab directory

dddesign = [];				% 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)
     ddjacobival = polyder(polyder(jacobi(i,j)));
     ddjacobival = ddjacobival / (rmax^2);	
     jdesignadd = 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 
			
     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

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

