[1] | 1 | if get('geometry2, 'version)=false then load("geometry2");
|
---|
| 2 | n:2;
|
---|
| 3 | /* Introduce simplifications */
|
---|
| 4 | u0:0;
|
---|
| 5 | v0:0;
|
---|
| 6 | u1:1;
|
---|
| 7 | v1:0;
|
---|
| 8 |
|
---|
| 9 | /* Define points Ai=[xi,yi] */
|
---|
| 10 | for i:0 thru n-1 do
|
---|
| 11 | A[i]:[concat(x,i),concat(y,i)];
|
---|
| 12 |
|
---|
| 13 | /* Define circle centers Oi=[ui,vi] */
|
---|
| 14 | for i:0 thru n-1 do
|
---|
| 15 | O[i]:[concat(u,i),concat(v,i)];
|
---|
| 16 |
|
---|
| 17 | /* Define radii and lengths */
|
---|
| 18 | for i:0 thru n-1 do (
|
---|
| 19 | R[i]:concat(R,i),
|
---|
| 20 | l[i]:concat(l,i),
|
---|
| 21 | p[i]:concat(p,i),
|
---|
| 22 | q[i]:concat(q,i)
|
---|
| 23 | );
|
---|
| 24 |
|
---|
| 25 | /* Ai is on the circle with center Oi and radius Ri */
|
---|
| 26 | eqns1:makelist(distance(O[i],A[i],R[i]),i,0,n-1);
|
---|
| 27 |
|
---|
| 28 |
|
---|
| 29 | shift(k):=if k=n-1 then 0 else k+1;
|
---|
| 30 | lshift(k):=if k=0 then n-1 else k-1;
|
---|
| 31 |
|
---|
| 32 |
|
---|
| 33 | /* Lagrange expression */
|
---|
| 34 | f: sum(l[i],i,0,n-1)-sum(p[i]*distance(O[i],A[i],r[i]),i,0,n-1)
|
---|
| 35 | -sum(q[i]*distance(A[i],A[shift(i)],l[i]),i,0,n-1);
|
---|
| 36 |
|
---|
| 37 | /* Lagrange equations */
|
---|
| 38 | eqns0:makelist(diff(f,l[i]),i,0,n-1);
|
---|
| 39 | eqns1:makelist(diff(f,concat(x,i)),i,0,n-1);
|
---|
| 40 | eqns2:makelist(diff(f,concat(y,i)),i,0,n-1);
|
---|
| 41 | eqns3:makelist(diff(f,p[i]),i,0,n-1);
|
---|
| 42 | eqns4:makelist(diff(f,q[i]),i,0,n-1);
|
---|
| 43 |
|
---|
| 44 | /* The rationale is that the entire circles would have to be orbits
|
---|
| 45 | and thus we are seeking for such parameters for which the circle equations
|
---|
| 46 | imply satisfaction of critical point equations */
|
---|
| 47 |
|
---|
| 48 | eqns:expand(append(eqns0,eqns1,eqns2));
|
---|
| 49 |
|
---|
| 50 |
|
---|
| 51 | vars:listofvars(makelist([p[i],q[i],l[i]],i,0,n-1));
|
---|
| 52 | params:listofvars(makelist([A[i],O[i],r[i]],i,0,n-1));
|
---|
| 53 | cover:expand([append(''eqns3,''eqns4),[]]);
|
---|
| 54 |
|
---|
| 55 | eqns:''eqns;
|
---|
| 56 | stringout("billiard2.out",eqns,vars,params,cover);
|
---|