FUNCTION USPF(X,A,N,XA,XB)
*     Returned is the function value at X of B-spline function.
*     A(N) is array of coefficients, obtained from USPCN.
*     XA and  XB are X-values of first and last function value 
*     Extrapolation is done for X outside the range XA ... XB.
      REAL A(N)
      R=FLOAT(N-1)*(X-XA)/(XB-XA)
      I=MAX(1,MIN(INT(R),N-3))
      Z=R-FLOAT(I)
      USPF
     +    =(A(I)+4.0*A(I+1)+A(I+2)+Z*(3.0*(-A(I)+A(I+2))
     +     +Z*(3.0*(A(I)-2.0*A(I+1)+A(I+2))
     +     +Z*(-A(I)+3.0*A(I+1)-3.0*A(I+2)+A(I+3)))))/6.0
      END