if get('geometry2, 'version)=false then load("geometry2");
if get('cone, 'version)=false then load("cone"); 
n:2;
/* Introduce simplifications */
u0:0;
v0:0;
u1:1;
v1:0;

/* Define points Ai=[xi,yi] */
for i:0 thru n-1 do
	A[i]:[concat(x,i),concat(y,i)];

/* Define circle centers Oi=[ui,vi] */
for i:0 thru n-1 do
	O[i]:[concat(u,i),concat(v,i)];

/* Define radii and lengths and Lagrange multipliers */
for i:0 thru n-1 do (
	R[i]:concat(R,i),
	l[i]:concat(l,i),
	p[i]:concat(p,i),
	q[i]:concat(q,i)
);

/* Ai is on the circle with center Oi and radius Ri */
eqns1:makelist(distance(O[i],A[i],R[i]),i,0,n-1);


shift(k):=if k=n-1 then 0 else k+1;
lshift(k):=if k=0 then n-1 else k-1;


/* Lagrange expression */
f: sum(l[i],i,0,n-1)-sum(p[i]*distance(O[i],A[i],r[i]),i,0,n-1)
	-sum(q[i]*distance(A[i],A[shift(i)],l[i]),i,0,n-1);

/* Lagrange equations */
eqns0:makelist(diff(f,l[i]),i,0,n-1);
eqns1:makelist(diff(f,concat(x,i)),i,0,n-1);
eqns2:makelist(diff(f,concat(y,i)),i,0,n-1);
eqns3:makelist(diff(f,p[i]),i,0,n-1);
eqns4:makelist(diff(f,q[i]),i,0,n-1);

eqns:expand(append(eqns0,eqns1,eqns2,eqns3,eqns4));

vars:listofvars(makelist([p[i],q[i],l[i]],i,0,n-1));
params:listofvars(makelist([A[i],O[i],r[i]],i,0,n-1));
cover:expand([[[],[]]]);

eqns:''eqns;
stringout("billiard3.out",eqns,vars,params,cover);