subroutine DBSPVT ( t, k, index, x, left, biatx, dbiatx ) c Derived from deBoor's DBSPVB by Eric Grosse, 4 May 1987, 23 Nov 1989 c Calling sequence is unchanged except for addition of dbiatx, c which is defined by dbiatx(i,m) = d B(left-k+i) / d t(left-k+m) c for 1<=i<=k and 1<=m<=2*k. c For my own sanity, jhigh was renamed k. It is still the spline order. parameter(jmax = 20) parameter(mmax = 40) integer index,k,left,i,j,jp1,m,mm DOUBLE PRECISION dbiatx(k,1), dsaved(mmax), dterm(mmax), denom DOUBLE PRECISION biatx(k),t(1),x,deltal(jmax),deltar(jmax),saved,t +erm data j/1/ save j,deltal,deltar go to (10,20), index 10 j = 1 biatx(1) = 1.D0 mm = 2*k do 1001 m=1,mm 1001 dbiatx(1,m) = 0.D0 if (j .ge. k) go to 99 c 20 jp1 = j + 1 deltar(j) = t(left+j) - x deltal(j) = x - t(left+1-j) saved = 0.D0 do 1002 m=1,mm 1002 dsaved(m) = 0.D0 do 26 i=1,j denom = 1/(deltar(i) + deltal(jp1-i)) term = biatx(i)*denom do 1003 m=1,mm dterm(m) = dbiatx(i,m)*denom if(i .eq.m-k) dterm(m) = dterm(m) - term*denom if(i-j.eq.m-k) dterm(m) = dterm(m) + term*denom 1003 continue biatx(i) = saved + deltar(i)*term do 1004 m=1,mm dbiatx(i,m) = dsaved(m) + deltar(i)*dterm(m) if(i.eq.m-k) dbiatx(i,m) = dbiatx(i,m) + term 1004 continue saved = deltal(jp1-i)*term do 1005 m=1,mm dsaved(m) = deltal(jp1-i)*dterm(m) if(i-j.eq.m-k) dsaved(m) = dsaved(m) - term 1005 continue 26 continue biatx(jp1) = saved do 1006 m=1,mm 1006 dbiatx(jp1,m) = dsaved(m) j = jp1 if (j .lt. k) go to 20 c 99 return end