function y=cusplval(a,b,c,d,knots,xx);
%% function y=cusplval(a,b,c,d,knots,xx);
%%
%%  Compute the values at the points xx of the cubic spline with knots 
%%  ``knots'' from the coefficients of its polynomial pieces
%%
%%  INPUTS:    (from output of natsp3_4.m or clsp3_5.m)
%%     knots  -  x_j
%%     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);
if length(a)~=n-1;
message = 'knots and coefficients not matched' 
end;
I=find(xx<knots(2));
if isempty(I)==0;
   y(I)=a(1)+(xx(I)-knots(1)).*(b(1)+...
         (xx(I)-knots(1)).*(c(1)+d(1)*(xx(I)-knots(1))));
end;clear I;
for j=2:n-2;
   I = find((xx>=knots(j))&(xx<knots(j+1)));
   if isempty(I)==0;
      y(I)=a(j)+(xx(I)-knots(j)).*(b(j)+...
          (xx(I)-knots(j)).*(c(j)+d(j)*(xx(I)-knots(j))));
   end;clear I;
end;
I=find( xx>=knots(n-1));
if isempty(I)==0;
  y(I)=a(n-1)+(xx(I)-knots(n-1)).*(b(n-1)+...
         (xx(I)-knots(n-1)).*(c(n-1)+d(n-1)*(xx(I)-knots(n-1))));
end;clear I;
