source: CGBLisp/examples/billiard2.macsyma@ 1

if get('geometry2, 'version)=false then load("geometry2");
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 */
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);

/* The rationale is that the entire circles would have to be orbits
and thus we are seeking for such parameters for which the circle equations
imply satisfaction of critical point equations */

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


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([append(''eqns3,''eqns4),[]]);

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