function fdir = fndir(f,direction)
%FNDIR Directional derivative of a function.
%
%   FNDIR(F,DIRECTION)  returns the (ppform of the) derivative in the given 
%   DIRECTION(s) of the function contained in F (and of the same order).
%   If the function in F is d-variate, then DIRECTION must be a (list of)
%   d-vector(s), i.e., of size [d,nd] for some nd. 
%
%   Assuming the function in F to be d-variate and m-vector-valued, FDIR 
%   describes the (m*nd)-valued function whose value at a point, reshaped 
%   as a matrix of size [m,nd], provides in its j-th column the directional 
%   derivative, of the function in F, in the direction DIRECTION(:,j), j=1:nd.
%
%   This command does not work for rational splines; for them, use FNTLR 
%   instead.
%
%   For example, if  f  describes a d-variate m-vector-valued function and 
%    x  is some point in its domain, then
%  
%      reshape(fnval(fndir(f,eye(d)),x),m,d)
%
%   is the Jacobian of that function at that point. 
%   As a related example, the next statements plot the gradients of (a good 
%   approximation to) the Franke function at a regular mesh: 
%
%      xx = linspace(-.1,1.1,13); yy = linspace(0,1,11);
%      [x,y] = ndgrid(xx,yy); z = franke(x,y);
%      pp2dir = fndir(csapi({xx,yy},z),eye(2));
%      grads = reshape(fnval(pp2dir,[x(:) y(:)].'),[2,length(xx),length(yy)]);
%      quiver(x,y,squeeze(grads(1,:,:)),squeeze(grads(2,:,:)))
%
%   See also FNDER, FNTLR, FNINT.

%   C. de Boor 25apr98
%   cb : 16may98 (don't reject B-form, convert to ppform instead)
%   cb : 23jan00 (correct help (dn --> nd); refer to FNTLR for rational splines)
%   cb : 10apr00 (introduce optional third argument for multidimensional case)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.7 $  $Date: 2000/05/09 02:26:49 $

if nargin<2, error('The DIRECTION(S) must be specified.'), end

if ~isstruct(f), f = fn2fm(f); end

if f.form(1)=='B', f = fn2fm(f,'pp'); end

switch f.form(1)
case 'p' % the function is in ppform:
   [breaks,coefs,l,k,d] = ppbrk(f);
   if iscell(breaks) % the function is multivariate
      m = length(k);
      dirsize = size(direction);
      if dirsize(1)~=m
         error(['DIRECTION(s) should be  ' num2str(m) '-vector(s).'])
      end
      sizec = [d,l.*k]; %size(coefs); 
      fdirc = zeros([sizec,dirsize(2)]);
      for i=m:-1:1
         dd = prod(sizec(1:m));
         if any(direction(i,:))~=0
            dpp = fndirp(ppmak(breaks{i},reshape(coefs,dd*l(i),k(i)),dd));
            %fdirc = fdirc + direction(i)*reshape(dpp.coefs,sizec);
            fdirc = fdirc + ...
                    reshape(dpp.coefs(:)*direction(i,:),[sizec,dirsize(2)]);
         end
         if m>1
            sizec(2:m+1) = sizec([m+1,2:m]);
            coefs = permute(coefs,[1,m+1,2:m]);
            fdirc = permute(fdirc,[1,m+1,2:m,m+2]);
         end
      end
      sizec(1) = d*dirsize(2);
      if m>1&dirsize(2)>1
         fdir = ppmak(breaks,reshape(permute(fdirc,[1,m+2,2:m+1]),sizec),sizec);
      else
         fdir = ppmak(breaks,fdirc,sizec);
      end
   else
      fdir = fncmb(fndirp(f),direction);
   end
case 'r'
   error('This command does not work for rational splines. Use FNTLR instead.')
otherwise
   error('F does not appear to describe a function.')
end

function fprime = fndirp(f)
%FNDIRP Differentiate a univariate function in ppform, but keep its order.
[breaks,coefs,l,k,d] = ppbrk(f);
tmp = [k-1:-1:1]; 
fprime = ppmak(breaks,[zeros(d*l,1),coefs(:,1:k-1).*(tmp(ones(d*l,1),:))],d);
