# The Generalized reflection restart; 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 Method LP := 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 Table X(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 Table X(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 Table X(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 Table X(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});