Project Euler 144: Investigating multiple reflections of a laser beam.

Project Euler 144: Investigating multiple reflections of a laser beam.

Written by on 13 April 2013

Topics: Project Euler

Problem 144 of Project Euler is once again a geometry problem, just like the previous. However, it is completely different. The problem reads

In laser physics, a “white cell” is a mirror system that acts as a delay line for the laser beam. The beam enters the cell, bounces around on the mirrors, and eventually works its way back out.

The specific white cell we will be considering is an ellipse with the equation 4x2 + y2 = 100

The section corresponding to -0.01  x  +0.01 at the top is missing, allowing the light to enter and exit through the hole.

The light beam in this problem starts at the point (0.0,10.1) just outside the white cell, and the beam first impacts the mirror at (1.4,-9.6).

Each time the laser beam hits the surface of the ellipse, it follows the usual law of reflection “angle of incidence equals angle of reflection.” That is, both the incident and reflected beams make the same angle with the normal line at the point of incidence.

In the figure on the left, the red line shows the first two points of contact between the laser beam and the wall of the white cell; the blue line shows the line tangent to the ellipse at the point of incidence of the first bounce.

The slope m of the tangent line at any point (x,y) of the given ellipse is: m = -4x/y

The normal line is perpendicular to this tangent line at the point of incidence.

The animation on the right shows the first 10 reflections of the beam.

How many times does the beam hit the internal surface of the white cell before exiting?

We will simply brute force our way through this problem, by calculating the laser beams path through the cell, and check if it hits the exit. In order to do that, we need to calculate how the laserbeam reflects. Once we know the angle of the reflecting beam, we can calculate the corresponding line, since we have the point of reflection. Once we have a line parameterization of the reflecting line, it is simply a matter of finding out where the line and the ellipse intersect. This will be the next point out laser beam hits. Confused yet? Don’t be. I will elaborate on it.  Lets start by finding the slope of the reflecting beam.

Slope of reflecting beam

Let us draw the situation and put some names on the different angles.

problem144

As you can see we have a horizontal line which we will base everything off. Then we have the ellipse tangent where the slope (which we will denote ) is already know from the problem description. Then we have the line AO, which we also know the slope of (denoted ). And then we have the line OB which is the one we are interested in finding the slope of (denoted ). So now we have already established a relation between the slopes and the angles.

One thing we can deduce from the image is that

 and  

We are interested in rewriting these relations into something based on the slopes of the lines.. Before continue let me just mention that we have

A relation which we will use in just a few seconds. For the first relation we have

For the second relation we get that

In this equation we can isolate such that we have

where tan(a) is given by the first equation.

So know we have the slop of the line OB, and since we know a point on the line (O) we can describe it as

where is the slope and d is the intercept. Which having a point becomes

Intersecting a line and an ellipse

Now that we have a full description of the line OB all we need is to figure out where the line intersects with the ellipse. The ellipse have the equation

Plugging the line equation into this gives us that

So all we have to do is solve this equation and take the one of the two solutions which is not , then we have the x value of B and the y value follows from the line equation.

Coding it all

For once I have actually managed to comment my C# code, so you should be able to read it based on the problem description. At least I wont comment further on it.

int result = 0;

double xA = 0.0;
double yA = 10.1;

double xO = 1.4;
double yO = -9.6;

while(xO > 0.01 || xO < -0.01 || yO < 0){

    //Calculate the slope of A
    double slopeA = (yO - yA) / (xO - xA);

    //Calculate the slope of the ellipse tangent
    double slopeO = -4*xO/yO;

    //Calculate the slope of B
    double tanA = (slopeA - slopeO)/(1+slopeA*slopeO);
    double slopeB = (slopeO- tanA)/ (1+ tanA*slopeO);

    //calculate intercept of line B
    double interceptB = yO - slopeB * xO;

    //solve the quadratic equation for finding
    // the intersection of B and the ellipse
    // a*x^2 + b*x + c = 0
    double a = 4 + slopeB*slopeB;
    double b = 2 * slopeB * interceptB;
    double c = interceptB * interceptB - 100;

    double ans1 = (-b + Math.Sqrt(b * b - 4 * a * c)) / (2 * a);
    double ans2 = (-b - Math.Sqrt(b * b - 4 * a * c)) / (2 * a);

    xA = xO;
    yA = yO;

    //Take the solution which is furtherst from x0
    xO = (Math.Abs(ans1 - xO) > Math.Abs(ans2 - xO)) ? ans1 : ans2;
    yO = slopeB * xO + interceptB;

    result++;
}

This gives me the answer in 0.02ms or so. So this is by far fast enough. As always you can find the full source code here.

12 Comments For This Post I'd Love to Hear Yours!

  1. Scott Dennison says:

    Hi.

    One suggestion and one question regarding this post, which I now mostly understand, after many failed attempts over the past few years.

    First the suggestion.
    When you are explaining tab(a-b), I thought that ‘a’ was referring to the ‘a’ in the diagram and earlier calculations, not just a generic variable in the tab(a-b) formula, and it wasn’t for a while that it twigged. It might be worth explaining that somewhere.

    Secondly, the question.
    I’m probably being dumb, but there again as I type this I’m still half asleep, but I can see how you extract a=(a1+pi)-a2 from the diagram, but not a=a0-a1. How do you do this?

  2. Kristian says:

    Hi Scott

    I take the first part as a huge compliment, so thank you. For the second part, I have changed it to p and q to avoid confusing now. Thanks for that suggestion.

    The question you have. Take a look at the illustration again. It is apparent in the right side of the illustration

    a is the angle between the ellipse tangent and the laser beam
    a0 is the line between the laser beam and the horizontal line
    a1 is the angle between the ellipse tangent and the horizontal line

    so if we take the a0 – a1 then we get

    (laser beam – horizontal) – (ellipse – horizontal) -> laserbeam – horizontal

    which is exactly the angle a.

  3. Scott Dennison says:

    Hi Kristian

    🙂

    My attempts in the past have all ended in failure, either from taking so long I killed it, giving me the wrong answer, or working only if I multiple something by -1 sometimes and not others, something I could never understand. I still have the code floating around somewhere (One of my old attempts is http://pastebin.com/BiHj1ACu). Yet this solution using tan(a-b) just works and is easy to understand.

    Using your explanation of a=a0-a1, I can now see where I was going wrong when trying to understand the diagram you provided.
    Firstly, I was only paying attention to the a on the left, not noticing a on the right which is the same angle due to the nature of the problem.
    Secondly, I was only seeing a0 as the portion between the laser beam and the ellipse tangent (which happens to actually be the right a), not realizing the mark on the diagram extended through the ellipse tangent to the horizontal line.

  4. Jean-Marie Hachey says:

    Table 1
    Discrete variations in the limits of the parameters j and k and its effects on the number of reflections of the laser beam.

    White cell equations:
    4*(x^2) + y^2 = 100 [1]
    j*(x^2) + y^2 = k [2]

    http://img4.hostingpics.net/pics/556817pe144ihtable1fig12.jpg

    As j is kept constant at 4, we observe that
    during the variation of k from 99 to 100.2, refl.
    will vary between 19 and 896. The value of 354
    reflections is reached when k is closer to 100.
    Similar type of observations when k=100 and j varies
    between 3.9 and 4.2. The value of 354 reflections is
    obtained when j is closer to 4.

    ___

    Sources:
    1) Project Euler 144
    Kristian’s algorithm
    http://www.mathblog.dk/files/euler/Problem144.cs
    2) Microsoft Visual C# 2010 Express
    3) Microsoft Office Excel (2007)

  5. Mayank says:

    //Take the solution which is furthest from x0
    xO = (Math.Abs(ans1 – xO) > Math.Abs(ans2 – xO)) ? ans1 : ans2;
    yO = slopeB * xO + interceptB;

    Why you took furthest ?

  6. Kristian says:

    There are two intersections between the oval and the line. One of these intersections is the one where the beam is hitting after the reflection. The other is the point of reflection.

    Due to numerical issues I cannot be sure that x0 does not vary slightly from the answer. So instead I just take the point which is farthest away.

  7. Oren says:

    Hi Kristian,

    You don’t need to solve a quadratic equation. Since you know (xO,yO) is a solution, you can factor out x-xO and remain with a linear equation for the other solution. I programmed it and it works.

    https://sites.google.com/site/orenlivne/

    Oren

  8. Kristian says:

    That sounds right. However, I cannot wrap my head around how to do it immidatly. Can you show how you did it?

  9. Oren says:

    Hi Kristian,

    Since OB passes through (xO,yO), write it as y=c(x-xO)+yO. Hence

    4x^2 + (c(x-xO)+yO)^2 = 100
    4x^2 + c^2(x-xO)^2 + 2c(x-xO)yO + yO^2 – 100 = 0
    But yO^2 – 100 = -4xO^2, since (xO,yO) also lies on the ellipse.
    4x^2 + c^2(x-xO)^2 + 2c(x-xO)yO – 4xO^2 = 0
    4(x^2-xO^2) + c^2(x-xO)^2 + 2c(x-xO)yO = 0
    4(x-xO)(x+xO) + c^2(x-xO)^2 + 2c(x-xO)yO = 0
    Now we can finally factor out x-xO.
    4(x+xO) + c^2(x-xO) + 2c yO = 0
    (c^2+4)x + 4 xO-c^2 xO + 2c yO = 0
    x = (xO(c^2-4) – 2c yO)/(c^2+4).

    Please double-check this computation as I have a tendency to make mistakes. Thank you very much for putting this website together. It saved me many sleepless nights over several of the Euler problems.

  10. Kristian says:

    Ah yes, now I see it. And that means we do indeed reduce x to be a linear equation based on values we already know. Nice one.

  11. Ayush Bahuguna says:

    I tried this problem instead in parametric form of an ellipse (x=acos(t),y=bsin(t)) and it’s MUCH simpler. No need of a quadratic equation. Finding the third parametric angle if first two parametric angles are known doesn’t even need inverse trigonometric functions. It’s almost elementary after you derive the expression.

  12. Axel says:

    It seems to me that you went of on a tangent here. (excuse the pun)
    This problem however is very well suited to be solved in vector form.

    Let u_k be the direction of the incomming ray and x_k be the point where it strikes the ellipse.
    u_0 = (1.4, -19.7), x_0 = (1.4, -9.6)
    Since we are told the formula for the slope of the tangent is -4*x/y it is easy to show that the direction of the normal vector as point x_k can be written as
    n_k = (-4*x_k[0], -x_k[1])
    By the law or reflection the vector differnce between the incoming and outgoing rays must be a multiple of the normal vector:
    u_(k+1) = u_k + p*n_k
    where p is an unknown scalar.

    The fact that both rays must make the same angle with the tangent line is equivalent to requirering that thier inner products with the normal vector be of same size but opposite sign. This can be used to solve for p:
    p = -2*dot(u_k, n_k)/dot(n_k, n_k)

    Next we need to find the new point where the outgoing ray strikes the ellipse. This point must lie on the line:
    x_(k+1) = x_k + q*u_(k+1)
    where q is an unknown scalar.
    Define a matrix:
    A = [[2,0],[0,1]]
    Note that the function f(x_k) = 4*x_k[0]^2 + x_k[1]^2
    can be written as: f(x_k) = dot(A*x_k, A*x_k)
    Therfore:
    f(x_(k+1)) = dot(A*(x_k + q*u_(k+1)), A*(x_k + q*u_(k+1)))
    = dot(A*x_k + q*A*u_(k+1), A*x_k + q*A*u_(k+1))
    = dot(A*x_k, A*x_k) + 2*q*dot(A*x_k, A*u_(k+1))
    + q^2*dot(A*u_k(k+1), A*u_(k+1))
    The first term is the function value of the current point, which should be unchanged. Therefore the sum last two terms must equal zero.
    The trivial solution, q = 0, is the current point.
    The non-trivial solution gives the new point:
    q = -2*dot(A*u_k(k+1), A*u_(k+1))/dot(A*x_k, A*u_(k+1))

    This became a lot of text, but for the code it boils down to needing just five lines inside the loop:
    calculate n_k at x_k
    calculate p
    update u_k
    calculate q
    update x_k

Leave a Comment Here's Your Chance to Be Heard!

You can use these html tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

You can use short tags like: [code language="language"][/code]

The code tag supports the following languages: as3, bash, coldfusion, c, cpp, csharp, css, delphi, diff,erlang, groovy, javascript, java, javafx, perl, php, plain, powershell, python, ruby, scala, sql, tex, vb, xml

You can use [latex]your formula[/latex] if you want to format something using latex

This site uses cookies.