I have the following code which calculates a list of points that a line would take after bouncing off various ellipses.
maxTime = 50; centers = {{0, 0}, {3, 0}};
radii = {{1, 1}, {2, 1}};
angles = {0, Pi/3};
ClearAll[x, y, circleIntersection, nextDirection, nextPt, whichCirc, \
tangentVector]
circleIntersection[direction_, pt_List,
circs_List] := {x, y} //. {ToRules[
Quiet@Reduce@
RegionMember[
RegionIntersection[RegionUnion[circs],
Line[{pt,
pt + {Cos[direction], Sin[direction]}*{maxTime,
maxTime}}]], {x, y}]]};
nextPt[direction_, pt_List] :=
With[{val =
First@PadLeft[
MinimalBy[EuclideanDistance[pt, #] &]@
Function[x, Select[x, Abs[EuclideanDistance[pt, #]] > .05 &]]@
N@circleIntersection[direction, pt,
MapThread[
TransformedRegion[Circle[#1, #2],
RotationTransform[#3, #1]] &, {centers, radii,
angles}]], 1, pt]},
If[ListQ[val], val,
pt + {Cos[direction], Sin[direction]}*{maxTime, maxTime}]];
whichCirc[direction_, pt_List] :=
First@PadLeft[
Flatten@Position[
RegionMember[#, nextPt[direction, pt]] & /@
MapThread[
TransformedRegion[Circle[#1, #2],
RotationTransform[#3, #1]] &, {centers, radii, angles}],
_?TrueQ, 1, 1], 1];
tangentVector[pt_List, pos_Integer] :=
Module[{a, b, c, d, e, m, n, s, r, dx, dy}, m = Cos[angles[[pos]]];
n = Sin[angles[[pos]]]; s = radii[[pos]][[2]];
r = radii[[pos]][[1]]; a = s^2*m^2 + r^2*n^2;
b = s^2*n^2 + r^2*m^2;
c = 2 (m*n*r^2*centers[[pos]][[2]] - centers[[pos]][[1]]*n^2*r^2 -
centers[[pos]][[1]]*m^2*s^2 - centers[[pos]][[2]]*m*n*s^2);
d = 2 (centers[[pos]][[1]]*m*n*r^2 - centers[[pos]][[2]]*m^2*r^2 -
centers[[pos]][[1]]*m*n*s^2 - centers[[pos]][[2]]*n^2*s^2);
e = -2 (m*n*r^2 - m*n*s^2); dx = 2*a*pt[[1]] + c + e*pt[[2]];
dy = 2*b*pt[[2]] + d + e*pt[[1]]; If[pos > 0, {-dy, dx}]];
nextDirection[direction_, pt_List, ptnow_List] :=
With[{pos = whichCirc[direction, ptnow]},
If[pos > 0 && ptnow != pt,
ArcTan[Sequence @@ (ReflectionTransform[tangentVector[pt, pos],
pt][ptnow] - pt)], direction]];
f1[{dir_, pt_, ___}] := {dir, nextPt[dir, pt], pt};
g1[{dir_, newpt_, pt_}] := {nextDirection[dir, newpt, pt], newpt, pt};
ellipsePtList[dir_, startPt_, mt_] :=
Module[{timeUpYet, len = 0},
timeUpYet[pt1_, pt2_] :=
Module[{dist = EuclideanDistance[pt1, pt2]}, len += dist;
Return[len < mt]];
Flatten[NestWhileList[g1[f1[#]] &, {dir, startPt},
timeUpYet[#1[[2]], #2[[2]]] &, 2][[All, 2 ;; 2]], 1]]
It allows me to do the following:
With[{dir = 0, pt = {0.6, 0.3}},
Graphics[{Red, Thick, Line[ellipsePtList[dir, pt, maxTime]], Blue,
MapThread[
GeometricTransformation[Circle[#1, #2],
RotationTransform[#3, #1]] &, {centers, radii, angles}]}]]
Which produces a very nice picture, but takes about 5 seconds on my computer to run. So I'm looking for ways to speed this code up. After some preliminary testing, it seems the whichCirc
function is the bottleneck (but I will take optimizations for any part of the code). Any and all suggestions are most welcome!
PS - sorry for the code formatting above, it's hard to paste in a lot of code and make it nice!
EDIT: there was some confusion over exactly what was being asked. I need this to work in the generality that there could be several ellipses, positioned in various ways. For example, changing the first few lines above to:
maxTime = 100; centers = {{0, 0}, {3, 0}, {0, 0}};
radii = {{1, 1}, {2, 1}, {6, 6}};
angles = {0, Pi/3, 0};
and then executing the same code at the bottom (with different start point):
With[{dir = 0, pt = {1.6, 0.5}},
Graphics[{Red, Thick, Line[ellipsePtList[dir, pt, maxTime]], Blue,
MapThread[
GeometricTransformation[Circle[#1, #2],
RotationTransform[#3, #1]] &, {centers, radii, angles}]}]]
gives the final picture: