[t,jhigh,k,indexmlv,x,ileft,vnikx,work,iwork]=bspvn(t,jhigh,k,indexmlv,x,ileft,vnikx,work,iwork);

function [t,jhigh,k,indexmlv,x,ileft,vnikx,work,iwork]=bspvn(t,jhigh,k,indexmlv,x,ileft,vnikx,work,iwork);
%***BEGIN PROLOGUE  BSPVN
%***PURPOSE  Calculate the value of all (possibly) nonzero basis
%            functions at X.
%***LIBRARY   SLATEC
%***CATEGORY  E3, K6
%***TYPE      SINGLE PRECISION (BSPVN-S, DBSPVN-D)
%***KEYWORDS  EVALUATION OF B-SPLINE
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Written by Carl de Boor and modified by D. E. Amos
%
%     Abstract
%         BSPVN is the BSPLVN routine of the reference.
%
%         BSPVN calculates the value of all (possibly) nonzero basis
%         functions at X of order MAX(JHIGH,(J+1)*(INDEX-1)), where
%         T(K) .LE. X .LE. T(N+1) and J=IWORK is set inside the routine
%         on the first call when INDEX=1.  ILEFT is such that T(ILEFT)
%         .LE. X .LT. T(ILEFT+1).  A call to INTRV(T,N+1,X,ILO,ILEFT,
%         MFLAG) produces the proper ILEFT.  BSPVN calculates using the
%         basic algorithm needed in BSPVD.  If only basis functions are
%         desired, setting JHIGH=K and INDEX=1 can be faster than
%         calling BSPVD, but extra coding is required for derivatives
%         (INDEX=2) and BSPVD is set up for this purpose.
%
%         Left limiting values are set up as described in BSPVD.
%
%     Description of Arguments
%         Input
%          T       - knot vector of length N+K, where
%                    N = number of B-spline basis functions
%                    N = sum of knot multiplicities-K
%          JHIGH   - order of B-spline, 1 .LE. JHIGH .LE. K
%          K       - highest possible order
%          INDEX   - INDEX = 1 gives basis functions of order JHIGH
%                          = 2 denotes previous entry with WORK, IWORK
%                              values saved for subsequent calls to
%                              BSPVN.
%          X       - argument of basis functions,
%                    T(K) .LE. X .LE. T(N+1)
%          ILEFT   - largest integer such that
%                    T(ILEFT) .LE. X .LT. T(ILEFT+1)
%
%         Output
%          VNIKX   - vector of length K for spline values.
%          WORK    - a work vector of length 2*K
%          IWORK   - a work parameter.  Both WORK and IWORK contain
%                    information necessary to continue for INDEX = 2.
%                    When INDEX = 1 exclusively, these are scratch
%                    variables and can be used for other purposes.
%
%     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  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  BSPVN
%
persistent imjp1 ipj jp1 jp1ml l vm vmprev ; 

if isempty(imjp1), imjp1=0; end;
if isempty(ipj), ipj=0; end;
if isempty(jp1), jp1=0; end;
if isempty(jp1ml), jp1ml=0; end;
if isempty(l), l=0; end;
if isempty(vm), vm=0; end;
if isempty(vmprev), vmprev=0; end;
%     DIMENSION T(ILEFT+JHIGH)
t_shape=size(t);t=reshape(t,1,[]);
vnikx_shape=size(vnikx);vnikx=reshape(vnikx,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
%     CONTENT OF J, DELTAM, DELTAP IS EXPECTED UNCHANGED BETWEEN CALLS.
%     WORK(I) = DELTAP(I), WORK(K+I) = DELTAM(I), I = 1,K
%***FIRST EXECUTABLE STATEMENT  BSPVN
if( k<1 )
%
%
xermsg('SLATEC','BSPVN','K DOES NOT SATISFY K.GE.1',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
vnikx_shape=zeros(vnikx_shape);vnikx_shape(:)=vnikx(1:numel(vnikx_shape));vnikx=vnikx_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
elseif( jhigh>k || jhigh<1 ) ;
xermsg('SLATEC','BSPVN','JHIGH DOES NOT SATISFY 1.LE.JHIGH.LE.K',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
vnikx_shape=zeros(vnikx_shape);vnikx_shape(:)=vnikx(1:numel(vnikx_shape));vnikx=vnikx_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
elseif( indexmlv<1 || indexmlv>2 ) ;
xermsg('SLATEC','BSPVN','INDEX IS NOT 1 OR 2',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
vnikx_shape=zeros(vnikx_shape);vnikx_shape(:)=vnikx(1:numel(vnikx_shape));vnikx=vnikx_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
else;
if( x<t(ileft) || x>t(ileft+1) )
xermsg('SLATEC','BSPVN','X DOES NOT SATISFY T(ILEFT).LE.X.LE.T(ILEFT+1)',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
vnikx_shape=zeros(vnikx_shape);vnikx_shape(:)=vnikx(1:numel(vnikx_shape));vnikx=vnikx_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
else;
if( indexmlv~=2 )
iwork = 1;
vnikx(1) = 1.0e0;
if( iwork>=jhigh )
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
vnikx_shape=zeros(vnikx_shape);vnikx_shape(:)=vnikx(1:numel(vnikx_shape));vnikx=vnikx_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
end;
%
while (1);
ipj = fix(ileft + iwork);
work(iwork) = t(ipj) - x;
imjp1 = fix(ileft - iwork + 1);
work(k+iwork) = x - t(imjp1);
vmprev = 0.0e0;
jp1 = fix(iwork + 1);
for l = 1 : iwork;
jp1ml = fix(jp1 - l);
vm = vnikx(l)./(work(l)+work(k+jp1ml));
vnikx(l) = vm.*work(l) + vmprev;
vmprev = vm.*work(k+jp1ml);
end; l = fix(iwork+1);
vnikx(jp1) = vmprev;
iwork = fix(jp1);
if( iwork>=jhigh )
break;
end;
%20
end;
end;
%
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
vnikx_shape=zeros(vnikx_shape);vnikx_shape(:)=vnikx(1:numel(vnikx_shape));vnikx=vnikx_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;
vnikx_shape=zeros(vnikx_shape);vnikx_shape(:)=vnikx(1:numel(vnikx_shape));vnikx=vnikx_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end %subroutine bspvn
%DECK BSQAD