function p = hermite_spline(x,f,d,t)

% function p = hermite_spline(x,f,d,t)
% values p(t) of hermite spline 
% with values f(k) and derivatives d(k) at x(k)

n = length(x);
delta = x(2:n)-x(1:n-1);
cm = f(2:n) - (delta/3).*d(2:n);
cp = f(1:n-1) + (delta/3).*d(1:n-1);

for k=1:n-1;
   c = [f(k); cp(k); cm(k); f(k+1)];
   ind = find((x(k)<=t) & (t<=x(k+1)));
   p(ind) = bezier_curve(c, (t(ind)-x(k))/delta(k));
end;


function p = bezier_curve(c,t)
p = c(1)*(1-t).^3+3*c(2).*t.*(1-t).^2+3*c(3).*t.^2.*(1-t) +c(4).*t.^3;



