function z = HConTriangleCVL1(P,F,alpha)
% P is a length-3 cell array of structures that identify 3 distinct points in
% the plane: (P{1}.x,P{1}.y), (P{2}.x,P{2}.y), and (P{3}.x,P{3}.y).
% F is a length-3 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:3
% alpha is a length-3 vector of nonnegatve numbers that sum to one.
% z is an estimate of f(xtilde,ytilde) where
% xtilde = alpha(1)*P{1}.x + alpha(2)*P{2}.x + alpha(3)*P{3}.x
% ytilde = alpha(1)*P{1}.y + alpha(2)*P{2}.y + alpha(3)*P{3}.y
%
% In other words, using the function and derivative information known at the
% vertices of triangle with vertices at P{1}, P{2}, and P{3}, z is an estimate
% of the value of the function at a point inside the triangle.

% Let Q be the point of interpolation.
% Let R be the intersection of the ray from P{1} thru Q and the line segment that
% connects P{2} and P{3}. Our plan is to use HConLineSeg to estimate f at R
%using data at P{2} and P{3} and then
% use HConLineSeg again to estimate f at Q using data at P{1} and R. It makes sense
% for P{1} to be the vertex >>furthest<< away from Q for then Q is closet to the
% side P{2}P{3}.

if alpha(2) == min(alpha)
temp=P{2}; P{2}=P{1}; P{1}=temp;
temp=F{2}; F{2}=F{1}; F{1}=temp;
temp=alpha(2); alpha(2)=alpha(1); alpha(1)=temp;
elseif alpha(3) == min(alpha)
temp=P{3}; P{3}=P{1}; P{1}=temp;
temp=F{3}; F{3}=F{1}; F{1}=temp;
temp=alpha(3); alpha(3)=alpha(1); alpha(1)=temp;
end

% Encode the points as 2-vectors...
P1 = [P{1}.x;P{1}.y];
P2 = [P{2}.x;P{2}.y];
P3 = [P{3}.x;P{3}.y];
% The point of interpolation...
Q = alpha(1)*P1 + alpha(2)*P2 + alpha(3)*P3;
% Find scalars t(1) and t(2) so that
% P1 + t(1)*(Q-P1) = P2 + t(2)*(P3-P2)
% This is where the "ray" from P1 thru Q intersects the line segment that % connects P2 and P3.
t = [ Q-P1 P2-P3]\(P2-P1);
%fprintf('\nt(1) = %10.3f t(2) = %10.3f\n',t(1),t(2))
% The point of intersection R = P2 + t(2)*(P3-P2) = P1 + t(1)*(Q-P1).
R = struct('x',P2(1) + t(2)*(P3(1)-P2(1)),'y',P2(2) + t(2)*(P3(2)-P2(2)));

% Estimate MyF and its partials at R using data at P2 and P3...
[fR,fxR,fyR] = HConLineSeg({P{2},P{3}},{F{2},F{3}},t(2));
% Estimate MyF at Q using data at P1 and R...
[fQ,fxQ,fyQ] = HConLineSeg({P{1},R},{F{1},struct('fvalue',fR,'fx',fxR,'fy',fyR)},t(1));
z = fQ;

% Because the 2nd call to HConLineSeg lives off of approximate endpoint
% derivative (that ar O(h^2) accurate), the overall method is about O(h^2).
% About -2 points if your method doesn't produce O(h^4) accuracy when Q is
% on a triangle edge.
% If you do the intersection computations wrong and/or overlook the case
% when Q is on the triangle boundary, then your computed z might be an NaN
% or an inf. Up to -2points if that happens