double precision function seval(n, u, x, y, b, c, d) integer n double precision u, x(n), y(n), b(n), c(n), d(n) c c this subroutine evaluates the cubic spline function c c seval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2 + d(i)*(u-x(i))**3 c c where x(i) .lt. u .lt. x(i+1), using horner's rule c c if u .lt. x(1) then i = 1 is used. c if u .ge. x(n) then i = n is used. c c input.. c c n = the number of data points c u = the abscissa at which the spline is to be evaluated c x,y = the arrays of data abscissas and ordinates c b,c,d = arrays of spline coefficients computed by spline c c if u is not in the same interval as the previous call, then a c binary search is performed to determine the proper interval. c integer i, j, k double precision dx data i/1/ if ( i .ge. n ) i = 1 if ( u .lt. x(i) ) go to 10 if ( u .le. x(i+1) ) go to 30 c c binary search c 10 i = 1 j = n+1 20 k = (i+j)/2 if ( u .lt. x(k) ) j = k if ( u .ge. x(k) ) i = k if ( j .gt. i+1 ) go to 20 c c evaluate spline c 30 dx = u - x(i) seval = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i))) return end