# How to simplify the calculations of reflecting a ray about an ellipse    I wrote the script that made these images several days ago, the segments each depict a ray of light, as the light hits the boundary of the ellipse, it is reflected by the ellipse according to the laws of reflection, and the reflected ray of light is again reflected by the ellipse, and the reflection of the reflection is again reflected by the ellipse...

The light keeps bouncing back and forth, over and over again, until the light has been reflected certain number of times.

Step by step on how I made these images.

First, you make an ellipse.

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$$

The above equation describes the points on the ellipse, it is also the boundary of the ellipse. (Assume the ellipse is horizontal).

Construct a right triangle from one semi-minor axis and one semi-major axis, let the acute angle adjacent to the semi-major axis be $$\alpha$$, then the relationship between a, b and $$\alpha$$ can be written as:

$$b = a \cdot tan(\alpha)$$

Rewrite the equation of the ellipse:

$$\frac{x^2}{a^2} + \frac{y^2}{a^2 \cdot tan(\alpha)^2} = 1$$

Pick a random point inside the ellipse, for all points inside by the ellipse, simply use this equation:

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} <= 1$$

Then the random point is given as:

$$(a \cdot cos(\alpha) \cdot m, b \cdot sin(\alpha) \cdot n)$$

Where $$\alpha$$ is in range $$[0, 2\pi]$$, and m, n are in range $$[0, 1]$$.

Then pick a random angle $$\beta$$, to construct a line passes the chosen point at angle $$\beta$$ with the x-axis.

I use the slope-intercept form of line equation ($$y = k \cdot x + c$$):

Let the chosen point be $$(x_0, y_0)$$, then the equation of the incident ray is:

$$y = tan(\beta) \cdot x + y_0 - x_0 \cdot tan(\beta)$$

But if $$\beta$$ is a multiple of $$\frac{\pi} {2}$$ things get complicated, because tan(0) = 0 and $$tan(\frac{\pi} {2})$$ is undefined.

So I use $$y = y_0$$ if the line is parallel to x-axis and $$x = x_0$$ if the line is perpendicular to x-axis.

Now to calculate the intersections between the incident ray and ellipse.

If:

$$c^2 < a^2 \cdot k^2 + b^2$$

Then there can be two intersections $$(x_0, y_0)$$ and $$(x_1, y_1)$$.




\begin{align} n_1 &= a^2 \cdot k^2 + b^2 \\ n_2 &= 2 \cdot a^2 \cdot k \cdot c \\ n_3 &= a^2 \cdot (c^2 - b^2) \\ n_4 &= \sqrt{(n_2^2 - 4 \cdot n_1 \cdot n_3)} \\ x_0 &= \frac{(-n_2 + n_4)} {(2 \cdot n_1)} \\ x_1 &= \frac{(-n_2 - n_4)} {(2 \cdot n_1)} \\ y_0 &= k \cdot x_0 + c \\ y_1 &= k \cdot x_1 + c \end{align}

For lines like $$x = x_0$$ and $$y = y_0$$ however:

First set:

If $$abs(x_0) <= a$$

$$y_0 = \sqrt{b^2 - x_0^2 \cdot \frac{b^2} {a^2}}$$

Intersections are $$(x_0, +y_0)$$ and $$(x_0, -y_0)$$.

Second set:

If $$abs(y_0) <= b$$

$$x_0 = \sqrt{a^2 - y_0^2 \cdot \frac{a^2} {b^2}}$$

Intersections are $$(+x_0, y_0)$$ and $$(-x_0, y_0)$$

Then choose the intersection to do further calculations, for $$\beta$$ in range $$[0, \frac{\pi} {2}]$$ and $$[\frac{3 \pi} {2}, 2 \pi]$$ I choose the right intersection, else I choose the left one.

Then I calculate the tangent of the ellipse at the intersection.

For a point $$(x_0, y_0)$$ on the ellipse given by a, b, the slope k of the line tangent to the ellipse at that point must satisfy:

$$k = \frac{-x_0 \cdot b^2} {a^2 \cdot y_0}$$

Then for almost all lines, the equation of the tangential line is:

$$y = k \cdot x + \frac{b^2} {y_0}$$

But if arctan(k) is a multiple of $$\frac{\pi} {2}$$ the above relationship breaks down and I instead fall back to constant form.

Then I calculate the normal (line perpendicular to that tangent passing through that intersection):

$$y = \frac{-x} {k} + y_0 - \frac{-x_0} {k}$$

Where $$(x_0, y_0)$$ is the intersection, and k is the slope of the normal.

But again the above relationship breaks down if the line is special. I won't show how I deal with exceptions here, I already have shown too many equations, you can see all the calculations in the code.

Then I calculate the signed angle formed by the incident ray and the normal:

Let $$k_1$$ be the slope of the incident ray, let $$k_2$$ be the slope of the normal:




\begin{align} \alpha_1 &= atan(k_1) \\ \alpha_2 &= atan(k_2) \\ \alpha_\delta &= \alpha_2 - \alpha_1 \\ \alpha_\delta &= (\alpha_\delta + \pi) \bmod 2 \pi - \pi \end{align}

Again, the above doesn't work if either of these lines are special, other calculations are required.

Then I calculate the reflected ray, simply by rotate the normal line about the intersection by $$\alpha_\delta$$ (assuming the previous calculations succeeded):

\begin{align} \alpha_3 &= \alpha_2 + \alpha_\delta \\ y &= tan(\alpha_3) \cdot x + y_0 - tan(\alpha_3) \cdot x_0 \end{align}

Then I calculated the intersections between the reflected ray and the ellipse, there will be two intersections, this time the intersection needed is the other intersection from the current one.

Then all above calculations are repeated recursively, until a certain number of iteration is reached.

How to simplify all calculations involved, and calculate the reflected ray in as few steps as possible, including all the edge cases?

I want it that way: given the coordinate of a point and an angle, calculate the intersection of the ray with the ellipse, and then calculate the tangent of the ellipse at that intersection, then calculate the reflected ray, all of these in as few steps as possible, using one set of equations without exceptions. Preferably the number of equations involved should be less than or equal to six.

• I think I have to clarify a few things, I am an unemployed individual, and I am a high school drop out, I was born in 1999, so all these questions are neither job-related nor homework-related, I did all of these just for the sake of it, for fun, because I genuinely found these things interesting, but sadly my family doesn't appreciate my work and even ridicule me for doing these things or anything I like that doesn't generate monetary revenue... May 31, 2022 at 8:15
• By the way I did not know about Mathjax before I started writing the question, and I learned Mathjax by clicking the edit button on one of the answers to my previous questions, to see the source code, I learned Mathjax just to ask this question, to show my calculations. It is tedious to edit Mathjax code, so the question contained a wrong equation before I corrected it, because I copied and pasted part of it... May 31, 2022 at 9:29
• I also like how my (but probably everyone's) brain tricks me (us) into seeing the pictures as 3D-shapes. Calculating what are these 3D-shapes exactly would make for an other great (separate) MSE-question. May 31, 2022 at 9:54
• Welcome to Math.SE! <> Not an answer, but you may be interested in mathematical billiards. May 31, 2022 at 16:10
• This may be of interest: math.stackexchange.com/questions/3619214/… May 31, 2022 at 17:44

I think one can simplify the computations using vectors. Let's start with a ray emanating from point $$A$$ on the ellipse and being reflected at another point $$B$$ on the ellipse. The reflected ray is then $$BA'$$, where $$A'$$ is the reflection of point $$A$$ about the normal $$BN$$.

If $$F$$ and $$G$$ are the foci, finding point $$N$$ on segment $$FG$$ is not difficult, because $$BN$$ is the bisector of $$\angle FBG$$, which entails $$FN:GN=FB:GB$$. Hence (remembering that $$FB+GB=2a$$, the major axis of the ellipse): $$N=\left(1-{FB\over2a}\right)F+{FB\over2a}G.$$ We can then find $$M$$, the projection of $$A$$ on line $$BN$$, setting $$M=B+t(B-N)$$ and using that the scalar product $$(B-M)\cdot(A-M)$$ vanishes. From there we can get $$A'=2M-A$$. One finds: $$A'=2B-A+2t(B-N), \quad\text{where:}\quad t=-{(B-A)\cdot(B-N)\over(B-N)\cdot(B-N)}.$$ Once we have $$A'$$ we can find the other intersection $$C$$ of line $$BA'$$ with the ellipse, and repeat the whole process to compute next reflection. • I agree using vectors is indeed the way to go, your equations are indeed much simpler than the other answers, but the reflective properties of ellipses is quite complex, it discriminates among cases where the incident ray is inside the circle defined by the foci, on the circle, passing through the foci, outside the circle and on the ellipse, do your equations work for all cases without failure? May 31, 2022 at 22:02
• @ΞένηΓήινος Point $N$ is always well defined by the above formula, because it works even in the border case when $B$ lies on line $FG$. And $t$ is always well defined too, because $B-N$ cannot vanish ($N$ is inside segment $FG$). So yes, the above construction works in any case. I tested it with GeoGebra and found no problems at all. May 31, 2022 at 22:24
• @Intelligentipauca Can we proof point $N$ is always inside segment $𝐹𝐺$? How about point $N$ is is on long-axis but is also outside segment $𝐹𝐺$? Jun 4, 2022 at 8:57
• @Mountain $BN$ is the bisector of $\angle FBG$, hence $N$ is forced to lie inside $FG$. You can also see that from the equation for $N$ written above. Jun 4, 2022 at 10:41
• Understand! Thanks @Intelligentipauca Jun 4, 2022 at 12:18

A common way to parametrise a ray is \vec{p} = \vec{p}_0 + t \vec{p}_\Delta \quad \iff \quad \left\lbrace ~ \begin{aligned} x &= x_0 + t x_\Delta \\ y &= y_0 + t y_\Delta \\ \end{aligned} \right . \tag{1}\label{1} with $$t \ge 0$$. Without loss of generality, we could choose $$x_\Delta^2 + y_\Delta^2 = 1$$, so that $$t$$ is also the distance along the ray. For example, the direction of the ray can also be described using angle $$\theta$$ (zero at positive $$x$$ axis, increases towards positive $$y$$ axis, and so on): \left\lbrace ~ \begin{aligned} x_\Delta &= \cos \theta \\ y_\Delta &= \sin \theta \\ \end{aligned} \right . \quad \iff \quad \theta = \operatorname{atan2}\left( y_\Delta, x_\Delta \right) \tag{2}\label{2} In this particular case, however, there is no need to do that, unless the length of the ray is somehow interesting. In general, the distance $$d$$ from point $$\vec{p}_0$$ to point $$\vec{p}$$ is $$d = t \sqrt{x_\Delta^2 + y_\Delta^2}$$

All real $$t$$, $$t \in \mathbb{R}$$, corresponds to a line a x + b y = c, \quad \left\lbrace ~ \begin{aligned} a &= - y_\Delta \\ b &= x_\Delta \\ c &= y_0 x_\Delta - x_0 y_\Delta \\ \end{aligned} \right . \tag{3}\label{3} For non-vertical lines ($$x_\Delta = b \ne 0$$), this corresponds to y = m x + d, \quad \left\lbrace ~ \begin{aligned} m &= \displaystyle \frac{y_\Delta}{x_\Delta} \\ d &= \displaystyle y_0 - \frac{y_\Delta}{x_\Delta} x_0 \\ \end{aligned} \right . \tag{4}\label{4} In a very real sense, equation $$\eqref{2}$$ is the most general one, but it describes a line and not a ray. The form in equation $$\eqref{1}$$ is most commonly used in raycasting (both 2D and 3D), as it is easy to convert to the other forms if needed, and intersections are easy to find by substituting equation $$\eqref{1}$$ to the implicit equation of the curve/surface and solving for $$t$$. If there are multiple solutions, the smallest positive $$t$$ is the first intersection in the direction of the ray.

If we use $$r_x$$ and $$r_y$$ for the semiaxes of the ellipse, $$\frac{x^2}{r_x^2} + \frac{y^2}{r_y^2} = 1 \tag{5}\label{5}$$ substituting $$\eqref{1}$$ into $$\eqref{5}$$ yields $$\frac{(x_0 + t x_\Delta)^2}{r_x^2} + \frac{(y_0 + t y_\Delta)^2}{r_y^2} = 1$$ which is obviously a quadratic equation in $$t$$. It simplifies to $$T_2 t^2 + 2 T_1 t + T = 0$$ where \begin{aligned} T_2 &= r_x^2 y_\Delta^2 + r_y^2 x_\Delta^2 \\ T_1 &= r_x^2 y_\Delta y_0 + r_y^2 x_\Delta x_0 \\ T_0 &= r_x^2 y_0^2 + r_y^2 x_0^2 - r_x^2 r_y^2 \\ \end{aligned} \tag{6a}\label{6a} and has solutions $$t = -\frac{T_1}{T_2} \pm \frac{\sqrt{T_1^2 - T_0 T_2}}{T_2} \tag{6b}\label{6b}$$ Numerically, assuming the initial starting point is within the ellipse, there should be two solutions for each following intersection, the larger corresponding to the next intersection, and the smaller one the previous intersection (and should be exactly zero, except that when using floating-point values, some rounding error can creep in, when the starting point for the current ray is on the ellipse).

When we use $$\eqref{1}$$ to describe the ray, reflection from a surface with normal $$\vec{n}$$, ($$\vec{n} = (x_n, y_n)$$ perpendicular to tangent(s), in the half-space towards the ray starting point) is \vec{p}_\Delta^\prime = \vec{p}_\Delta - 2 \frac{\vec{p}_\Delta \cdot \vec{n}}{\vec{n} \cdot \vec{n}} \vec{n} \quad \iff \quad \left\lbrace \begin{aligned} x_\Delta^\prime &= x_\Delta - 2 \frac{x_\Delta x_n + y_\Delta y_n}{x_n^2 + y_n^2} x_n \\ y_\Delta^\prime &= y_\Delta - 2 \frac{x_\Delta x_n + y_\Delta y_n}{x_n^2 + y_n^2} y_n \\ \end{aligned} \right . \tag{7}\label{7}

In parametrised form, we can describe the same ellipse as $$\eqref{5}$$ as \left\lbrace ~ \begin{aligned} x &= r_x \cos \varphi \\ y &= r_y \sin \varphi \\ \end{aligned} \right . where $$\varphi$$ is the angular parameter. In 2D, the tangent of any such parametrised curve is $$\left(\frac{d x}{d \varphi}, \frac{d y}{d \varphi}\right)$$. If the curve is counterclockwise like the above one is, then the inwards normal is $$\vec{n} = (x_n, y_n) = \left(-\frac{d y}{d \varphi}, \frac{d x}{d \varphi}\right)$$: \left\lbrace ~ \begin{aligned} x_n &= - r_y \cos\varphi = -\frac{r_y}{r_x} x \\ y_n &= - r_x \sin\varphi = -\frac{r_x}{r_y} y \\ \end{aligned} \right. \tag{8}\label{8} Note that the right side is only valid if $$(x, y)$$ is on the ellipse. Substituting $$\eqref{8}$$ into $$\eqref{7}$$ we get the reflection at point $$(x, y)$$ on the ellipse for a ray with direction $$(x_\Delta, y_\Delta)$$: \left\lbrace ~ \begin{aligned} x_\Delta^\prime &= x_\Delta - 2 \frac{ x_\Delta x^2 r_y^4 + y_\Delta x y r_x^2 r_y^2 }{x^2 r_y^4 + y^2 r_x^4} \\ y_\Delta^\prime &= y_\Delta - 2 \frac{ y_\Delta y^2 r_x^4 + x_\Delta x y r_x^2 r_y^2 }{x^2 r_y^4 + y^2 r_x^4} \\ \end{aligned} \right . \tag{9}\label{9}

To recap, we have a ray starting at $$(x_0, y_0)$$ with direction $$(x_\Delta, y_\Delta)$$, initially $$(\cos\theta, \sin\theta)$$ where $$\theta$$ is the direction angle. We assume $$(x_0, y_0)$$ is within the ellipse.

We use $$\eqref{6a}$$ and $$\eqref{6b}$$ to find the next intersection with the ellipse, by choosing the larger of the two solutions $$t$$.

The intersection on the ellipse is then $$x = x_0 + t x_\Delta$$, $$y = y_0 + t y_\Delta$$. We use $$\eqref{9}$$ to calculate the new direction $$(x_\Delta^\prime, y_\Delta^\prime)$$.

For the next ray, we set $$x_0 \gets x$$, $$y_0 \gets y$$, $$x_\Delta \gets x_\Delta^\prime$$, $$y_\Delta \gets y_\Delta^\prime$$, and repeat.

Here is an example Python 3 implementation:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from math import sqrt, sin, cos, radians, isfinite

def toFloat(value):
try:
f = float(value)
if not isfinite(f):
f = None
except ValueError:
f = None
return f

def intersections(rx2, ry2, x0, y0, xd, yd):
t2 = rx2*yd*yd + ry2*xd*xd
t1 = rx2*y0*yd + ry2*x0*xd
t0 = rx2*y0*y0 + ry2*x0*x0 - rx2*ry2
try:
ca = t1 / t2
cb = sqrt(t1*t1 - t0*t2) / t2
return cb - ca, -cb - ca
except ValueError:
return None

def ellipsereflect(rx2, ry2, x, y, xd, yd):
ca = x*x*ry2*ry2
cb = y*y*rx2*rx2
cc = x*y*rx2*ry2
cd = (ca + cb) / 2
try:
return xd - (xd * ca + yd * cc) / cd, yd - (xd * cc + yd * cb) / cd
except ValueError:
return None, None

if __name__ == '__main__':
from sys import argv, stdout, stderr, exit

if len(argv) != 7:
stderr.write("\n")
stderr.write(f"Usage: {argv} [ -h | --help ]\n")
stderr.write(f"       {argv} RX RY X0 Y0 ANGLE N > IMAGE.svg\n")
stderr.write("\n")
stderr.write("Where RX and RY are the ellipse semiaxis lengths,\n")
stderr.write("the first ray starts at (X0, Y0) in direction ANGLE,\n")
stderr.write("with ANGLE in degrees (0 to 360), and we trace N\n")
stderr.write("reflections.  The output is an SVG image.\n")
stderr.write("\n")
exit(1)

rx = toFloat(argv)
if (rx is None) or (rx <= 0):
stderr.write(f"{argv}: Invalid semiaxis length (RX).\n")
exit(1)

ry = toFloat(argv)
if (ry is None) or (ry <= 0):
stderr.write(f"{argv}: Invalid semiaxis length (RY).\n")
exit(1)

x0 = toFloat(argv)
if x0 is None:
stderr.write(f"{argv}: Invalid ray starting point X coordinate (X0).\n")
exit(1)

y0 = toFloat(argv)
if y0 is None:
stderr.write(f"{argv}: Invalid ray starting point Y coordinate (Y0).\n")
exit(1)

if x0*x0*ry*ry + y0*y0*rx*rx >= rx*rx*ry*ry:
exit(1)

theta = toFloat(argv)
if theta is None:
stderr.write(f"{argv}: Invalid ray starting angle in degrees (ANGLE).\n")
exit(1)
else:
xd = cos(theta)
yd = sin(theta)

try:
count = int(argv)
except ValueError:
count = 0
if count < 1:
stderr.write(f"{argv}: Invalid number of reflections to calculate (N).\n")
exit(1)

stdout.write('<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n')
stdout.write('<svg xmlns="http://www.w3.org/2000/svg" version="1.1" viewBox="%.0f %.0f %.0f %.0f">\n' % (round(-rx-4), round(-ry-4), round(2*rx+8), round(2*ry+8)))
# stdout.write('<rect x="%.0f" y="%.0f" width="%.0f" height="%.0f" fill="#ffffff" stroke="none" />\n' % (round(-rx-4), round(-ry-4), round(2*rx+8), round(2*ry+8)))
stdout.write('<ellipse cx="0" cy="0" rx="%.3f" ry="%.3f" fill="#ffffff" stroke="#000000" />\n' % (rx, ry))
stdout.write('<path d="M%.3f%+.3f' % (x0, y0))

x, y = x0, y0
for i in range(0, count):
t1, t2 = intersections(rx*rx, ry*ry, x, y, xd, yd)
t = max(t1, t2)

x += t * xd
y += t * yd
stdout.write(' L%.3f%+.3f' % (x, y))
xd, yd = ellipsereflect(rx*rx, ry*ry, x, y, xd, yd)

stdout.write('" fill="none" stroke="#0000ff" />\n')
stdout.write('</svg>\n')


For an example, run python3 example.py 500 200 0 50 0 158 > example.svg or python3 example.py 500 200 0 0 45 170 > example.svg, and view the generated example.svg in a browser for example. This draws an ellipse with semiaxes $$500$$ and $$200$$, with the initial ray starting at $$0,50$$ or $$0,0$$ at angle $$0$$ (along positive $$x$$ axis) or $$45°$$, and the reflections up to the $$158$$th or $$170$$th intersection with the ellipse.

Here's my go at it. First of all since we are only interested in drawing the result we can take $$b=1$$. At each iteration specify ray which consists of a point $$(x_0,y_0)$$ as well as a normal vector $$(n_x,n_y)$$ which determines the direction of the ray. An iteration consists of intersecting the ray with the ellipse and then reflecting the normal vector in the new point.

We can determine the intersection by dividing $$x$$ and $$n_x$$ by $$a$$, which makes the ellipse a circle, then intersecting the ray with the circle and then scaling $$x$$ and $$n_x$$ with $$a$$ again. The circle intersection can be found using $$\begin{equation} \begin{cases} x^2+y^2=1\\ y-y_0=\frac{n_y}{n_x}(x-x_0) \end{cases} \end{equation}$$ which results in \begin{align} x=\frac{ x_0(n_y^2-n_x^2)-2 {n_x} {n_y} y_0}{{n_x}^2+{n_y}^2}\\ y=\frac{ y_0(n_x^2-n_y^2)-2 {n_x} {n_y} x_0}{{n_x}^2+{n_y}^2} \end{align}

Now we need to perform the reflection of the normal. First we need a normal to the ellipse using the tangent. We can find the tangent by first parametrizing $$(x_0,y_0)=(a\cos\theta,\sin\theta)$$ and then taking the $$\theta$$-derivative. This gives a tangent $$(-a\sin\theta,\cos\theta)$$ from which we can construct the (inward facing) normal $$(-\cos\theta,-a\sin\theta)$$. We can remove the need for trigonometric functions by writing it again as $$(-x_0/a,-a y_0)$$. Call this normal $$\vec N=(N_x,N_y)$$. We can then decompose our normal vector in components parallel and perpendicular to $$\vec N$$, i.e. $$\vec n=\vec n_\parallel+n_\perp$$. The reflected vector is then \begin{align} \vec n'&=\vec n_\perp-\vec n_\parallel\\ &=\vec n-2\vec n_\parallel\\ &=\vec n-\frac{\vec N\cdot\vec n}{\lvert{N}\rvert^2}\vec N \end{align}

Now we can iterate this as long as we want. It might be beneficial to stay in the "circle coordinates", i.e. once we divide the $$x$$-axis by $$a$$ we do not scale up again until we want to draw the points. We would then also have to scale $$\vec N$$.

I notice you are using $$y=kx+c$$ as the general equation of a line. As you realise, this breaks down for vertical lines, where $$k$$ would have to be infinite.

The most general equation for a line is $$rx+sy=t$$. This has the minor disadvantage that there are three coefficients rather than two, and that the representation is not unique - $$r,s,t$$ can be scaled by the same factor without changing the line. If the latter bothers you, you can always scale things to make $$t=1$$ or $$t=0$$ to make the representation unique. $$t$$ then becomes more or less a flag to indicate whether the line goes through the origin or not.

Anyway, I would recommend using this more general equation $$rx+sy=t$$ with three coefficients $$(r,s,t)$$ for any line in the plane without needing special cases.

• You haven't shown the actual calculations, I had done all those calculations by myself, and demonstrated all of them, if you want your answer to be accepted, I would expect you to also show the calculations, so please do some actual calculations. May 31, 2022 at 9:50
• Otherwise I recommend you to just delete your answer. It doesn't really answer my question, after all. May 31, 2022 at 9:54
• @ΞένηΓήινος I know, but this was too long to be just a comment on your post. Still, I think you'll find things simplify considerably if you consistently use the more general equation so that you can avoid denominators that may or may not be zero. May 31, 2022 at 9:54

I assume that the ellipse is given by the equation

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1.$$

We are given a starting point $$\begin{pmatrix} x_0 \\ y_0 \end{pmatrix}$$ inside the ellipse and the ray direction vector $$\begin{pmatrix} u \\ v \end{pmatrix}$$.

1. The point of the intersection is

$$\tag{1} \begin{pmatrix} x_1 \\ y_1 \end{pmatrix} = \begin{pmatrix} x_0 \\ y_0 \end{pmatrix} + t \cdot \begin{pmatrix} u \\ v \end{pmatrix},$$

where $$t \in \mathbb{R}$$, satisfying

$$\frac{(x_0+tu)^2}{a^2} + \frac{(y_0+tv)^2}{b^2} = 1.$$

Short calculations yield

$$\tag{2} t = \frac{\sqrt{\frac{u^2}{a^2} + \frac{v^2}{b^2} - \left( \frac{uy_0 - vx_0}{ab} \right)^2} - \left( \frac{ux_0}{a^2} + \frac{vy_0}{b^2} \right)}{\frac{u^2}{a^2} + \frac{v^2}{b^2}},$$

which is the positive solution of the quadratic equation because the ray only moves forward.

2. The normal vector $$N$$ to the ellipse at the point $$\begin{pmatrix} x_1 \\ y_1 \end{pmatrix}$$ is the gradient of the function $$F \begin{pmatrix} x \\ y \end{pmatrix} = \frac{x^2}{a^2} + \frac{y^2}{b^2}$$, i.e.

$$\tag{3} N = \nabla F \begin{pmatrix} x_1 \\ y_1 \end{pmatrix} = \begin{pmatrix} \frac{2x_1}{a^2} \\ \frac{2y_1}{b^2} \end{pmatrix}.$$

The new ray direction $$\begin{pmatrix} u' \\ v' \end{pmatrix}$$ is given by the formula

$$\tag{4} \begin{pmatrix} u' \\ v' \end{pmatrix} = \begin{pmatrix} u' \\ v' \end{pmatrix} - 2 P_N \begin{pmatrix} u' \\ v' \end{pmatrix} = \begin{pmatrix} u' \\ v' \end{pmatrix} - 2 \cdot \frac{\left< \begin{pmatrix} u' \\ v' \end{pmatrix}, N \right>}{\left< N, N \right>} \cdot N.$$

So ultimately you need four equations, $$(1)$$ - $$(4)$$, to compute the new starting point $$\begin{pmatrix} x_1 \\ y_1 \end{pmatrix}$$ and the new ray direction $$\begin{pmatrix} u' \\ v' \end{pmatrix}$$.

Suppose the ellipse is decided by

$$\begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}a \cos \phi\\b \sin \phi\end{pmatrix}$$

And $$a >= b$$, $$a^2 + b^2 = c^2$$ Given points $$\phi_0$$, $$\phi_1$$, $$\phi_2$$ decided by above equation and the reflection relationship, we will list the steps of the algorithm

1. inputs: $$\phi_0$$, $$\phi_1$$

2. calculate $$c_0, s_0, c_1, s_1$$

$$c_0 = \cos \phi_0, s_0 = \sin \phi_0$$

$$c_1 = \cos \phi_1, s_1 = \sin \phi_1$$

1. calculate $$t_1$$

$$t_1 = \frac{ab(1 - c_0c_1 - s_0s_1)}{a^2c_0s_1 - b^2s_0c_1 - c^2 c_1s_1}$$

1. calculate $$u_1, v_1, w_1$$

$$u_1 = abc_1 - a^2t_1s_1$$

$$v_1 = abs_1 + b^2t_1c_1$$

$$w_1 = ab - b^2t_1c_1$$

1. calculate $$c_2, s_2$$

$$c_2 = \frac{2u_1w_1}{u_1^2 + v_1^2} - c_1$$

$$s_2 = \frac{2v_1w_1}{u_1^2 + v_1^2} - s_1$$

1. correction of the error

$$c_2 = \frac{c_2}{\sqrt{c_2^2 + s_2^2}}$$

$$s_2 = \frac{s_2}{\sqrt{c_2^2 + s_2^2}}$$

1. calculate $$\phi_2$$ from $$c_2$$ and $$s_2$$, we can use arctan2 from many numerical libs

The proof of above algorithm is elementary.

The most interesting part is that

(1) If the ray pass the one foci, then it will pass another foci, and then is trapped inside the ellipse

(2) If the ray pass beside the two foci, it will always beside the two foci

(3) If the ray intersect the interval decided by the two foci, then it will alway intersect the interval

Ellipse can trap the ray makes me astonished at the time.

I also want to know whether some kind of ergodicity holds for (2) and (3)

• Unfortunately the calculations are not much simpler than my own. May 31, 2022 at 11:42
• Yes, thanks for your contribution May 31, 2022 at 11:43