function [zStar,fxStar,fyStar] = HConLineSeg(P,F,tStar)
% P is a length-2 cell array of structures that identify 2 distinct points in
% the plane: (P{1}.x,P{1}.y) and (P{2}.x,P{2}.y).
% F is a length-2 cell array of structures that encode the value of a
% function f(x,y) at these points together with the value of the x and y
% partial derivatives, e.g., F{k}.fvalue, F{k}.fx, and F{k}.fy , k=1:2
% tStar is a scalar that satisfies 0<=tStar<=1.
% zStar = c(tStar) where c(t) is the piecewise cubic hermite interpolant of
% the function f(P{1}.x+t*(P{2}.x-P{1}.x),P{1}.y+t*(P{2}.y-P{1}.y)) at
% t=0 and t=1.
% fxStar and fyStar are estimates of f's partial derivatives in the x and y directions
% at tStar

u0 = P{1}.x; v0 = P{1}.y; f0 = F{1}.fvalue; fx0 = F{1}.fx; fy0 = F{1}.fy;
u1 = P{2}.x; v1 = P{2}.y; f1 = F{2}.fvalue; fx1 = F{2}.fx; fy1 = F{2}.fy;

% Let ftilde(t) = f(u0+t*(u1-u0),v0+t*(v1-v0)).
% To invoke Hermite cubic interpolation, we need the value of ftilde and
% its derivatives at t=0 and t=1.

% Note that ftilde(0) = f0 and ftilde(1) = f1.

% By chain rule, ftilde'(0) is given by
fp0 = fx0*(u1-u0) + fy0*(v1-v0);

% By chain rule, ftilde'(1) is given by
fp1 = fx1*(u1-u0) + fy1*(v1-v0);

% Compute the Hermite cubic interpolant
% c(t) = a1 + a2*t +a3*t^2 + a4*t^2*(t-1) ...
[a1,a2,a3,a4] = HCubic(0,f0,fp0,1,f1,fp1);

% Evaluate at tStar...
zStar = ((a4*(tStar-1) + a3)*tStar + a2)*tStar + a1;

% Since c'(t) = (3t^2-2t)a4+2t*a3 + a2 = 3*a4*t^2 + 2*(a3-a4)*t + a2

cpStar = (3*a4*tStar + 2*(a3-a4))*tStar + a2;

fxStar = fx0+tStar*(fx1-fx0);
fyStar = fy0+tStar*(fy1-fy0);

% Getting the right c(t) = 3 pts. If you do, error is O(h^4)
% A method for fxStar and fyStar that is at least as good as linear
% interpolation, which is O(h^2) is worth 2 points. You start losing points
% if your approximation fails to be exact at the endpoints.