function y=nakspval(a,b,c,d,knots,xx);
%% function y=nakspval(a,b,c,d,knots,xx);
%%
%%  Compute the values at the points xx of the NOT-A-KNOT cubic spline  
%%  intepolating at the ``knots'' from the coefficients of its polynomial pieces
%%
%%  INPUTS:    (from output of naknotsp.m)
%%     knots  -  x_j, including the ``inactive'' ones
%%     a - column vector of constant terms for the spline on [x_j,x_j-1] 
%%     b - column vector of coefficients of x for the spline on [x_j,x_j-1] 
%%     c - column vector of coefficients of x^2 for the spline on [x_j,x_j-1] 
%%     d - column vector of coefficients of x^3 for the spline on [x_j,x_j-1] 
%%     xx - points at which spline is to be evaluated
%%  OUTPUT:
%%     y - value of the spline
n=length(knots)-1;
I=find(xx<knots(3));
if isempty(I)==0;
   y(I)=a(1)+(xx(I)-knots(2)).*(b(1)+...
         (xx(I)-knots(2)).*(c(1)+d(1)*(xx(I)-knots(2))));
end;clear I;
for j=3:n-2;
   I = find((xx>=knots(j))&(xx<knots(j+1)));
   if isempty(I)==0;
      y(I)=a(j-1)+(xx(I)-knots(j)).*(b(j-1)+...
          (xx(I)-knots(j)).*(c(j-1)+d(j-1)*(xx(I)-knots(j))));
   end;clear I;
end;
I=find( xx>=knots(n-1));
if isempty(I)==0;
  y(I)=a(n-2)+(xx(I)-knots(n-1)).*(b(n-2)+...
         (xx(I)-knots(n-1)).*(c(n-2)+d(n-2)*(xx(I)-knots(n-1))));
end;clear I;
