function [t,p,c] = hermite(x,f,d,m) 
% HERMITE: evaluation of a Hermite spline
% p = hermite_spline(x,f,d,m) 
% m+1 values p(j) per interval at equally spaced points t(j) and Bernstein coefficients c(1:4,:)
% of a Hermite spline with data f(k), d(k) at x(k)

% convert to column vectors (if necessary)
x = x(:); f = f(:); d = d(:);

% evaluation points
h = diff(x); 
t = ones(m+1,1)*x(1:end-1)' + ([0:m]'/m)*h'; 
t = t(:);

% Bernstein coefficients
c = [f(1:end-1), f(1:end-1)+d(1:end-1).*h/3, ...
   f(2:end)-d(2:end).*h/3, f(2:end)]';

% evaluation 
b = bernstein(3,[0:m]'/m);
p = b*c; 
p = p(:);
 
