function C6()
% Tests Partition
N =63;
V1 = MakePoint(4,-2);
V2 = MakePoint(2,4);
V3 = MakePoint(-7,1);
S = MakeTriangle(V1,V2,V3);
close all
figure
PlotTriangle(S)
axis equal off
hold on
T = PartitionCVL(S,N);
for k=1:N
PlotTriangle(T(k))
pause(.1)
end
hold off
% Check for equal arae
goal = Area(S)/N;
err = 0;
for k=1:N
err = max(err,abs(goal-Area(T(k))));
end
title(sprintf('N = %1d max|Area(T(k)) - Area(S)/N| = %7.3e',N,err),'fontsize',18)

function T = PartitionCVL(S,N)
% S is a triangle and N is a positive integer.
% T is a structure array of triangles with the property
% that T(1),...,T(N) partition S and each has the same area.
if N==1
% Nothing to do--the base case
T = S;
else
% Roughly halve N...
N1 = floor(N/2);
N2 = N-N1;
mu = N1/(N1+N2);
% Split S into S1 and S2 where the the Area(S1) = (N1/N)*Area(S)
% and Area(S2) = (N2/N)*Area(S)..
% Since we "bisect" the largest vertex angle, need the side lengths..
a = sqrt((S.P1.x-S.P2.x)^2+(S.P1.y-S.P2.y)^2);
b = sqrt((S.P2.x-S.P3.x)^2+(S.P2.y-S.P3.y)^2);
c = sqrt((S.P3.x-S.P1.x)^2+(S.P3.y-S.P1.y)^2);
if a>=b && a>=c
% Side opposite P3 is the longest
Q = MakePoint((1-mu)*S.P1.x+mu*S.P2.x,(1-mu)*S.P1.y+mu*S.P2.y);
S1 = MakeTriangle(Q,S.P3,S.P1);
S2 = MakeTriangle(Q,S.P3,S.P2);
elseif b>=a && b>=c
% Side opposite P1 is the longest
Q = MakePoint((1-mu)*S.P2.x+mu*S.P3.x,(1-mu)*S.P2.y+mu*S.P3.y);
S1 = MakeTriangle(Q,S.P1,S.P2);
S2 = MakeTriangle(Q,S.P1,S.P3);
else
% Side opposite P2 is the longest
Q = MakePoint((1-mu)*S.P3.x+mu*S.P1.x,(1-mu)*S.P3.y+mu*S.P1.y);
S1 = MakeTriangle(Q,S.P2,S.P3);
S2 = MakeTriangle(Q,S.P2,S.P1);
end
T1 = PartitionCVL(S1,N1);
T2 = PartitionCVL(S2,N2);
T = T1;
for k = length(T1)+1:length(T1)+length(T2)
T(k) = T2(k-length(T1));
end
end

function PlotTriangle(T)
% Displays triangle T in th ecurrent figure window.
% Assumes that hold is on
plot([T.P1.x T.P2.x T.P3.x T.P1.x],[T.P1.y T.P2.y T.P3.y T.P1.y],'k','linewidth',2)

function A = Area(T)
% A is the area of triangle T
a = sqrt((T.P1.x-T.P2.x)^2+(T.P1.y-T.P2.y)^2);
b = sqrt((T.P2.x-T.P3.x)^2+(T.P2.y-T.P3.y)^2);
c = sqrt((T.P3.x-T.P1.x)^2+(T.P3.y-T.P1.y)^2);
s = (a+b+c)/2;
A = sqrt(s*(s-a)*(s-b)*(s-c));

function P = MakePoint(x,y)
% P is a 2-field structure with the property that
% (P.x,P.y) = (x,y).
P = struct('x',x,'y',y);

function T = MakeTriangle(V1,V2,V3)
% V1, V2, and V3 are points and T is a triangle that has
% those points as vertices.
T = struct('P1',V1,'P2',V2,'P3',V3);