Cómo simplificar los cálculos de la reflexión de un rayo sobre una elipse

Question

He escrito el guión que hizo estas imágenes hace varios días, los segmentos representan cada uno un rayo de luz, como la luz golpea el límite de la elipse, es reflejada por la elipse de acuerdo con las leyes de la reflexión, y el rayo de luz reflejado es de nuevo reflejado por la elipse, y la reflexión de la reflexión es de nuevo reflejada por la elipse...

La luz sigue rebotando de un lado a otro, una y otra vez, hasta que la luz se ha reflejado un determinado número de veces.

Paso a paso sobre cómo hice estas imágenes.

Primero, se hace una elipse.

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

La ecuación anterior describe los puntos de la elipse, también es el límite de la elipse. (Supongamos que la elipse es horizontal).

Construir un triángulo rectángulo a partir de un semieje menor y un semieje mayor, siendo el ángulo agudo adyacente al semieje mayor $\alpha$ entonces la relación entre a, b y $\alpha$ se puede escribir como:

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

Reescribe la ecuación de la elipse:

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

Elige un punto al azar dentro de la elipse, para todos los puntos dentro de la elipse, simplemente usa esta ecuación:

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

Entonces el punto aleatorio está dado como

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

Donde $\alpha$ está en el rango $[0, 2\pi]$ y m, n están en el rango $[0, 1]$ .

Entonces elige un ángulo al azar $\beta$ para construir una línea que pase por el punto elegido con un ángulo $\beta$ con el eje x.

Utilizo la forma de intersección de pendientes de la ecuación de la línea ( $y = k \cdot x + c$ ):

Sea el punto elegido $(x_0, y_0)$ , entonces la ecuación del rayo incidente es:

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

Pero si $\beta$ es un múltiplo de $\frac{\pi} {2}$ las cosas se complican, porque tan(0) = 0 y $tan(\frac{\pi} {2})$ es indefinido.

Así que uso $y = y_0$ si la línea es paralela al eje x y $x = x_0$ si la línea es perpendicular al eje x.

Ahora hay que calcular las intersecciones entre el rayo incidente y la elipse.

Si:

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

Entonces puede haber dos intersecciones $(x_0, y_0)$ y $(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}$$

Para líneas como $x = x_0$ y $y = y_0$ sin embargo:

Primer set:

Si $abs(x_0) <= a$

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

Las intersecciones son $(x_0, +y_0)$ y $(x_0, -y_0)$ .

Segundo set:

Si $abs(y_0) <= b$

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

Las intersecciones son $(+x_0, y_0)$ y $(-x_0, y_0)$

A continuación, elija la intersección para realizar más cálculos, para $\beta$ en el rango $[0, \frac{\pi} {2}]$ y $[\frac{3 \pi} {2}, 2 \pi]$ Elijo la intersección de la derecha, si no, elijo la de la izquierda.

Luego calculo la tangente de la elipse en la intersección.

Para un punto $(x_0, y_0)$ en la elipse dada por a, b, la pendiente k de la recta tangente a la elipse en ese punto debe satisfacer:

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

Entonces para casi todas las líneas, la ecuación de la línea tangencial es:

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

Pero si arctan(k) es un múltiplo de $\frac{\pi} {2}$ la relación anterior se rompe y en su lugar vuelvo a la forma constante.

Luego calculo la normal (línea perpendicular a esa tangente que pasa por esa intersección):

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

Donde $(x_0, y_0)$ es la intersección, y k es la pendiente de la normal.

Pero de nuevo la relación anterior se rompe si la línea es especial. No mostraré aquí cómo trato las excepciones, ya he mostrado demasiadas ecuaciones, puedes ver todos los cálculos en el código.

Luego calculo el ángulo con signo formado por el rayo incidente y la normal:

Dejemos que $k_1$ sea la pendiente del rayo incidente, y que $k_2$ sea la pendiente de la 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}$$

De nuevo, lo anterior no funciona si alguna de estas líneas es especial, se requieren otros cálculos.

Luego calculo el rayo reflejado, simplemente girando la línea normal sobre la intersección por $\alpha_\delta$ (suponiendo que los cálculos anteriores hayan tenido éxito):

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

Luego calculé las intersecciones entre el rayo reflejado y la elipse, habrá dos intersecciones, esta vez la intersección que se necesita es la otra intersección desde la actual.

A continuación, todos los cálculos anteriores se repiten de forma recursiva, hasta alcanzar un determinado número de iteraciones.

¿Cómo simplificar todos los cálculos implicados y calcular el rayo reflejado en el menor número de pasos posible, incluyendo todos los casos límite?

---

Lo quiero así: dada la coordenada de un punto y un ángulo, calcular la intersección del rayo con la elipse, y luego calcular la tangente de la elipse en esa intersección, luego calcular el rayo reflejado, todo esto en el menor número de pasos posible, utilizando un conjunto de ecuaciones sin excepciones. Preferiblemente, el número de ecuaciones involucradas debe ser menor o igual a seis.

Answer 1

Creo que se pueden simplificar los cálculos utilizando vectores. Empecemos con un rayo que emana del punto $A$ en la elipse y que se refleja en otro punto $B$ en la elipse. El rayo reflejado es entonces $BA'$ , donde $A'$ es la reflexión del punto $A$ sobre lo normal $BN$ .

Si $F$ y $G$ son los focos, encontrando el punto $N$ en el segmento $FG$ no es difícil, porque $BN$ es la bisectriz de $\angle FBG$ lo que implica $FN:GN=FB:GB$ . Por lo tanto (recordando que $FB+GB=2a$ el eje mayor de la elipse): $$ N=\left(1-{FB\over2a}\right)F+{FB\over2a}G. $$ Entonces podemos encontrar $M$ la proyección de $A$ en línea $BN$ , ajuste $M=B+t(B-N)$ y utilizando que el producto escalar $(B-M)\cdot(A-M)$ desaparece. De ahí podemos obtener $A'=2M-A$ . Uno encuentra: $$ A'=2B-A+2t(B-N), \quad\text{where:}\quad t=-{(B-A)\cdot(B-N)\over(B-N)\cdot(B-N)}. $$ Una vez que tenemos $A'$ podemos encontrar la otra intersección $C$ de la línea $BA'$ con la elipse, y repetir todo el proceso para calcular la siguiente reflexión.

Answer 2

Una forma común de parametrizar un rayo es $$\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}$$ con $t \ge 0$ . Sin pérdida de generalidad, podríamos elegir $x_\Delta^2 + y_\Delta^2 = 1$ para que $t$ es también la distancia a lo largo del rayo. Por ejemplo, la dirección del rayo también puede describirse mediante el ángulo $\theta$ (cero en positivo $x$ eje, aumenta hacia el positivo $y$ eje, etc.): $$\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}$$ Sin embargo, en este caso concreto, no es necesario hacerlo, a no ser que la longitud del rayo sea de algún modo interesante. En general, la distancia $d$ desde el punto $\vec{p}_0$ para señalar $\vec{p}$ es $$d = t \sqrt{x_\Delta^2 + y_\Delta^2}$$

Todo lo real $t$ , $t \in \mathbb{R}$ corresponde a una línea $$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}$$ Para las líneas no verticales ( $x_\Delta = b \ne 0$ ), esto corresponde a $$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}$$ En un sentido muy real, la ecuación $\eqref{2}$ es el más general, pero describe una línea y no un rayo. La forma en la ecuación $\eqref{1}$ es la más utilizada en el raycasting (tanto en 2D como en 3D), ya que es fácil de convertir a las otras formas si es necesario, y las intersecciones son fáciles de encontrar sustituyendo la ecuación $\eqref{1}$ a la ecuación implícita de la curva/superficie y resolviendo para $t$ . Si hay varias soluciones, la más pequeña positiva $t$ es la primera intersección en la dirección del rayo.

Si utilizamos $r_x$ y $r_y$ para los semiejes de la elipse, $$\frac{x^2}{r_x^2} + \frac{y^2}{r_y^2} = 1 \tag{5}\label{5}$$ sustituyendo $\eqref{1}$ en $\eqref{5}$ rinde $$\frac{(x_0 + t x_\Delta)^2}{r_x^2} + \frac{(y_0 + t y_\Delta)^2}{r_y^2} = 1$$ que obviamente es una ecuación cuadrática en $t$ . Se simplifica a $$T_2 t^2 + 2 T_1 t + T = 0$$ donde $$\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}$$ y tiene soluciones $$t = -\frac{T_1}{T_2} \pm \frac{\sqrt{T_1^2 - T_0 T_2}}{T_2} \tag{6b}\label{6b}$$ Numéricamente, asumiendo que el punto inicial de partida está dentro de la elipse, debería haber dos soluciones para cada intersección siguiente, la más grande correspondiente a la siguiente intersección, y la más pequeña a la intersección anterior (y debería ser exactamente cero, excepto que cuando se usan valores de punto flotante, puede aparecer algún error de redondeo, cuando el punto de partida de la semirrecta actual está en la elipse).

Cuando usamos $\eqref{1}$ para describir el rayo, reflejo de una superficie con normalidad $\vec{n}$ , ( $\vec{n} = (x_n, y_n)$ perpendicular a la(s) tangente(s), en el semiespacio hacia el punto de partida del rayo) es $$\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}$$

En forma parametrizada, podemos describir la misma elipse como $\eqref{5}$ como $$\left\lbrace ~ \begin{aligned} x &= r_x \cos \varphi \\ y &= r_y \sin \varphi \\ \end{aligned} \right .$$ donde $\varphi$ es el parámetro angular. En 2D, la tangente de cualquier curva parametrizada es $\left(\frac{d x}{d \varphi}, \frac{d y}{d \varphi}\right)$ . Si la curva es en sentido contrario a las agujas del reloj como la anterior, entonces la normal hacia adentro es $\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}$$ Tenga en cuenta que el lado derecho sólo es válido si $(x, y)$ está en la elipse. Sustituyendo $\eqref{8}$ en $\eqref{7}$ obtenemos la reflexión en el punto $(x, y)$ en la elipse para un rayo con dirección $(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}$$

---

Para recapitular, tenemos un rayo que comienza en $(x_0, y_0)$ con dirección $(x_\Delta, y_\Delta)$ , inicialmente $(\cos\theta, \sin\theta)$ donde $\theta$ es el ángulo de dirección. Suponemos que $(x_0, y_0)$ está dentro de la elipse.

Utilizamos $\eqref{6a}$ y $\eqref{6b}$ para encontrar la siguiente intersección con la elipse, eligiendo la mayor de las dos soluciones $t$ .

La intersección en la elipse es entonces $x = x_0 + t x_\Delta$ , $y = y_0 + t y_\Delta$ . W $\eqref{9}$ para calcular la nueva dirección $(x_\Delta^\prime, y_\Delta^\prime)$ .

Para el siguiente rayo, ponemos $x_0 \gets x$ , $y_0 \gets y$ , $x_\Delta \gets x_\Delta^\prime$ , $y_\Delta \gets y_\Delta^\prime$ y repite.

Aquí hay un ejemplo de implementación de Python 3:

#!/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')

Para un ejemplo, ejecute python3 example.py 500 200 0 50 0 158 > example.svg ou python3 example.py 500 200 0 0 45 170 > example.svg y ver los datos generados example.svg en un navegador, por ejemplo. Esto dibuja una elipse con semiejes $500$ y $200$ con el rayo inicial que comienza en $0,50$ ou $0,0$ en ángulo $0$ (a lo largo de la positiva $x$ eje) o $45°$ y los reflejos hasta el $158$ o $170$ de la intersección con la elipse.

Answer 3

Aquí está mi intento. En primer lugar, como sólo nos interesa dibujar el resultado, podemos tomar $b=1$ . En cada iteración especificar el rayo que consiste en un punto $(x_0,y_0)$ así como un vector normal $(n_x,n_y)$ que determina la dirección del rayo. Una iteración consiste en intersecar el rayo con la elipse y luego reflejar el vector normal en el nuevo punto.

Podemos determinar la intersección dividiendo $x$ y $n_x$ por $a$ lo que convierte la elipse en un círculo, luego se interseca el rayo con el círculo y luego se escala $x$ y $n_x$ con $a$ de nuevo. La intersección del círculo se puede encontrar utilizando \begin{equation} \begin{cases} x^2+y^2=1\\ y-y_0=\frac{n_y}{n_x}(x-x_0) \end{cases} \end{equation} lo que resulta en \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}

Ahora tenemos que realizar la reflexión de la normal. Primero necesitamos una normal a la elipse utilizando la tangente. Podemos encontrar la tangente parametrizando primero $(x_0,y_0)=(a\cos\theta,\sin\theta)$ y luego tomar el $\theta$ -derivada. Esto da una tangente $(-a\sin\theta,\cos\theta)$ a partir de la cual podemos construir la normal (orientada hacia el interior) $(-\cos\theta,-a\sin\theta)$ . Podemos eliminar la necesidad de las funciones trigonométricas escribiéndolo de nuevo como $(-x_0/a,-a y_0)$ . Llama a esto normal $\vec N=(N_x,N_y)$ . Podemos entonces descomponer nuestro vector normal en componentes paralelas y perpendiculares a $\vec N$ es decir $\vec n=\vec n_\parallel+n_\perp$ . El vector reflejado es entonces \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}

Ahora podemos iterar esto todo lo que queramos. Podría ser beneficioso permanecer en las "coordenadas del círculo", es decir, una vez que dividimos el $x$ -eje por $a$ no volvemos a escalar hasta que queramos dibujar los puntos. Entonces también tendríamos que escalar $\vec N$ .

Answer 4

Me he dado cuenta de que está utilizando $y=kx+c$ como la ecuación general de una recta. Como se puede comprobar, esto se rompe en el caso de las líneas verticales, en las que $k$ tendría que ser infinito.

La ecuación más general de una recta es $rx+sy=t$ . Esto tiene la pequeña desventaja de que hay tres coeficientes en lugar de dos, y que la representación no es única - $r,s,t$ pueden ser escalados por el mismo factor sin cambiar la línea. Si esto último le molesta, siempre puede escalar las cosas para hacer $t=1$ ou $t=0$ para que la representación sea única. $t$ se convierte entonces más o menos en una bandera para indicar si la línea pasa por el origen o no.

De todos modos, yo recomendaría utilizar esta ecuación más general $rx+sy=t$ con tres coeficientes $(r,s,t)$ para cualquier línea del plano sin necesidad de casos especiales.

Answer 5

Supongamos que la elipse se decide por

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

Y $a >= b$ , $a^2 + b^2 = c^2$

Puntos dados $\phi_0$ , $\phi_1$ , $\phi_2$ decidido por la ecuación anterior y la relación de reflexión, enumeraremos los pasos del algoritmo

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

2. calcular $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. calcular $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. calcular $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. calcular $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. corrección del 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. calcular $\phi_2$ de $c_2$ y $s_2$ podemos utilizar arctan2 desde muchas librerías numéricas

La prueba del algoritmo anterior es elemental.

Lo más interesante es que

(1) Si el rayo pasa por uno de los focos, entonces pasará por otro y quedará atrapado dentro de la elipse

(2) Si el rayo pasa al lado de los dos focos, siempre estará al lado de los dos focos

(3) Si el rayo interseca el intervalo decidido por los dos focos, entonces siempre intersectará el intervalo

La elipse puede atrapar el rayo me deja asombrado en ese momento.

También quiero saber si se cumple algún tipo de ergodicidad para (2) y (3)