function p = Horner2(x, f, c)
%***************************************************************
% p = Horner2(x, f, c)
%
% This function evaluates a polynomial at c; this polynomial has
% coefficients from the Hermite interpolant vector f that were
% sampled at value(s) of x (most x-values are used twice).  The
% evaluation is recursive and uses Horner's Rule.  This function
% is intented to complement HermInt.m, but only one interpolant
% vector (f) can be used at a time.  Call this function in a FOR
% loop if multiple interpolant vectors from the output of
% HermInt.m will be used.
%
%***************************************************************

[m1, n1] = size(c); [m2, n2] = size(x);
[m3, n3] = size(f);
p = zeros(m1, n1);
if ((m2 > 1) & (n2 > 1)) | ((m3 > 1) & (n3 > 1))
    p(:) = NaN;         % If x or f is a matrix, exit
elseif (~(m3 - 1) & (n3 > (2*n2+1))) | ...
       (~(n3 - 1) & (m3 > (2*m2+1)))
    p(:) = NaN;         % If size of x (counting points twice)
                        % is less than (length(f) - 1), exit
else
    if (m2 > 1)         % Cause x to repeat number pairs
        x = reshape([x'; x'], 2*m2, 1); % x is a column vector
    else
        x = reshape([x; x], 1, 2*n2);   % x is a row vector
    end
    n = length(f);      % Counter spans the length of f
    p(:) = f(n);        % Initial value in execution
    for k = (n-1):-1:1  % Recursive loop using Horner's Rule
        p = p.*(c-x(k)) + f(k);
    end
end