# How to simulate the true reflective movement of a particle bouncing around in an ellipse?

Please help me to simulate the movement of a particle inside a region with elliptical walls such that particle is reflected from the walls and continues to move.

A friend was able write code to simulate a particle represented by a Disk bouncing around inside a square, but we can't do it for an ellipse.

x = 0.5;
y = 0.5;

vx = 1;
vy = Pi/2;

step = 0.01;

Animate[
x = x + vx*step;
y = y + vy*step;
If[Abs[x - 1] <= radius || Abs[x] <= radius , vx = -vx];
If[Abs[y - 1] <= radius || Abs[y] <= radius, vy = -vy];
Graphics[{
Cyan, Rectangle[{0, 0}, {1, 1}],
Point[{0.0, 0.0}], Point[{1.0, 1.0}]
}],
{t, 0, Infinity}
]


## Edit V10!

This is simple example what we can now do in real time!

R = RegionUnion @@ Table[Disk[{Cos[i], Sin[i]}, .4], {i, 0, 2 Pi, Pi/6.}];
R2 = RegionBoundary@DiscretizeRegion@R;

go[] := (While[r > .105, x += v; r = RegionDistance[R2, x]; Pause[.01]]; bounce[];)

bounce[] := With[{normal = Normalize[x - RegionNearest[R2, x]]},
If[break, Abort[]];
v = .01 Normalize[v - 2 v.normal normal];
x = x + v;
r = RegionDistance[R2, x]; go[]
]

x = {1, 0.};
pos = {x};
break = False;
v = .01 Normalize@{2, 1.};
r = RegionDistance[R2, x];

RegionPlot[R2, Epilog -> Dynamic@Disk[x, .1], AspectRatio -> Automatic]
Button["break at edge", break = True;]
go[]


This is an example, not perfect but nice enough to start.

## V9

Unfortunately I don't have time to explain now. But take a look at wikipedia ellips site, tangent line part especially.

DynamicModule[{u = 0, t0, imp, v1, x0 = {0, .49}, v0 = {.5, -1.0}, t, a = 1, b = .5,
c, f1, f2},
DynamicWrapper[
Graphics[{ Thick, Scale[Circle[], {a, b}], AbsolutePointSize@7, Dynamic@Point[x0],
Dashed, Thin, Dynamic@Line[{{x0, imp}, {imp, imp + Normalize@v1},
{imp - normal, imp + normal}}]
}, PlotRange -> 1.1, ImageSize -> 500, Frame -> True],
Refresh[
If[(#/a)^2 + (#2/b)^2 & @@ x0 < 1,
x0 += v0;,
x0 = imp + v1; v0 = v1; rec]
, TrackedSymbols :> {}, UpdateInterval -> .001]]
,
Initialization :> (
c = Sqrt[a^2 - b^2]; v0 = Normalize[v0]/100; f1 = {-c, 0}; f2 = {c, 0};

rec := ({t0, imp} = {t, x0 + t v0
} /. Quiet@NSolve[(#/a)^2 + (#2/b)^2 & @@ (x0 + t v0) == 1. &&
t > 0, t, Reals][[1]];
normal = Normalize[Normalize[imp - f1] + Normalize[imp - f2]];

v1 = Normalize[v0 - 2 normal (v0.normal)]/100;(*bounce*));

rec)]


• @Kuba, what is imp ? Dec 18, 2013 at 21:31
• @Spizhen imp = x0 + t v0 where t is from NSolve. So it is next impact position.
– Kuba
Dec 18, 2013 at 21:33
• @Kuba, and what happens if the a space between the variables simply placed? Here, for example imp = x0 + t v0 I do not know about this. Maybe it's a feature version 9? Dec 18, 2013 at 21:44
• @Kuba, space between t and v0. what is it? Dec 19, 2013 at 6:22
• @Spizhen t is the parameter, time if you want, x0 and v0 are vectors, of position and velocity. imp = x0 + t v0 means that starting from x0 with velocity v0 the ball will hit the border after t time.
– Kuba
Dec 19, 2013 at 6:57

My goal was quite ambitious. I wanted to create a way to let any rigid body bounce elastically against any other surface in MMA V9. To do this I use "masks" for the object and the environment. These masks are black and white images. White indicates that this is where the object/surface is, black is empty space. I can calculate the overlap between the object and the surface using Mathematica's image processing functions. Using the overlap I can calculate the normal of the surface. After that it's simple physics to change the velocity of the object accordingly. The code looks like this:

obj[mask_] := Graphics[{
},
PlotRange -> {{0, 500}, {0, 500}},
ImageSize -> {500, 500},
Background -> Black
]

forceVector[obj_, env_, center_] := N@Normalize[Plus @@ (center - # & /@ PixelValuePositions[ImageMultiply[obj, env], 1])]

step[{pt_, v_}] := Module[{f, nv},
f = forceVector[obj[Disk[pt, 20]], ColorNegate@obj[Disk[{250, 250}, {100, 200}]], pt] /. (0. -> {0, 0});
nv = If[v.f < 0, v - 2 v.f f, v];
{pt + nv, nv}
]

pts = NestList[step, {{250, 250}, {1, 2}}, 1000];

frames = Graphics[{
Black, Rectangle[{0, 0}, {500, 500}],
White, Disk[{250, 250}, {100, 200}],
Orange, Disk[#, 20]
},
PlotRange -> {{0, 500}, {0, 500}},
ImageSize -> {500, 500}
] & /@ pts[[All, 1]];

ListAnimate[frames]


Here's a gif with a reduced number of frames:

One can play with the velocity of the disk as well as the number of frames to get a longer path without as many calculations. This method is not very fast.

If you don't have the time/computing power to pre-compute the position list, you can still view the simulation using the code below. It will probably be very slow on many computers though (which is why I chose to pre-compute positions):

DynamicModule[{pt = {250, 250}, v = {6, 2}, f},
Dynamic[
f = forceVector[obj[Disk[pt, 20]],
ColorNegate@obj[Disk[{250, 250}, {100, 200}]],
pt] /. (0. -> {0, 0});
If[v.f < 0, v = v - 2 v.f f];
pt = pt + v;
Graphics[{
Black, Rectangle[{0, 0}, {500, 500}],
White, Disk[{250, 250}, {100, 200}],
Orange, Disk[pt, 20]
},
PlotRange -> {{0, 500}, {0, 500}},
ImageSize -> {500, 500}
]
]
]

• Your method is very long build. Even if the decrease in the value of two times. Is this normal? How long should I wait? 10 minutes is not ready. P.S. Hardware: MacBook Air 13" (mid 2013) Dec 17, 2013 at 20:13
• @Spizhen Yeah, it's kind of slow because the way it detects collision is very expensive. I use an iMac and it certainly didn't take me ten minutes to execute that code but it may very well be the case on Macbook Air. Dec 17, 2013 at 20:23
• whether there is a primitive example of solving tasks? Maybe a circle instead of an ellipse. Tomorrow morning need to show my teacher =( Dec 17, 2013 at 20:32
• @Spizhen whoah, aren't you ashamed writing such comments? have you tried anything? C. E., interesting approach ;) +1.
– Kuba
Dec 17, 2013 at 20:35
• @Spizhen I added a version that you can show without having to do anything in advance. You still need the function definition for forceVector and obj from earlier. Dec 18, 2013 at 0:37
• @Kuba has provided an excellent solution. Here we follow his idea and use another approach like WhenEvent to get the particle tracing.

• r[t] is the trace of the point.

Needs["OpenCascadeLink"];
SeedRandom[1];
reflect[vector_,
normal_] = -(vector - 2 (vector - Projection[vector, normal])) //
Simplify;
RegionUnion @@
Table[Ball[{Cos[i], Sin[i], 0}, .4], {i, 0, 2 Pi, Pi/6.}]];
"ShapeSurfaceMeshOptions" -> {"AngularDeflection" -> .1}];
R = BoundaryMeshRegion[bm];
R2 = RegionBoundary[R];
dist = RegionDistance[R2];
proj = RegionNearest[R2];
pt0 = RandomPoint[R, 1][[1]];
v0 = {1., 1., 0.};
d0 = 0.01*Norm[v0];
sol = NDSolveValue[{r''[t] == {0, 0, 0}, r[0] == pt0, r'[0] == v0,
WhenEvent[dist@r[t] <= d0,
r'[t] -> reflect[r'[t], r[t] - proj@r[t]]]}, r[t], {t, 0, 100},
MaxStepSize -> 0.01];
ani = Animate[
Show[Graphics3D[{{FaceForm[Opacity[.2], Darker@Cyan], EdgeForm[],
R2}}, Boxed -> False],
ParametricPlot3D[sol, {t, 0, c}, Mesh -> {{c}},
MeshStyle -> {AbsolutePointSize[8], Red},
Method -> {"BoundaryOffset" -> False},
PlotStyle -> {Opacity[.9], White}, PlotPoints -> 400,
PerformanceGoal -> "Quality", PlotRange -> All] /. Line -> Arrow,
ViewPoint -> {1, 1, .8}, Background -> Gray], {c, $MachineEpsilon, 100}, AnimationRate -> 1, ControlPlacement -> Bottom]  • Test another type of bounce and another obstacle surfaces. SeedRandom[1]; reflect[vector_, normal_] = -(vector - 2 (vector - Projection[vector, normal])) // Simplify; surf = ExampleData[{"Geometry3D", "UtahTeapot"}, "Region"]; cube = TransformedRegion[Cuboid @@ Transpose@RegionBounds[surf], ScalingTransform[{1.2, 1.2, 1.2}, RegionCentroid[surf]]] // DiscretizeRegion // RegionBoundary; reg = RegionUnion[surf, cube]; dist = RegionDistance[reg]; proj = RegionNearest[reg]; pt0 = RegionCentroid[reg]; v0 = {1., 1., 0.2}; d0 = 0.01*Norm[v0]; sol = NDSolveValue[{r''[t] == {0, 0, -9.8}, r[0] == pt0, r'[0] == v0, WhenEvent[dist@r[t] <= d0, r'[t] -> reflect[r'[t], r[t] - proj@r[t]]]}, r[t], {t, 0, 20}, MaxStepSize -> 0.001] ani = Animate[ Show[Graphics3D[{{FaceForm[Opacity[.2], Cyan], EdgeForm[], surf}, {FaceForm[Opacity[.3], LightGray], EdgeForm[], cube}}, Boxed -> False], ParametricPlot3D[sol, {t, 0, c}, Mesh -> {{c}}, MeshStyle -> {AbsolutePointSize[8], Red}, Method -> {"BoundaryOffset" -> False}, PlotStyle -> {Opacity[.9], White}, PlotPoints -> 400, PerformanceGoal -> "Quality", PlotRange -> All] /. Line -> Arrow, Background -> LightGray], {c,$MachineEpsilon, 20},
AnimationRate -> 1, ControlPlacement -> Bottom]
`