[t,jhigh,indexmlv,x,ileft,vnikx]=bsplvn(t,jhigh,indexmlv,x,ileft,vnikx);

function [t,jhigh,indexmlv,x,ileft,vnikx]=bsplvn(t,jhigh,indexmlv,x,ileft,vnikx);
persistent deltam deltap firstCall i imjp1 ipj j jp1 jp1ml l vm vmprev ; if isempty(firstCall),firstCall=1;end; 

if isempty(deltam), deltam=zeros(1,20); end;
if isempty(deltap), deltap=zeros(1,20); end;
if isempty(vm), vm=0; end;
if isempty(vmprev), vmprev=0; end;
if isempty(i), i=0; end;
if isempty(imjp1), imjp1=0; end;
if isempty(ipj), ipj=0; end;
if isempty(j), j=0; end;
if isempty(jp1), jp1=0; end;
if isempty(jp1ml), jp1ml=0; end;
if isempty(l), l=0; end;
%***BEGIN PROLOGUE  BSPLVN
%***SUBSIDIARY
%***PURPOSE  Subsidiary to FC
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (BSPLVN-S, DFSPVN-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
% Calculates the value of all possibly nonzero B-splines at *X* of
%  order MAX(JHIGH,(J+1,INDEX-1)) on *T*.
%
%***SEE ALSO  FC
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   780801  DATE WRITTEN
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%***end PROLOGUE  BSPLVN
t_shape=size(t);t=reshape(t,1,[]);
vnikx_shape=size(vnikx);vnikx=reshape(vnikx,1,[]);
if firstCall,   j=[1];  end;
if firstCall, [deltam([1:20]) ]=[0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]; end;
if firstCall, [ deltap([1:20])]=[0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  BSPLVN
if( indexmlv~=2 )
j = 1;
vnikx(1) = 1.;
if( j>=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;
return;
end;
end;
%
while (1);
ipj = fix(ileft + j);
deltap(j) = t(ipj) - x;
imjp1 = fix(ileft - j + 1);
deltam(j) = x - t(imjp1);
vmprev = 0.;
jp1 = fix(j + 1);
for l = 1 : j;
jp1ml = fix(jp1 - l);
vm = vnikx(l)./(deltap(l)+deltam(jp1ml));
vnikx(l) = vm.*deltap(l) + vmprev;
vmprev = vm.*deltam(jp1ml);
end; l = fix(j+1);
vnikx(jp1) = vmprev;
j = fix(jp1);
if( j>=jhigh )
break;
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;
end %subroutine bsplvn
%DECK BSPPP