ProjectEuler 144: Investigating multiple reflections of a laser beam

tl;dr; you can stop reading, if you aren't interested in programming or math.

In this post, I want to show you importance of programming in math field. For a year I was afraid to solve one geometrical problem. So today decided to overcome my fear.

I don't want to make you afraid of this post. So let's start with something simple.

1. Draw an ellipse

We are given an ellipse with the equation: 4*x^2+y^2 = 100

So we want to visualise the ellipse. In the end, we want to plot ellipse as shown in the following picture:

There is a great python plotting library for that: matplotlib.

Firstly, notice that x needs to be between -5..5 inclusive, otherwise 4*x^2 + y^2 will be greater that 100.

If we know x, then y = sqrt(100-4*x^2) or y = -sqrt(100-4*x^2)

Let's take 5 points from [-5..5] and draw line between them.

import numpy as np
import matplotlib.pyplot as plt
import math

plt.figure(figsize=(4,5))

X = np.linspace(-5, 5, 5) 
Y = [math.sqrt(100-4*x*x) for x in X]
plt.plot(X, Y, color="red")
plt.show()

We will get something like that:

Have you noticed 5 points? Let's draw more points, let's say 100 points.

import numpy as np
import matplotlib.pyplot as plt
import math

plt.figure(figsize=(4,5))

X = np.linspace(-5, 5, 100) 
Y = [math.sqrt(100-4*x*x) for x in X]
plt.plot(X, Y, color="red")
plt.show()

Have you noticed which number has changed in the code?

Looks like we have drawn upper part of the desired ellipse, lower part can be drawn similarly:

import numpy as np
import matplotlib.pyplot as plt
import math

plt.figure(figsize=(4,5))

X = np.linspace(-5, 5, 100) 
Y = [math.sqrt(100-4*x*x) for x in X]
plt.plot(X, Y, color="red")
Y = [-math.sqrt(100-4*x*x) for x in X]
plt.plot(X, Y, color="blue") 
plt.show()

Let's try to draw a line between points (0, 10.1) and (1.4, -9.6):

import numpy as np
import matplotlib.pyplot as plt
import math

plt.figure(figsize=(4,5))

X = np.linspace(-5, 5, 100) 
Y = [math.sqrt(100-4*x*x) for x in X]
plt.plot(X, Y, color="red")
Y = [-math.sqrt(100-4*x*x) for x in X]
plt.plot(X, Y, color="blue") 
plt.plot([0, 1.4], [10.1, -9.6], color="black")
plt.show()

2. Problem Statement

Now it's time to solve the real problem. You can read the statement from here: ProjectEuler 144.

Shortly, you are given the equation of an ellipse: 4*x^2+y^2=100. The laser beam enters the ellipse at point (0, 10.1) and then touches point (1.4, -9.6), after that reflects, and so on... See the animation below or read the full statement:

3. Reflecting one segment

If we will learn to reflect the laser beam at one point, then we solved the problem, because after that we can just simulate it.

The picture above gives some draft ideas about how we are going to reflect the laser beam, i.e:

a. Find the normal vector(yellow segment) to tangent line(brown) at point q.

b. Find angle between laser beam and normal vector (i.e to find measure of the angle A). A is negative if laser beam is located in the left side of the normal vector.

c. Rotate the laser beam to 2*A angles, in other words rotate the point p around the point q to an angle 2*A.

d. Finally, find a point t (look to the picture), which will be located in the line formed by points q and rotated p.

3.a. Finding the normal vector

Let's recall lectures from calculus, specially gradient vector. Don't be afraid, it is just vector formed by partial derivatives at points x and y. In our case, gradient vector or normal vector will be (8*x, 2*y). I am not going to prove it, but I am going to use benefits of python to verify visually for correctness. Let's draw normal vector at random points.

Seems like something reasonable to continue.

3.b. Finding an angle between two vectors

From linear algebra we can remember that vector product of two vectors (x1, y1) and (x2, y2) is equal to x1*y2-y2*x1.

On the other hand, it equals  to  sqrt(x1^2+y1^2)*sqrt(x2^2+y2^2)*sin(A),

i.e

x1*y2-y2*x1=sqrt(x1^2+y1^2)*sqrt(x2^2+y2^2)*sin(A)

sin(A) = (x1*y2-y2*x1)/(sqrt(x1^2+y1^2)*sqrt(x2^2+y2^2))

A = arcsin((x1*y2-y2*x1)/(sqrt(x1^2+y1^2)*sqrt(x2^2+y2^2)))

def veclen(u):
    return (u[0] * u[0] + u[1] * u[1]) ** 0.5

def angle(u, v):
    dot = u[0] * v[1] - u[1] * v[0]
    sin = dot / veclen(u) / veclen(v)
    return math.asin(sin)

3.C. Rotation of a vector

You can read it from here. To be honest, I just copied the formula from wikipedia :D.

def rotate_vector(v, angle, origin):
    x, y = v
    x = x - origin[0]
    y = y - origin[1]
    cos_theta = math.cos(angle)
    sin_theta = math.sin(angle)
    nx = x*cos_theta - y*sin_theta
    ny = x*sin_theta + y*cos_theta
    nx = nx + origin[0]
    ny = ny + origin[1]
    return np.array([nx, ny])

After this step we will get this picture, now the last step is to find point t.

3.D Finding the point t.

Now we have points q and rotated p. We want to find point t, such that it belongs to a line formed by points (q, rotated p) and belongs to the ellipse.

One way to solve, is to parameterise the line equation. Lets denote rotated p as pr.

x=(pr.x-q.x) * k + q.x

y=(pr.y-q.y) * k + q.y

This is parametrised equation of the line formed by point points q and pr.

We want to find k, such that:

4*((pr.x-q.x) * k + q.x) ^ 2 + ((pr.y-q.y) * k + q.y) ^ 2 = 100

In a result, we can simplify it to a form a*k^2+b*k+c=0 and solve quadratic equation.

4. Final draw

First few reflections are drawn in the picture below, seems like the solution worked out.

You can find full code here: pastebin

By the way, we might draw nothing and solve the problem from projecteuler (it just needs a number). But by drawing, it is easier to see that you are in the right direction.