How to reflect from an ellipse while solving a system of differential equations?

While solving a system of differential equations to plot trajectories, we can use NDSolve and WhenEvent like shown here.

I am trying to solve this for a system of equations (Hamiltons Equations) for a free particle, so that the reflection is from a surface of ellipse (like elliptical pool table). Where r and th are r and theta coordinates.

In the r, theta cordinates, the equation of ellipse reads as: The initial conditions I have given (in the variable of time i.e. r(t), theta(t), etc) are r(t=0)= 2.29 (which is r(theta=Pi)), theta(t=0) = Pi, pr(t=0)=1, pth(t=0) = 0. Where pr and pth are radial and angular momentum given by Hamilton's Equations

This is how I am implementing it, but the reflection is occurring from a circle and not an ellipse. Please help.

The code is:

b=1;
eps=0.9;
m=1;
R=b/(Sqrt[1-(eps*Cos[th])^2]);
H = pr[t]^2/(2*m) + pth[t]^2/(2*m*r[t]^2);
xxx = NDSolve[{r'[t] == D[H, pr[t]],th'[t] == D[H, pth[t]], pr'[t] == -D[H, r[t]], pth'[t] == -D[H, th[t]],r == 2.29, th == Pi, pr == 1, pth == 0,WhenEvent[r[t]^2 == (b/Sqrt[1 - (eps*Cos[th[t]])^2])^2, pr[t] -> -pr[t]]}, {r[t], th[t], pr[t], pth[t]},{t, 0,100}]
Show[PolarPlot[{xxx[][][]}, {t, 0, 100}],PolarPlot[b/(Sqrt[1 - (eps*Cos[th])^2]), {th, 0, 2*Pi},PlotStyle -> Red]]

The output I get is a very pretty but wrong image: The trajectories must lie inside the ellipse (red), because the reflection condition given by WhenEvent defines the boundary of the circle.

What am I doing wrong?

• I guess this would be simpler in Cartesian coordinates... – Henrik Schumacher Jun 18 '18 at 9:38
• @HenrikSchumacher May be, but the next extension of the problem which I wish to solve will be easier to be solved in Polar Coordinates. I wish to stick to this coordinate system! – Nishchal Dwivedi Jun 19 '18 at 9:23
• Then you have to express the normal of the boundary ellipsoid and apply a reflection of the speed vector with respect to this normal in the WhenEvent. – Henrik Schumacher Jun 19 '18 at 9:29

Here the interpolating functions generated are of $r(t)$ and $\theta(t)$. We need to do a ParametricPlot of $r(t)cos(\theta(t))$ vs $r(t)sin(\theta(t))$ and it will give the correct graph. • Yes, PolarPlot specifically is looking to plot the radius as a function of $\theta$, not a parametric curve. – KraZug Jun 19 '18 at 10:49