Bouncing Photon

This is a very non-trivial problem, so I'll break it down into pieces. The general outline of my approach is this:

1) Record the position and velocity values
2) Determine the travel distance to the next collision point with one of the circles
3) Update the position (corresponding to the next collision point)
4) Update the velocity to reflect the bounce
5) Go back to Step 1 and repeat this process until enough bounces have happened to account for at least 10 seconds of time passage.

Some of these steps require some additional explanation

Derivation for Step 2:

2a) Start at point $\vec{P} = (x,y)$ with velocity $\vec{v}$ . Note that for this problem, $\vec{v}$ is a unit-vector.
2b) Equation of circle with center $(m,n):$ $(x-m)^{2} + (y-n)^{2} = R^{2}$
2c) If we start at $\vec{P}$ with velocity $\vec{v}$ and travel a distance $\alpha$ , this can be represented as:

$(x + \alpha v_x - m)^{2} + (y + \alpha v_y - n)^{2} = R^{2} \\ \alpha^{2} v_x^{2} + 2v_x(x-m)\alpha + (x-m)^{2} + \alpha^{2} v_y^{2} + 2v_y(y-n)\alpha + (y-n)^{2} - R^{2} = 0$

$a = v_x^{2} + v_y^{2} \\ b = 2v_x(x-m)+2v_y(y-n) \\ c = (x - m)^{2} + (y - n)^{2} - R^{2}$

$\alpha = \frac{-b \pm \sqrt{b^{2} -4ac}}{2a}$

2d) Search over a space of integer (m,n) pairs within the range of (-10,10), and find the minimum positive $\alpha$ value, which is the travel distance to the boundary of the nearest circle. Also store the (m,n) coordinates so that the normal and tangential vectors can be calculated for Step 4.

Processing for Step 3

$x_{new} = x_{old} + \alpha v_x \\ y_{new} = y_{old} + \alpha v_y$

Derivation for Step 4:

Calculate the x and y coordinates of the normal and tangential vectors:

$n_x = x - m \\ n_y = y - n \\ t_x = -n_y \\ t_y = n_x$

Resolve the velocity just prior to the bounce into a linear combination of the normal and tangential vectors:

$v_x = \sigma n_x + \gamma t_x \\ v_y = \sigma n_y + \gamma t_y$

Rearranging the second equation gives: $\gamma = \frac{v_y - \sigma n_y}{t_y}$

Substituting into the first equation and solving for $\sigma$ gives:

$\sigma = \frac{v_x - \frac{v_y t_x}{t_y}}{n_x - \frac{n_y t_x}{t_y} }$

Plugging the value of $\sigma$ back into the $\gamma$ equation allows us to solve for that parameter. In order to model the bounce, we reverse the normal component of the velocity and preserve the tangential component as follows:

$v_{x_{new}} = -\sigma n_x + \gamma t_x \\ v_{y_{new}} = -\sigma n_y + \gamma t_y$

Python Code:

import math

x = 0.5
y = 0.1

vx = 1.0
vy = 0.0

R = 1.0/3.0

m_store = 0
n_store = 0

for q in range(0,20):  # For 20 bounces

    print x,y

    alpha_min = 10.0**12.0


    for m in range(-10,10):     # Find location of next bounce
        for n in range(-10,10):

            a = vx**2.0 + vy**2.0
            b = 2.0*vx*(x-m) + 2.0*vy*(y-n)
            c = (x-m)**2.0 + (y-n)**2.0 - R**2.0

            alpha_plus = 10.0**12.0
            alpha_minus = 10.0**12.0

            if b**2.0 > 4.0*a*c:

                alpha_plus = (-b + math.sqrt(b**2.0 - 4.0*a*c))/(2.0*a)
                alpha_minus = (-b - math.sqrt(b**2.0 - 4.0*a*c))/(2.0*a)


            if alpha_plus < 0:
               alpha_plus = 10.0**12.0

            if alpha_minus < 0:
                alpha_minus = 10.0**12.0

            alpha_test = min(alpha_plus,alpha_minus)

            if (alpha_test < alpha_min) and (alpha_test > 0.01):
                alpha_min = alpha_test
                m_store = m
                n_store = n


    x = x + alpha_min * vx      # Update coordinates
    y = y + alpha_min * vy

    # Re-calculate velocity

    nvec_x = x - m_store
    nvec_y = y - n_store

    tvec_x = -nvec_y
    tvec_y = nvec_x

    sigma = (vx - vy*tvec_x/tvec_y) / (nvec_x - nvec_y*tvec_x/tvec_y)

    gamma = (vy - sigma*nvec_y)/tvec_y

    vx = -sigma * nvec_x + gamma * tvec_x
    vy = -sigma * nvec_y + gamma * tvec_y

Code Output for 20 Bounces (these are the x and y coordinate pairs associated with each collision):

0.5 0.1
0.682020266194 0.1
-0.669499718785 1.04336675256
-0.123873465869 1.6905384102
-0.147815682201 1.29876685761
-0.797941821159 1.73488945047
-0.279026615343 0.182360245232
-0.699476804778 0.144211373493
-0.158068075392 0.706528375046
-0.274764202113 0.188721340468
-0.676257918623 -0.07938624476
-0.275111646835 -0.81178551361
-0.785422442156 -0.744916518017
-0.917073039951 -0.322853264515
-0.838751893192 -0.708263545024
0.693386880249 -0.130765079085
-0.726578772593 -0.809338143526
-0.182097464933 -0.279198181183
-0.744349303805 -1.78609854454
1.87688240854 -1.30976308654

I then imported the data into Excel. Knowing the present and previous coordinates, as well as the velocity magnitude (equal to one), it is simple to associate a time value with each collision. A final linear interpolation gave the coordinates associated with a time value of 10 seconds. At 10 seconds, the distance from the origin is approximately 0.995263.