# The Generalized reflectionrestart;E1 := solve({V[y] - Q[y] = k[3]*(V[x]-Q[x]),
V[y] - T[y] = k[1]*(V[x] - T[x])}, {V[x], V[y]}):
assign(E1):
E2 := solve({C[y] - P[y] = k[2]*(C[x] - P[x]),
C[y] - T[y] = k[1]*(C[x] - T[x])}, {C[x], C[y]}):
assign(E2):M[x] := P[x] + 2*(C[x] - P[x]): M[y] := P[y] + 2*(C[y] - P[y]):
k[3] := (M[y] - Q[y])/(M[x] - Q[x]):
T[y] := Y(t): T[x] := X(t):
k[2] := -1/k[1]:
k[1] := diff(Y(t), t)/diff(X(t),t):
Aimerr := normal((T[x] - V[x])^2 + (T[y] - V[y])^2):DA := op(1, numer(Aimerr)):
DB := denom(Aimerr):
DA1 := factor(select(has, DA, {diff(X(t), t)^2,
diff(Y(t), t)^2})):
Da2 := select(has, DA, diff(X(t), t)):
DA2 := factor(select(has, Da2, diff(Y(t), t))):
DA - expand(DA2 + DA1);DAA := DA1 + DA2;# The Shortest Trajectory MethodLP := sqrt((X(t) - P[x])^2 + (Y(t) - P[y])^2):
LQ := sqrt((X(t) - Q[x])^2 + (Y(t) - Q[y])^2):
lP := diff(LP,t): lQ := diff(LQ,t):
dll := numer(lP)*denom(lQ) + numer(lQ)*denom(lP);# The Circular Billiard TableX(t) := cos(t): Y(t) := sin(t): Q[y]:= 0:DA := simplify(DAA);SolDA := solve(DA,t);P[x] := 1/2: P[y] := 1/2: Q[x] := -3/5:
DA := simplify(DA):
SolDA := solve(DA, t):FSolDA := evalf(allvalues([SolDA])); for i from 1 by 1 to nops(FSolDA) do:
tA[i] := FSolDA[i];
DBi := evalf(subs(t = tA[i], DB));
if DBi = 0 then print(`Root`[i]*` unusable!`); fi;
print('i' = i, ` `*'tA' = tA[i],` ==> `*'DB' = DBi);
od:P[x] := 'P[x]': P[y] := 'P[y]': Q[x] := 'Q[x]': Q[y] := 0:
dl := simplify(dll);# Soldl := solve(dl,t):eq := op(1, dl) + op(2, dl) = -op(3, dl):
eq := expand(map(u -> u^2, eq)):
eq := lhs(eq) - rhs(eq) = 0:
Soldl := solve(eq, t):P[x] := 1/2: P[y] := 1/2: Q[x] := -3/5:
eq;Soldl := solve(eq,t):FSoldl := evalf(allvalues([Soldl]));DAd := diff(DA, t): dld := diff(dl, t):
tamax := fsolve(DAd = 0, t, 6..7):
tlmax := fsolve(dld = 0, t, 6..7):
Dan := DA/abs(subs(t = tamax, DA)):
dln := dl/abs(subs(t = tlmax, dl)):
with(plots):
p1 := plot(Dan, t = 0..2*Pi, linestyle=3):
p2 := plot(dln, t = 0..2*Pi, linestyle=0):
display({p1,p2});plot([Dan, dln, t = 0..2*Pi]);read `ResultPlot.map`:
ResultPlot(P[x], P[y], Q[x], Q[y], X(t), Y(t),
Dan, 0, 2*Pi, tA, 4, -1..2, -1.25..1); # The Elliptical Billiard TableX(t) := cos(t): Y(t) := 4/5*sin(t):
P[x] := -3/5: P[y] := -3/10: Q[x] := 1/2: Q[y] := 1/2:
DA := simplify(DAA);SolDA := [solve(DA,t)]: PZ6:=op(1,selectfun(SolDA[1],RootOf)[]);Index:=op(2,selectfun(SolDA[1],RootOf)[]);SolDa:=[subs(PZ6='PZ6',Index='i',SolDA[1]),`i=1..6`];irreduc(PZ6);tA[1] := fsolve(DA, t, -3..-2);
tA[2] := fsolve(DA, t, -1..0);
tA[3] := fsolve(DA, t, 1..2);
tA[4] := fsolve(DA, t, 2..3);
tl[1] := fsolve(dll, t, -3..-2);
tl[2] := fsolve(dll, t, -1..0);
tl[3] := fsolve(dll, t, 1..2);
tl[4] := fsolve(dll, t, 2..3);Dan := DA/abs(evalf(subs( t = fsolve(diff(DA, t),
t, -2..-1), DA))):
dln := dll/abs(evalf(subs( t = fsolve(diff(dll, t),
t, -2..-1), dll))):
p1:=plot(Dan, t = 0..2*Pi,linestyle=0):
p2:=plot(dln, t = 0..2*Pi,linestyle=3):
display({p1,p2});
ResultPlot(P[x], P[y], Q[x], Q[y], X(t), Y(t),
Dan, -Pi, Pi, tA, 4, -1.5..2, -0.8..0.8);# The Snail Billiard TableX(t) := t*cos(t): Y(t) := t*sin(t):P[x] := -6: P[y] := -12: Q[x] := 7: Q[y] := 5:DA := simplify(DAA):Dan := DA/abs(evalf(subs(t = fsolve(diff(DA, t),
t, 15..17), DA))): dln := dll/abs(evalf(subs(t = fsolve(diff(dll, t),
t, 15..17), dll))): p1:=plot(Dan, t = 4*Pi..6*Pi,linestyle=0): p2:=plot(dln, t = 4*Pi..6*Pi,linestyle=3):display({p1,p2});p1:=plot(Dan, t = 12.5..12.85,linestyle=0):
p2:=plot(dln, t = 12.5..12.85,linestyle=3):
display({p1,p2});tA[1] := fsolve(DAA, t, 12.5..12.6);
tA[2] := fsolve(DAA, t, 12.7..12.8);
tA[3] := fsolve(DAA, t, 14..15);
tA[4] := fsolve(DAA, t, 16..17);ResultPlot(P[x], P[y], Q[x], Q[y], X(t), Y(t),
dln, 4*Pi, 6*Pi, tA, 4, -20..30, -31..26);# The Star Billiard TableX(t) := cos(t)*(1 + sin(5*t)/5):
Y(t) := sin(t)*(1 + sin(5*t)/5):
P[x] := -4/5: P[y] := 1/5: Q[x] := 3/5: Q[y] := -4/5:
plot({DAA, dll}, t = 0..2*Pi);Dan := DAA/evalf(abs(subs(t = fsolve(diff(DAA, t),
t, 1..1.5), DAA))): dln := dll/evalf(abs(subs(t = fsolve(diff(dll, t),
t, 1.75..1.85), dll))):p1:=plot(Dan,t = 0..2*Pi,linestyle=0):
p2:=plot(dln,t = 0..2*Pi,linestyle=3):
display({p1,p2});
plot([Dan, dln, t = 0..2*Pi]);tA[1] := fsolve(DAA = 0, t, 0.3..0.4);
tA[2] := fsolve(DAA = 0, t, 0.9..1.0);
tA[3] := fsolve(DAA = 0, t, 1.5..1.6);
tA[4] := fsolve(DAA = 0, t, 2.3..2.4);
tA[5] := fsolve(DAA = 0, t, 3.3..3.5);
tA[7] := fsolve(DAA = 0, t, 4.0..4.2);
tA[8] := fsolve(DAA = 0, t, 4.7..4.9);
tA[9] := fsolve(DAA = 0, t, 5.3..5.5);
tA[10] := fsolve(DAA = 0, t, 5.5..5.7);ResultPlot(P[x], P[y], Q[x], Q[y], X(t), Y(t),
Dan, 0, 2*Pi, tA, 10, -1.4..1.8, -1.6..2.2);T1[x] := evalf(subs(t = tA[1], X(t))):
T1[y] := evalf(subs(t = tA[1], Y(t))):
eq1 := X(t) = Q[x] + u*(T1[x] - Q[x]):
eq2 := Y(t) = Q[y] + u*(T1[y] - Q[y]):
fsolve({eq1, eq2}, {u, t}, {u = 0..0.5, t=3*Pi/2..2*Pi});fsolve({eq1, eq2}, {u, t}, {u = 0.33..0.8, t = 5.82..6.2});