How to simplify the calculations of reflecting a ray about an ellipse

Question

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.

Answer 1

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.

Answer 2

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[0]} [ -h | --help ]\n")
        stderr.write(f"       {argv[0]} 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[1])
    if (rx is None) or (rx <= 0):
        stderr.write(f"{argv[1]}: Invalid semiaxis length (RX).\n")
        exit(1)

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

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

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

    if x0*x0*ry*ry + y0*y0*rx*rx >= rx*rx*ry*ry:
        stderr.write("Ray does not start within the ellipse.\n")
        exit(1)

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

    try:
        count = int(argv[6])
    except ValueError:
        count = 0
    if count < 1:
        stderr.write(f"{argv[6]}: 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.

Answer 3

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

Answer 4

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.

Answer 5

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

Answer 6

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

3. 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} $$

4. 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$$

5. 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$$

6. 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}}$$

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