function z = HConTriangleCVL2(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.

%Check if Q is on an edge..
if alpha(1)<=10^-6
% On the P2-P3 edge
[z,fxStar,fyStar] = HConLineSeg({P{2},P{3}},{F{2},F{3}},alpha(3));
return
elseif alpha(2)<=10^-6
% On the P3-P1 edge
[z,fxStar,fyStar] = HConLineSeg({P{3},P{1}},{F{3},F{1}},alpha(1));
return
elseif alpha(3)<=10^-6
% On the P1-P2 edge
[z,fxStar,fyStar] = HConLineSeg({P{1},P{2}},{F{1},F{2}},alpha(2));
return
end
% 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.
% 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);
% 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));

% Let (x0,y0) be the coordinates of the point Q...
x0 = Q(1);
y0 = Q(2);
% and let x and y be column-4 vectors
x = [P{1}.x; P{2}.x; P{3}.x; R.x];
y = [P{1}.y; P{2}.y; P{3}.y; R.y];
% Our goal is to interpolate the F, FX, and Fy values at the triangle vertices
% (x(1),y(1)), (x(2),y(2)), and (x(3),x(4)) and just the estimate of F at
% the point R = (x(4),y(4)). We will use the following 10-parameter
% interpolant to do the job:
%
% c(x,y) = a(1) + a(2)(x-x0) + a(3)(y-y0) +
% a(4)(x-x0)^2 + a(5)(x-x0)(y-y0) + a(6)(y-y0)^2 +
% a(7)(x-x0)^3 + a(8)(x-x0)^2(y-y0) + a(9)(x-x0)(y-y0)^2 +a(10)(y-y0)^3
%
% This will lead to a 10-by-10 linear system M*a = rhs.
M = zeros(10,10);
rhs = zeros(10,1);
% We first assemble the equations associated with interpolating the
% F-values...

M(1:4,:) = [ones(4,1) x-x0 y-y0 ...
(x-x0).^2 (x-x0).*(y-y0) (y-y0).^2 ...
(x-x0).^3 (x-x0).^2.*(y-y0) (x-x0).*(y-y0).^2 (y-y0).^3];
rhs(1:4) = [F{1}.fvalue; F{2}.fvalue; F{3}.fvalue; fR];

% Noting that the partial of c(x,y) with respect to x is
%
% cx(x,y) = a(2)(x-x0) + 2*a(4)(x-x0) + a(5)(y-y0) +
% 3*a(7)(x-x0)^2 + 2*a(8)(x-x0)(y-y0) + a(9)(y-y0)^2
%
% we assemble the equations associated with interpolating the Fx values at
% the vertices..
M(5:7,:) = [zeros(3,1) ones(3,1) zeros(3,1) ...
2*(x(1:3)-x0) (y(1:3)-y0) zeros(3,1) ...
3*(x(1:3)-x0).^2 2*(x(1:3)-x0).*(y(1:3)-y0) 3*(y(1:3)-y0).^2 zeros(3,1)];
rhs(5:7) = [F{1}.fx;F{2}.fx ; F{3}.fx];

% Noting that the partial of c(x,y) with respect to y is
%
% cy(x,y) = a(3)(y-y0) + a(5)(x-x0) + 2*a(6)(y-y0) +
% a(8)(x-x0)^2 + 2*a(9)(x-x0)(y-y0) +3*a(10)(y-y0)^2
%
% we assemble the equations associated with interpolating the Fx values at
% the vertices..

M(8:10,:) = [zeros(3,1) zeros(3,1) ones(3,1) ...
zeros(3,1) (x(1:3)-x0) 2*(y(1:3)-y0) ...
zeros(3,1) (x(1:3)-x0).^2 2*(x(1:3)-x0).*(y(1:3)-y0) 3*(y(1:3)-y0).^2];

rhs(8:10) = [F{1}.fy;F{2}.fy;F{3}.fy];
% Solve the linear system
a = M\rhs;
% Note that we estimate F and Q via c(x0,y0) = a(1)...
z = a(1);

% This procedure is O(h^4) based on the test script.
% -1 or -2 if your "10th condition" leads to O(h^2) accuracy.
% -1 for sparse commenting.