[t,ad,n,k,nderiv,x,inev,svalue,work]=bspev(t,ad,n,k,nderiv,x,inev,svalue,work);

function [t,ad,n,k,nderiv,x,inev,svalue,work]=bspev(t,ad,n,k,nderiv,x,inev,svalue,work);
%***BEGIN PROLOGUE  BSPEV
%***PURPOSE  Calculate the value of the spline and its derivatives from
%            the B-representation.
%***LIBRARY   SLATEC
%***CATEGORY  E3, K6
%***TYPE      SINGLE PRECISION (BSPEV-S, DBSPEV-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
%         BSPEV is the BSPLEV routine of the reference.
%
%         BSPEV calculates the value of the spline and its derivatives
%         at X from the B-representation (T,A,N,K) and returns them
%         in SVALUE(I),I=1,NDERIV, T(K) .LE. X .LE. T(N+1).  AD(I) can
%         be the B-spline coefficients A(I), I=1,N if NDERIV=1.  Other-
%         wise AD must be computed before hand by a call to BSPDR (T,A,
%         N,K,NDERIV,AD).  If X=T(I),I=K,N, right limiting values are
%         obtained.
%
%         To compute left derivatives or left limiting values at a
%         knot T(I), replace N by I-1 and set X=T(I), I=K+1,N+1.
%
%         BSPEV calls INTRV, BSPVN
%
%     Description of Arguments
%         Input
%          T       - knot vector of length N+K
%          AD      - vector of length (2*N-NDERIV+1)*NDERIV/2 containing
%                    the difference table from BSPDR.
%          N       - number of B-spline coefficients
%                    N = sum of knot multiplicities-K
%          K       - order of the B-spline, K .GE. 1
%          NDERIV  - number of derivatives, 1 .LE. NDERIV .LE. K.
%                    NDERIV=1 gives the zero-th derivative = function
%                    value
%          X       - argument, T(K) .LE. X .LE. T(N+1)
%          INEV    - an initialization parameter which must be set
%                    to 1 the first time BSPEV is called.
%
%         Output
%          INEV    - INEV contains information for efficient process-
%                    ing after the initial call and INEV must not
%                    be changed by the user.  Distinct splines require
%                    distinct INEV parameters.
%          SVALUE  - vector of length NDERIV containing the spline
%                    value in SVALUE(1) and the NDERIV-1 derivatives
%                    in the remaining components.
%          WORK    - work vector of length 3*K
%
%     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  BSPVN, INTRV, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800901  DATE WRITTEN
%   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)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  BSPEV
%
persistent i id iwork jj kp1 kp1mn l left ll mflag summlv ; 

if isempty(i), i=0; end;
if isempty(id), id=0; end;
if isempty(iwork), iwork=0; end;
if isempty(jj), jj=0; end;
if isempty(kp1), kp1=0; end;
if isempty(kp1mn), kp1mn=0; end;
if isempty(l), l=0; end;
if isempty(left), left=0; end;
if isempty(ll), ll=0; end;
if isempty(mflag), mflag=0; end;
if isempty(summlv), summlv=0; end;
%     DIMENSION T(N+K)
t_shape=size(t);t=reshape(t,1,[]);
ad_shape=size(ad);ad=reshape(ad,1,[]);
svalue_shape=size(svalue);svalue=reshape(svalue,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
%***FIRST EXECUTABLE STATEMENT  BSPEV
if( k<1 )
%
%
xermsg('SLATEC','BSPEV','K DOES NOT SATISFY K.GE.1',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
ad_shape=zeros(ad_shape);ad_shape(:)=ad(1:numel(ad_shape));ad=ad_shape;
svalue_shape=zeros(svalue_shape);svalue_shape(:)=svalue(1:numel(svalue_shape));svalue=svalue_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
elseif( n<k ) ;
xermsg('SLATEC','BSPEV','N DOES NOT SATISFY N.GE.K',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
ad_shape=zeros(ad_shape);ad_shape(:)=ad(1:numel(ad_shape));ad=ad_shape;
svalue_shape=zeros(svalue_shape);svalue_shape(:)=svalue(1:numel(svalue_shape));svalue=svalue_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
elseif( nderiv<1 || nderiv>k ) ;
xermsg('SLATEC','BSPEV','NDERIV DOES NOT SATISFY 1.LE.NDERIV.LE.K',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
ad_shape=zeros(ad_shape);ad_shape(:)=ad(1:numel(ad_shape));ad=ad_shape;
svalue_shape=zeros(svalue_shape);svalue_shape(:)=svalue(1:numel(svalue_shape));svalue=svalue_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
else;
id = fix(nderiv);
[t,dumvar2,x,inev,i,mflag]=intrv(t,n+1,x,inev,i,mflag);
if( x>=t(k) )
if( mflag~=0 )
if( x>t(i) )
xermsg('SLATEC','BSPEV','X IS NOT IN T(K).LE.X.LE.T(N+1)',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
ad_shape=zeros(ad_shape);ad_shape(:)=ad(1:numel(ad_shape));ad=ad_shape;
svalue_shape=zeros(svalue_shape);svalue_shape(:)=svalue(1:numel(svalue_shape));svalue=svalue_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
while (1);
if( i==k )
xermsg('SLATEC','BSPEV','A LEFT LIMITING VALUE CANNOT BE OBTAINED AT T(K)',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
ad_shape=zeros(ad_shape);ad_shape(:)=ad(1:numel(ad_shape));ad=ad_shape;
svalue_shape=zeros(svalue_shape);svalue_shape(:)=svalue(1:numel(svalue_shape));svalue=svalue_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
else;
i = fix(i - 1);
if( x~=t(i) )
break;
end;
end;
end;
end;
%
% *I* HAS BEEN FOUND IN (K,N) SO THAT T(I) .LE. X .LT. T(I+1)
%     (OR .LE. T(I+1), IF T(I) .LT. T(I+1) = T(N+1) ).
kp1mn = fix(k + 1 - id);
kp1 = fix(k + 1);
[t,kp1mn,k,dumvar4,x,i,dumvar7,dumvar8,iwork]=bspvn(t,kp1mn,k,1,x,i,work(sub2ind(size(work),max(1,1)):end),work(sub2ind(size(work),max(kp1,1)):end),iwork);   dumvar7i=find((work(sub2ind(size(work),max(1,1)):end))~=(dumvar7));dumvar8i=find((work(sub2ind(size(work),max(kp1,1)):end))~=(dumvar8));   work(1-1+dumvar7i)=dumvar7(dumvar7i); work(kp1-1+dumvar8i)=dumvar8(dumvar8i); 
jj =fix(fix(((n+n-id+2).*(id-1))./2));
%     ADIF(LEFTPL,ID) = AD(LEFTPL-ID+1 + (2*N-ID+2)*(ID-1)/2)
%     LEFTPL = LEFT + L
while (1);
left = fix(i - kp1mn);
summlv = 0.0e0;
ll = fix(left + jj + 2 - id);
for l = 1 : kp1mn;
summlv = summlv + work(l).*ad(ll);
ll = fix(ll + 1);
end; l = fix(kp1mn+1);
svalue(id) = summlv;
id = fix(id - 1);
if( id==0 )
%
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
ad_shape=zeros(ad_shape);ad_shape(:)=ad(1:numel(ad_shape));ad=ad_shape;
svalue_shape=zeros(svalue_shape);svalue_shape(:)=svalue(1:numel(svalue_shape));svalue=svalue_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
else;
jj = fix(jj -(n-id+1));
kp1mn = fix(kp1mn + 1);
[t,kp1mn,k,dumvar4,x,i,dumvar7,dumvar8,iwork]=bspvn(t,kp1mn,k,2,x,i,work(sub2ind(size(work),max(1,1)):end),work(sub2ind(size(work),max(kp1,1)):end),iwork);   dumvar7i=find((work(sub2ind(size(work),max(1,1)):end))~=(dumvar7));dumvar8i=find((work(sub2ind(size(work),max(kp1,1)):end))~=(dumvar8));   work(1-1+dumvar7i)=dumvar7(dumvar7i); work(kp1-1+dumvar8i)=dumvar8(dumvar8i); 
continue;
end;
break;
%20
end;
end;
xermsg('SLATEC','BSPEV','X IS NOT IN T(K).LE.X.LE.T(N+1)',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
ad_shape=zeros(ad_shape);ad_shape(:)=ad(1:numel(ad_shape));ad=ad_shape;
svalue_shape=zeros(svalue_shape);svalue_shape(:)=svalue(1:numel(svalue_shape));svalue=svalue_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
ad_shape=zeros(ad_shape);ad_shape(:)=ad(1:numel(ad_shape));ad=ad_shape;
svalue_shape=zeros(svalue_shape);svalue_shape(:)=svalue(1:numel(svalue_shape));svalue=svalue_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end %subroutine bspev
%DECK BSPLVD