function fprime = fnder(f,dorder)
%FNDER Differentiate a function.
%
%   FNDER(F)
%
%   returns the (representation of the) first derivative of the
%   univariate function contained in F (and in the same form).  
%
%   FNDER(F,DORDER) returns the DORDER-th derivative, with DORDER expected
%   to be of the form [d1,...,dm] in case the function in F is m-variate,
%   and, for each i=1,..,m,  di  an integer to indicate that the function
%   in F is to be differentiated di-fold with respect to its i-th argument.
%   Here, di may be negative, resulting in di-fold integration with respect
%   to the i-th argument.
%
%   This command does not work for rational splines; for them, use FNTLR 
%   instead.
%
%   For example,
%
%   fnval( fnder( sp, 2), 3.14 );
% 
%   gives the value at 3.14 of the function in sp, while
%
%   sp0 = fnint( fnder( sp ) );
%
%   gives a function that differs from sp by some constant only (namely, by
%   its value at 0).
%
%   See also FNDIR, FNTLR, FNINT.

%   C. de Boor latest change: Feb.25, 1989
%   cb 20 oct 94 : remove extraneous argument from call to  spmak
%   cb  4 nov 94 : include differentiation of bivariate tensor product splines
%                  in both B-form and ppform.
%   cb : 9 may '95 (use .' instead of '), 29aug96 (include BB-form)
%   cb : 5 oct 97 (replace use of ANS, to help compilation)
%   cb : 26oct97 (treat also multivariate functions correctly;
%                 permit DORDER to be negative to do integration)
%   cb : 05may98 (standardize help; accept earlier versions of forms)
%   cb : 11dec98 (replicate via repmat)
%   cb : 23jan00 (refer to FNTLR for rational splines)
%   cb : 09apr00 (deal with fact that matlab drops trailing singleton dims)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.11 $

if nargin<2, dorder=1; end

if ~isstruct(f), f = fn2fm(f); 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);
      if length(dorder)~=m
         error(['DORDER should be a ' num2str(m) '-vector.']), end
      sizec = [d,l.*k]; %size(coefs);
      for i=m:-1:1
         dd = prod(sizec(1:m));
         dpp = fnderp(ppmak(breaks{i},reshape(coefs,dd*l(i),k(i)),dd), ...
                      dorder(i));
         breaks{i} = dpp.breaks; sizec(m+1) = dpp.pieces*dpp.order;
         coefs = reshape(dpp.coefs,sizec);
         if m>1
             coefs = permute(coefs,[1,m+1,2:m]);
             sizec(2:m+1) = sizec([m+1,2:m]);
         end
      end
      fprime = ppmak(breaks,coefs,sizec);
   else
      fprime = fnderp(f,dorder);
   end

case 'B' % the function is in B-form or BB-form;
                           % omit trivial B-spline terms.
   [knots,coefs,n,k,d]=spbrk(f);
   if iscell(knots)       % the function is multivariate
      m = length(knots);
      if length(dorder)~=m
         error(['DORDER should be a ' num2str(m) '-vector.']), end
      sizec = [d,n];% size(coefs);
      for i=m:-1:1
         dsp = fnderb(spmak(knots{i},...
            reshape(coefs,prod(sizec(1:m)),sizec(m+1))),dorder(i));
         knots{i} = dsp.knots; sizec(m+1) = dsp.number; 
         coefs = reshape(dsp.coefs,sizec); 
         if m>1
            coefs = permute(coefs,[1,m+1,2:m]);
            sizec(2:m+1) = sizec([m+1,2:m]);
         end
      end
      fprime = spmak(knots,coefs,sizec);
   else
      fprime = fnderb(f,dorder);
   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 = fnderp(f,dorder)
%FNDERP Differentiate a univariate function in ppform.
[breaks,coefs,l,k,d]=ppbrk(f);
if k<=dorder
   fprime=ppmak([breaks(1) breaks(l+1)],zeros(d,1));
elseif dorder<0    % we are to integrate
   fprime = f;
   for j=1:(-dorder)
      fprime = fnint(fprime);
   end
else
   knew=k-dorder;
   for j=k-1:-1:knew
      coefs=coefs.*repmat([j:-1:j-k+1],d*l,1);
   end
   fprime=ppmak(breaks,coefs(:,1:knew),d);
end

function fprime = fnderb(f,dorder)
%FNDERB Differentiate a univariate function in B-form.

[t,a,n,k,d]=spbrk(f);
if k<=dorder
   fprime=spmak(t,zeros(d,n));
elseif dorder<0    % we are to integrate
   fprime = f;
   for j=1:(-dorder)
      fprime = fnint(fprime);
   end
else
   knew=k-dorder;
   for j=k-1:-1:knew
      tt=t(j+1+[0:n])-t(1:n+1); z=find(tt>0); nn=length(z);
      temp=(diff([zeros(1,d);a.'; zeros(1,d)])).';
      a=temp(:,z)./repmat(tt(z)/j,d,1);
      t=[t(z) t(n+2:n+j+1)]; n=nn;
   end
   fprime=spmak(t,a);
end
