function v=derppval(pp,xx)
%PPVAL  Evaluate piecewise polynomial.
%   V = PPVAL(PP,XX) returns the value at the points XX of the 
%   piecewise polynomial contained in PP, as constructed by SPLINE 
%   (or MKPP).
%
%   Revised by B.A. Broughton on 10/01/1999 to calculate the derivative
%   for a pre-fitted cubic spline in the points xx
%
%   See also SPLINE.

%   Carl de Boor 7-2-86
%   Revised 10-14-97 CB to speed up locating points in mesh (as in
%   PPUAL) and to handle vector-valued functions.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.8 $  $Date: 1997/11/21 23:41:00 $

[mx,nx] = size(xx); lx = mx*nx; xs = reshape(xx,1,lx);
% if necessary, sort xx 
tosort=0;
if any(diff(xs)<0)
   tosort=1;[xs,ix]=sort(xs);
end

% take apart pp
[x,c,l,k,d]=unmkpp(pp);

% for each data point, compute its breakpoint interval
[ignored,index] = sort([x(1:l) xs]);
index = max([find(index>l)-(1:lx);ones(1,lx)]);

% now go to local coordinates ...
xs = xs-x(index);

if d>1 % ... replicate xs and index in case pp is vector-valued ...
   xs = reshape(xs(ones(d,1),:),1,d*lx);
   index = d*index; temp = [-d:-1].';
   index = reshape(1+index(ones(d,1),:)+temp(:,ones(1,lx)), d*lx, 1 );
end

% ... and apply nested multiplication to get derivative:
   v = 3*c(index,1).';
   for i=2:k-1
      v = xs.*v + (4-i)*c(index,i).';
   end
v = reshape(v,d,lx);
if tosort>0, v(:,ix) = v; end
v = reshape(v,d*mx,nx);