[ppvalresult,ldc,c,xi,lxi,k,ideriv,x,inppv]=ppval(ldc,c,xi,lxi,k,ideriv,x,inppv);

function [ppvalresult,ldc,c,xi,lxi,k,ideriv,x,inppv]=ppval(ldc,c,xi,lxi,k,ideriv,x,inppv);
ppvalresult=[];
persistent dx fltk i j ndummy ppval ; 

if isempty(ppvalresult), ppvalresult=0; end;
%***BEGIN PROLOGUE  PPVAL
%***PURPOSE  Calculate the value of the IDERIV-th derivative of the
%            B-spline from the PP-representation.
%***LIBRARY   SLATEC
%***CATEGORY  E3, K6
%***TYPE      SINGLE PRECISION (PPVAL-S, DPPVAL-D)
%***KEYWORDS  B-SPLINE, DATA FITTING, INTERPOLATION, SPLINES
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Written by Carl de Boor and modified by D. E. Amos
%
%     Abstract
%         PPVAL is the PPVALU function of the reference.
%
%         PPVAL calculates (at X) the value of the IDERIV-th
%         derivative of the B-spline from the PP-representation
%         (C,XI,LXI,K).  The Taylor expansion about XI(J) for X in
%         the interval XI(J) .LE. X .LT. XI(J+1) is evaluated, J=1,LXI.
%         Right limiting values at X=XI(J) are obtained.  PPVAL will
%         extrapolate beyond XI(1) and XI(LXI+1).
%
%         To obtain left limiting values (left derivatives) at XI(J),
%         replace LXI by J-1 and set X=XI(J),J=2,LXI+1.
%
%     Description of Arguments
%         Input
%          LDC     - leading dimension of C matrix, LDC .GE. K
%          C       - matrix of dimension at least (K,LXI) containing
%                    right derivatives at break points XI(*).
%          XI      - break point vector of length LXI+1
%          LXI     - number of polynomial pieces
%          K       - order of B-spline, K .GE. 1
%          IDERIV  - order of the derivative, 0 .LE. IDERIV .LE. K-1
%                    IDERIV=0 gives the B-spline value
%          X       - argument, XI(1) .LE. X .LE. XI(LXI+1)
%          INPPV   - an initialization parameter which must be set
%                    to 1 the first time PPVAL is called.
%
%         Output
%          INPPV   - INPPV contains information for efficient process-
%                    ing after the initial call and INPPV must not
%                    be changed by the user.  Distinct splines require
%                    distinct INPPV parameters.
%          PPVAL   - value of the IDERIV-th derivative at X
%
%     Error Conditions
%         Improper input is a fatal error
%
%***REFERENCES  Carl de Boor, Package for calculating with B-splines,
%                 SIAM Journal on Numerical Analysis 14, 3 (June 1977),
%                 pp. 441-472.
%***ROUTINES CALLED  INTRV, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800901  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  PPVAL
%
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(ndummy), ndummy=0; end;
if isempty(dx), dx=0; end;
if isempty(fltk), fltk=0; end;
xi_shape=size(xi);xi=reshape(xi,1,[]);
c_shape=size(c);c=reshape([c(:).',zeros(1,ceil(numel(c)./prod([ldc])).*prod([ldc])-numel(c))],abs(ldc),[]);
%***FIRST EXECUTABLE STATEMENT  PPVAL
ppvalresult = 0.0e0;
if( k<1 )
xermsg('SLATEC','PPVAL','K DOES NOT SATISFY K.GE.1',2,1);
elseif( ldc<k ) ;
%
%
xermsg('SLATEC','PPVAL','LDC DOES NOT SATISFY LDC.GE.K',2,1);
elseif( lxi<1 ) ;
xermsg('SLATEC','PPVAL','LXI DOES NOT SATISFY LXI.GE.1',2,1);
elseif( ideriv<0 || ideriv>=k ) ;
xermsg('SLATEC','PPVAL','IDERIV DOES NOT SATISFY 0.LE.IDERIV.LT.K',2,1);
else;
i = fix(k - ideriv);
fltk = i;
[xi,lxi,x,inppv,i,ndummy]=intrv(xi,lxi,x,inppv,i,ndummy);
dx = x - xi(i);
j = fix(k);
while( true );
ppvalresult =(ppvalresult./fltk).*dx + c(j,i);
j = fix(j - 1);
fltk = fltk - 1.0e0;
if( fltk<=0.0e0 )
break;
end;
end;
end;
xi_shape=zeros(xi_shape);xi_shape(:)=xi(1:numel(xi_shape));xi=xi_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',xi); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(7)), assignin('caller','FUntemp',x); evalin('caller',[inputname(7),'=FUntemp;']); end
if csnil&&~isempty(inputname(4)), assignin('caller','FUntemp',lxi); evalin('caller',[inputname(4),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',ldc); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(5)), assignin('caller','FUntemp',k); evalin('caller',[inputname(5),'=FUntemp;']); end
if csnil&&~isempty(inputname(8)), assignin('caller','FUntemp',inppv); evalin('caller',[inputname(8),'=FUntemp;']); end
if csnil&&~isempty(inputname(6)), assignin('caller','FUntemp',ideriv); evalin('caller',[inputname(6),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',c); evalin('caller',[inputname(2),'=FUntemp;']); end
end %function ppval
%DECK PROC