function [y,yp] = hermiteval(rho,t,tau,sin,w,D)
%
% [y,yp] = hermiteval( rho, t, tau, s, w, D )
% (c)  Robert M. Corless, 2007
% Evaluate a Hermite Interpolant and its derivative
% at a potentially large number of points, t.
% y^(j)(tau(i))/j! = rho(k) where k = s(1)+s(2)+...s(i-1)+j+1
% confluencies s(i)
% generalized barycentric weights w
% differentiation matrix D (Taylor form if rho is)
%
% based on a code fragment in Berrut and Trefethen
% exact node value requests replaced with input data
% or (possibly) computed derivatives if t=tau(i) and s(i)=1
% 

[n,dummy] = size(tau(:));
[m,dummy] = size(t(:));   % expect m >> n

if max(size(sin))==1,
    s = sin*ones(n,1);
else
    s = sin;
end;
smax = max(s);

numer = zeros(m,1);
denom = zeros(m,1);
exact = zeros(m,1); % which t-values are exactly nodes

temp = zeros(m,smax);

if nargout>1,
    dert   = D*rho(:);  % D must match rho's form; Taylor or not
    numerd = zeros(m,1);
end;

ik = 0;
for i=1:n,
    tdiff = t(:) - tau(i);
    % B-T trick had a typo, fixed now!  i, not 1
    exact(tdiff==0) = i;  % nonzero is true
    % Schneider-Werner barycentric Hermite form
    for j=1:s(i),
        temp(:,j) = w(i,j)./(tdiff.^j);
        rhoik = zeros(m,1);
        if nargout>1,
            rhoikd = zeros(m,1);
        end;
        for k=j:-1:1,
            rhoik = rho(ik+k) + tdiff.*rhoik;
            if nargout>1,
                rhoikd = dert(ik+k) + tdiff.*rhoikd;
            end;
        end;
        numer  = numer + temp(:,j).*rhoik;
        if nargout>1,
            numerd = numerd + temp(:,j).*rhoikd;
        end;
        denom = denom + temp(:,j);
    end;
    ik = ik + s(i);
end;
y = numer./denom;
if nargout>1,
    yp= numerd./denom;
end;


% for replacing node-poles, flatten data for easy access.

ik = 1+cumsum(s(:)');
ik = [1,ik(1:n-1)]; % [1,s(1)+1,s(1)+s(2)+1,...,] points at y-vals
rhoflat = rho(ik); 

jj=find(exact);  % Berrut-Trefethen trick

y(jj) = rhoflat(exact(jj));
if nargout>1,
    dertflat = dert(ik);
    yp(jj) = dertflat(exact(jj));
end;