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.
 A: 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$.
A: 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.
A: 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

*

*inputs:  $\phi_0$, $\phi_1$


*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$$


*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}
$$


*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$$


*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$$


*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}}$$


*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)
A: 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}$.

*

*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.


*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}$.
A: 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.

