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