if get('geometry2, 'version)=false then load(geometry2);
/* Points: A,B,C, X,Y,Z.
Cosines and sines of angles at A B and C are:
COSA, COSB, COSC, SINA, SINB, SINC
CXBC, CYAC, CYCA, CZAB, CZBA, CZBA.
Lengths:
	XB, XC, YA, YC, ZA, ZB
	AB, AC, BC,
	XY, XZ, YZ 
X,Y,Z are characterized by the conditions
that 
	angle XBC is 1/3 of ABC
	angle XCB is 1/3 of ACB
	angle YAC is 1/3 of BAC
	angle YCA is 1/3 of BCA
	angle ZAB is 1/3 of CAB
	angle ZBA is 1/3 of CBA
Conclusion is that XYZ form an equilateral triangle
		XY=XZ=YZ
*/

A:[0,0]; B:[1,0]; C:[u1,u2];
X:[x1,x2]; Y:[x3,x4]; Z:[x5,x6];

vars:append(listofvars([X,Y,Z]),[XB,XC,YA,YC,ZA,ZB],
	[AB,AC,BC], [XY,XZ,YZ],
	[COSA, COSB, COSC, SINA, SINB, SINC],
	[CXBC, CXCB, CYAC, CYCA, CZAB, CZBA],
	[SXBC, SXCB, SYAC, SYCA, SZAB, SZBA],
	[V]);
	
params:listofvars([A,B,C]);
cover:[[],[u2]];


eqns:[
	/* Angles: cosine law */
	XB^2+BC^2-2*XB*BC*CXBC-XC^2,
	XC^2+BC^2-2*XC*BC*CXCB-XB^2,
	YA^2+AC^2-2*YA*AC*CYAC-YC^2,
	YC^2+AC^2-2*YC*AC*CYCA-YA^2,
	ZA^2+AB^2-2*ZA*AB*CZAB-ZB^2,
	ZB^2+AB^2-2*ZB*AB*CZBA-ZA^2,
	CXBC^2+SXBC^2-1,
	CXCB^2+SXCB^2-1,
	CYAC^2+SYAC^2-1,
	CZAB^2+SZAB^2-1,
	CZBA^2+SZBA^2-1,
	/* Main angles: cosine law */
	AB^2+AC^2-2*AB*AC*COSA-BC^2,
	AB^2+BC^2-2*AB*BC*COSB-AC^2,
	AC^2+BC^2-2*AC*BC*COSC-AB^2,
	SINA^2+COSA^2-1,		
	SINB^2+COSB^2-1,		
	SINC^2+COSC^2-1,		
	/* Distances of X,Y,Z from vertices */
	distance(X,B,XB),
	distance(X,C,XC),
	distance(Y,A,YA),
	distance(Y,C,YC),
	distance(Z,A,ZA),
	distance(Z,B,ZB),
	/* Lengths of sides of main triangle */
	distance(A,B,AB),
	distance(A,C,AC),
	distance(B,C,BC),
	/* Lengths of the trisector triangle */
	distance(X,Y,XY),
	distance(X,Z,XZ),
	distance(Y,Z,YZ),
	/* Trisector relations */
	4*CXBC^3-3*CXBC-COSB,
	4*CXCB^3-3*CXCB-COSC,
	4*CYAC^3-3*CYAC-COSA,
	4*CYCA^3-3*CYCA-COSC,
	4*CZAB^3-3*CZAB-COSA,
	4*CZBA^3-3*CZBA-COSB,
	/* sine rules */
	4*CXBC^2*SXBC-3*SXBC-SINB,
	4*CXCB^2*SXCB-3*SXCB-SINC,
	4*CYAC^2*SXCB-3*SYAC-SINA,
	4*CYCA^2*SYCA-3*SYCA-SINC,
	4*CZAB^2*SZAB-3*SZAB-SINA,
	4*CZBA^2*SZBA-3*SZBA-SINB,
	/* Conclusion: XY=XZ */
	1-v*(XY-XZ) 
];

eqns:expand(eqns);
cover:expand(cover);
stringout("trisector.out",eqns,vars,params,cover);