3
$\begingroup$

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}]}]]

output

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:

enter image description here

|improve this question
$\endgroup$
  • 2
    $\begingroup$ related: 38927 $\endgroup$ – Kuba♦ Mar 12 '15 at 20:21
  • $\begingroup$ Also related (and unanswered): (73512) $\endgroup$ – Mr.Wizard♦ Mar 12 '15 at 21:03
  • $\begingroup$ @Kuba: I looked at your code, and it definitely helps. But it is different than what I'm trying to do: I'd like a list of points (up to some length) visited by the bouncing particle, and I'd also like to do the computations for when the particle hits the boundary, not when it gets sufficiently close. $\endgroup$ – Steve D Mar 12 '15 at 22:45
  • $\begingroup$ The title of the question contains "in ellipses" - but that's misleading if you want to allow exterior reflections. $\endgroup$ – Jens Mar 13 '15 at 3:47
  • $\begingroup$ Also related (but still just one curve): 63690 $\endgroup$ – Michael E2 Mar 13 '15 at 3:48
5
$\begingroup$

Update

The original approach for one ellipse (below) may be adapted for several:

eq = With[{p = {x, y} - {x0, y0}},
    (RotationMatrix[-t0].p).{{1/a^2, 0}, {0, 1/b^2}}.(RotationMatrix[-t0].p)] - 1;
sub[pt_] := Thread[{x, y} -> pt];
dir[t0_] := {Cos[t0], Sin[t0]};

ClearAll[next, cuts];
Block[{a, b, t0, x0, y0, x1, y1, α, x, y, t, ellipses}, 
  next = With[{sol = Simplify[
        t /. Solve[eq == 0 /. sub[{x1, y1} + t dir[α]], t], 
        TimeConstraint -> 0.1]}, 
    With[{x2 = x /. sub[{x1, y1} + t dir[α]], 
      y2 = y /. sub[{x1, y1} + t dir[α]]}, 
     With[{θ = ArcTan[D[eq, x], D[eq, y]] /. sub[{x2, y2}]},
      (* definitions using above algebra *)
      cuts[sect_][ell_] := sect /. Thread[{a, b, t0, x0, y0} -> ell];
      Function @@ Hold[{ellipses, x1, y1, α},
        With[{sects = cuts[sol] /@ ellipses},
         t = First@Sort@Select[Chop@Flatten[sects], Positive]; {x2, y2, 
           Mod[2 θ - α + Pi, 2 Pi]} /. 
          Thread[{a, b, t0, x0, y0} -> ellipses~Part~First@FirstPosition[sects, t]
         ]
        ]
      ]]]
  ];

OP's new example:

centers = {{0, 0}, {3, 0}, {0, 0}};
radii = {{1, 1}, {2, 1}, {6, 6}};
angles = {0, Pi/3, 0};
ells = Flatten[{radii, List /@ angles, centers}, {{2}, {1, 3}}];

With[{dir = 0, pt = {1.6, 0.5}}, 
  ptsdir = NestList[
    next[ells, Sequence @@ #] &, {Sequence @@ pt, dir}, 19];
  pts = ptsdir[[All, 1 ;; 2]]; 
  Graphics[{Red, Thick, Line[pts], Blue, 
    MapThread[
     GeometricTransformation[Circle[#1, #2], 
       RotationTransform[#3, #1]] &, {centers, radii, angles}]}]
  ] // AbsoluteTiming

Mathematica graphics

With 100 points:

Mathematica graphics

Original answer

I used a cartesian equation of a rotated ellipse and the angle of the normal to compute the reflection. Replace Function by Compile if you want more speed, but the figure below is computed in a little over 0.05 sec, most of which time was spent computing the plot of the ellipse.

Given the ellipse x^2/a^2 + y^2/b^2 == 1 rotated by an angle t0, an initial point {x0, y0}, and a direction α, the function next returns the next intersection and direction in a list {x1, y1, α1}. 

eq = (RotationMatrix[-t0].{x, y}).{{1/a^2, 0}, {0, 1/b^2}}.(RotationMatrix[-t0].{x, y}) - 1;
sub[pt_] := Thread[{x, y} -> pt];
dir[t0_] := {Cos[t0], Sin[t0]};

Block[{a, b, t0, x0, y0, α, x, y, t}, 
  next = (* preliminary algebra *)
   With[{sol = t /. Solve[eq == 0 /. sub[{x0, y0} + t dir[α]], t ] // Simplify},
    With[{x1 = x /. sub[{x0, y0} + t dir[α]], 
      y1 = y /. sub[{x0, y0} + t dir[α]]},
     With[{θ = ArcTan[D[eq, x], D[eq, y]] /. sub[{x1, y1}]}, (* angle of normal *)
       (* function definition *)
      Function @@ Hold[
        {a, b, t0, x0, y0, α},
        t = First@Select[sol, # > 1.*^-8 &];       (* tolerance could be ~1.*^-14 *)
        {x1, y1, Mod[2 θ - α + Pi, 2 Pi]}          (* new x, y, reflected angle *)
        ]
      ]]]
  ];

Here ptsdir contains a list of {x, y, theta} and pts contains a list of the points.

Block[{a = 4, b = 2, t0 = Pi/6},
 ptsdir = NestList[next[a, b, t0, Sequence @@ #] &, {2., 2., 0.}, 170];
 pts = ptsdir[[All, 1 ;; 2]]; 
 ContourPlot[eq == 0, {x, -5, 5}, {y, -5, 5},
  Epilog -> {Red, Thickness[0.001], Line[pts]}]
 ]
|improve this answer
$\endgroup$
  • $\begingroup$ I think the example I gave in the question makes it somewhat unclear. But I have several ellipses that need to be reflected off of (potentially) . $\endgroup$ – Steve D Mar 13 '15 at 3:27
  • $\begingroup$ @SteveD Yeah, I saw that in your comment to Jens' answer, but I had already solved it this way. I think you could adapt this approach to pick the nearest intersection (smallest positive t). The equations are all quadratic, so you can solve them beforehand as I did the one. Then knowing which ellipse was intercepted, you could compute the reflected angle. Given the image in your question, it's not clear what sort of configurations we'd have to worry about. $\endgroup$ – Michael E2 Mar 13 '15 at 3:38
  • $\begingroup$ @Jens I suppose my answer will work on any algebraic equation for which the solutions returned by Solve will always contain exactly one positive solution (and for which the gradient points outward). Otherwise, one could use t = First@Sort@Select[Chop@sol, Quiet@TrueQ[# > 0] &]. Still, I wouldn't expect it to work as well for complicated equations. (+1 to you, too, already.) $\endgroup$ – Michael E2 Mar 13 '15 at 4:28
  • $\begingroup$ Wow, that is a serious speed improvement! I'm going to accept this answer because it works in the generality I wanted, but I also would like to point out that I've noticed my NestWhileList was actually adding to the slow performance of my code. And I noticed that because you used NestList instead. Originally I had been stopping the line when it passed a certain length threshold, but it does seem to slow things up. Anyway, thanks again! $\endgroup$ – Steve D Mar 14 '15 at 22:14
0
$\begingroup$

I changed the following functions based on what @Kuba did in his related post:

maxTime = 50; centers = {{0, 0}, {3, 0}};
radii = {{1, 1}, {2, 1}};
angles = {0, Pi/3};
circs = MapThread[
   TransformedRegion[Circle[#1, #2], 
     RotationTransform[#3, #1]] &, {centers, radii, angles}];
r = RegionUnion @@ circs;
region = RegionBoundary@DiscretizeRegion@r;
ClearAll[x, y, circleIntersection, nextDirection, nextPt, whichCirc, \
tangentVector]
circleIntersection[direction_, pt_List] := 
  p /. Quiet@
    NSolve[p \[Element] 
      RegionIntersection[r, 
       Line[{pt, 
         pt + {Cos[direction], Sin[direction]}*{maxTime, maxTime}}]], 
     p];

and

nextDirection[direction_, pt_List, ptnow_List] := 
  Module[{vec, v = {Cos[direction], Sin[direction]}, 
    normal = Normalize[pt - RegionNearest[region, pt]]}, 
   vec = Normalize[v - 2 v.normal normal] + ptnow; 
   If[pt == ptnow, direction, ArcTan[First@vec, Last@vec]]];

I saw about a 5x speed-up versus my original code, but something is now wrong: the trajectories are being computed incorrectly (something is off in the tangent vectors). Guidance is most appreciated.

|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.