<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/laser-beams-and-elliptical-billiards-euler-problem-144/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://r.prevos.net/euler-problem-144/"> The Devil is in the Data</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Euler problem 144 has been one of the most fun to solve. The underlying problem is the pathway of the reflection of a laser inside an ellipse-shaped mirror. Before I delve into this problem, I like to share this delightful video from <a rel="nofollow" href="http://www.numberphile.com/" target="_blank">Numberphile</a> in which Alex Bellos demonstrates an elliptical billiards table. The billiards problem is mathematically equivalent to the laser problem. The reflection rule optics is the same as the bouncing rule in mechanics, but instead of using light, we use a billiard ball.</p>
<p><iframe allowfullscreen="true" src="https://www.youtube.com/embed/4KHCuXN2F3I?version=3&rel=1&%23038;fs=1&%23038;autohide=2&%23038;showsearch=0&%23038;showinfo=1&%23038;iv_load_policy=1&%23038;wmode=transparent" width="450" type="text/html"></iframe></p>
<p>This article outlines my solution to Euler problem 104 and simulates the elliptical pool table in the R language. You can find the code on the <a rel="nofollow" href="https://github.com/pprevos/r.prevos.net/blob/master/ProjectEuler/problem144.R" target="_blank">GitHub repository</a> for this website.</p>
<h2>Euler Problem 144 Definition</h2>
<p>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 <img alt="4x^2 + y^2= 100 " title="4x^2 + y^2= 100 " src="https://s0.wp.com/latex.php?latex=4x%5E2+%2B+y%5E2%3D+100+&bg=ffffff&%23038;fg=000&%23038;s=0" class="latex">. The section corresponding to <img alt="-0.01 leq times leq +0.01 " title="-0.01 leq times leq +0.01 " src="https://s0.wp.com/latex.php?latex=-0.01+%5Cleq+%5Ctimes+%5Cleq+%2B0.01+&bg=ffffff&%23038;fg=000&%23038;s=0" class="latex"> at the top is missing, allowing the light to enter and exit through the hole.</p>
<p><img src="https://i2.wp.com/projecteuler.net/project/images/p144_1.gif?resize=268%2C240&ssl=1" alt="Euler Problem 144" width="268" data-recalc-dims="1" height="240"></p>
<p>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 the 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 <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&%23038;fg=000&%23038;s=0" class="latex"> of the tangent line at any point <img alt="(x,y)" title="(x,y)" src="https://s0.wp.com/latex.php?latex=%28x%2Cy%29&bg=ffffff&%23038;fg=000&%23038;s=0" class="latex"> of the given ellipse is <img alt="m = -4x/y" title="m = -4x/y" src="https://s0.wp.com/latex.php?latex=m+%3D+-4x%2Fy&bg=ffffff&%23038;fg=000&%23038;s=0" class="latex">. The normal line is perpendicular to this tangent line at the point of incidence.</p>
<p>How many times does the beam hit the internal surface of the white cell before exiting?</p>
<h2>Proposed Solution to Euler Problem 144</h2>
<p>The first step was to rewrite the equation to use functions to generalise the problem. The general Cartesian equation for an ellipse is:</p>
<p><img alt="frac{x^2}{a^2} + frac{y^2}{b^2} = 1, a &lt; b " title="frac{x^2}{a^2} + frac{y^2}{b^2} = 1, a &lt; b " src="https://s0.wp.com/latex.php?latex=%5Cfrac%7Bx%5E2%7D%7Ba%5E2%7D+%2B+%5Cfrac%7By%5E2%7D%7Bb%5E2%7D+%3D+1%2C+a+%3C+b+&bg=ffffff&%23038;fg=000&%23038;s=2" class="latex"></p>
<p>The length of the axes for this problem is <img alt="a =5" title="a =5" src="https://s0.wp.com/latex.php?latex=a+%3D5&bg=ffffff&%23038;fg=000&%23038;s=0" class="latex"> and <img alt="b = 10" title="b = 10" src="https://s0.wp.com/latex.php?latex=b+%3D+10&bg=ffffff&%23038;fg=000&%23038;s=0" class="latex">. While the Project Euler description gives the formula for the slope of the tangent to the ellipse, I have generalised the code to reuse it for the elliptical billiards table. The slope of the tangent to an ellipse at point <img alt="(x,y)" title="(x,y)" src="https://s0.wp.com/latex.php?latex=%28x%2Cy%29&bg=ffffff&%23038;fg=000&%23038;s=0" class="latex"> is:</p>
<p><img alt="m=-frac{b^2x}{a^2y}" title="m=-frac{b^2x}{a^2y}" src="https://s0.wp.com/latex.php?latex=m%3D-%5Cfrac%7Bb%5E2x%7D%7Ba%5E2y%7D&bg=ffffff&%23038;fg=000&%23038;s=2" class="latex"></p>
<p>This first code snippet defines functions to draw an ellipse and calculate the bouncing angle. The last part of the code bounces the laser inside the cell until it exits through the top.</p>
<p>I sourced the formula to find the intersection between a line and an ellipse from <a rel="nofollow" href="http://www.ambrsoft.com/TrigoCalc/Circles2/Ellipse/EllipseLine.htm" target="_blank">Ambrsoft</a>. The equation has two possible solutions, one of which is the same as the original point.</p>
<pre>
plot_ellipse &lt;- function(a, b, colour = NA, line = &quot;black&quot;) {
    plot.new()
    plot.window(xlim = c(-a, a), ylim = c(-b, b), asp = 1)
    par(mar = rep(0,4))
    x &lt;- seq(-a, a, length = 200)
    y &lt;- sqrt(b^2 - (b^2 / a^2) * x^2)
    lines(x, y, col = line)
    lines(x, -y, col = line)
    polygon(x, y, col = colour, border = NA)
    polygon(x, -y, col = colour, border = NA)
}

bounce &lt;- function(coords) {
    x &lt;- coords$x
    y &lt;- coords$y
    ## Tangent to ellipse
    t &lt;- -(b^2 / a^2) * (x[2] / y[2])
    ## Deflection on sloping mirror y = mx + c
    dydx &lt;- diff(y) / diff(x)
    m &lt;- tan(pi - atan(dydx) + 2 * atan(t))
    c &lt;- y[2] - m * x[2]
    ## Determine intersection point
    ## Source: http://www.ambrsoft.com/TrigoCalc/Circles2/Ellipse/EllipseLine.htm
    x[1] &lt;- x[2]
    y[1] &lt;- y[2]
    x2 &lt;- (-a^2 * m * c + c(-1, 1) * (a * b * sqrt(a^2 * m^2 + b^2 - c^2))) /
          (a^2 * m^2 + b^2)
    x[2] &lt;- ifelse(round(x[1] / x2[1], 6) == 1, x2[2], x2[1])
    y[2] &lt;- m * x[2] + c
    return(data.frame(x, y))
}

# Initial conditions
a &lt;- 5
b &lt;- 10
x1 &lt;- 0
y1 &lt;- 10.1
x2 &lt;- 1.4
y2 &lt;- -9.6
answer &lt;- 0
plot_ellipse(a, b)
points(c(0,0), c(-c, c), pch = 19)
## Bounce laser breams
laser &lt;- data.frame(x = c(x1, x2), y = c(y1, y2))
while((laser$x[2] &lt; -0.01 | laser$x[2] > 0.01) | laser$y[2] &lt; 0) { ## Escape?
    lines(laser$x, laser$y, col = &quot;red&quot;, lwd = .5)
    laser &lt;- bounce(laser)
    answer &lt;- answer + 1
}
print(answer)
</pre>
<p>The result of this code is a pretty image of all the laser beams that have bounced around the mirror, which looks like the evil <a rel="nofollow" href="http://lotr.wikia.com/wiki/Eye_of_Sauron" target="_blank">Eye of Sauron</a> in Lord of the Rings.<br>
<a rel="nofollow" href="https://i0.wp.com/r.prevos.net/wp-content/uploads/sites/11/2018/05/problem144.png?ssl=1" target="_blank"><img data-orig-file="https://i0.wp.com/r.prevos.net/wp-content/uploads/sites/11/2018/05/problem144.png?fit=319%2C569&amp;ssl=1" data-orig-size="319,569" data-medium-file="https://i0.wp.com/r.prevos.net/wp-content/uploads/sites/11/2018/05/problem144.png?fit=168%2C300&amp;ssl=1" data-attachment-id="1245" data-image-meta="{" data-image-description="" src="https://i0.wp.com/r.prevos.net/wp-content/uploads/sites/11/2018/05/problem144.png?resize=319%2C569&ssl=1" data-permalink="https://r.prevos.net/euler-problem-144/problem144/" data-large-file="https://i0.wp.com/r.prevos.net/wp-content/uploads/sites/11/2018/05/problem144.png?fit=319%2C569&amp;ssl=1" alt="Solution to Euler Problem 144" data-image-title="problem144" sizes="(max-width: 319px) 100vw, 319px" width="319" data-comments-opened="1" class="wp-image-1245" data-recalc-dims="1" height="569"></a></p>
<h2>Elliptical Pool Table</h2>
<p>We can use the solution to Euler problem 144 to play billiards on an elliptical billiards table. To close the article, we return to the elliptical pool table demonstrated by Alex Bellos. This code draws the pool table to the dimensions mentioned in the video. We know that the table has an eccentricity of <img alt="e = 0.43" title="e = 0.43" src="https://s0.wp.com/latex.php?latex=e+%3D+0.43&bg=ffffff&%23038;fg=000&%23038;s=0" class="latex"> and a long axis of <img alt="a = 130" title="a = 130" src="https://s0.wp.com/latex.php?latex=a+%3D+130&bg=ffffff&%23038;fg=000&%23038;s=0" class="latex"> cm. The code defines the short axis (<img alt="b" title="b" src="https://s0.wp.com/latex.php?latex=b&bg=ffffff&%23038;fg=000&%23038;s=0" class="latex">) and the distance of the focal points from the centre.</p>
<p>The code selects a random starting point and angle of the shot. The code first determines whether the line passes through the pocket. If this is not the case, the algorithm then finds the place where the ball hits and keeps bouncing it until it falls into the pocket or the ball bounces 100 times.</p>
<p>Elliptical billiard tables have four possible outcomes. Any ball the pass through a focal point will fall into the pocket, ending the simulation. Any ball that passes outside the focal points will bounce around, and the combined trajectories form an ellipse. When the ball moves between the foci, the result is a hyperbola. Lastly, there are some unique circumstances which result in a regular polygon.</p>
<p><a rel="nofollow" href="https://r.prevos.net/?attachment_id=1290" target="_blank"><img data-orig-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2018/05/epilliptical_billiards.png?fit=3208%2C949&amp;ssl=1" data-orig-size="3208,949" data-medium-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2018/05/epilliptical_billiards.png?fit=300%2C89&amp;ssl=1" data-attachment-id="1290" data-image-meta="{" data-image-description="" src="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2018/05/epilliptical_billiards.png?resize=584%2C173&ssl=1" data-permalink="https://r.prevos.net/euler-problem-144/epilliptical_billiards/" data-large-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2018/05/epilliptical_billiards.png?fit=584%2C173&amp;ssl=1" alt="Elliptical billiards: Euler problem 144" data-image-title="epilliptical_billiards" sizes="(max-width: 584px) 100vw, 584px" width="450" data-comments-opened="1" class="wp-image-1290" data-recalc-dims="1"></a>If simulations are not enough for you, then head over to the <a rel="nofollow" href="http://www.instructables.com/id/Elliptical-Pool-Table/" target="_blank">Instructables</a> website to find out how you can construct an elliptical billiards table. There is even a patent for an <a rel="nofollow" href="https://patents.google.com/patent/US5226644" target="_blank">elliptical pocket billiard table</a>, with the pockets at the edge.</p>
<pre>
## Elliptical pool table
## https://www.youtube.com/watch?time_continue=54&amp;v=4KHCuXN2F3I
e &lt;- 0.43
a &lt;- 130
b &lt;- a * sqrt((1 + e) * (1 - e)) # a > b
f &lt;- sqrt(a^2 - b^2)
plot_ellipse(a, b, &quot;darkgreen&quot;, NA)
points(-f, 0, pch = 19, cex = 2)
points(f, 0, pch = 19, col = &quot;grey&quot;)

## Simulate random shot
angle &lt;- runif(1, 0, 2 * pi)
x1 &lt;- runif(1, -a, a)
ymax &lt;- sqrt(b^2 - (b^2 / a^2) * x1^2)
y1 &lt;- runif(1, -ymax, ymax)

## First shot
m &lt;- tan(angle)
c &lt;- y1 - m * x1
x2 &lt;- (-a^2 * m * c + c(-1, 1) * (a * b * sqrt(a^2 * m^2 + b^2 - c^2))) / (a^2 * m^2 + b^2)
y2 &lt;- m * x2 + c
x2 &lt;- x2[which(((x2 - x1) &lt; 0) == (cos(angle) &lt; 0))]
y2 &lt;- y2[which(((y2 - y1) &lt; 0) == (sin(angle) &lt; 0))]
shot &lt;- (data.frame(x = c(x1, x2), y = c(y1, y2)))

## Bounce ball
for (i in 1:100){
    dydx &lt;- diff(shot$y) / diff(shot$x)
    if (all.equal(dydx, (shot$y[1] - 0) / (shot$x[1] - -f)) == TRUE) {
        shot[2, ] &lt;- c(-f, 0)
    }
    lines(shot$x, shot$y, col = &quot;yellow&quot;, lwd = 1)
    if (shot[2,2] == 0) break
    shot &lt;- bounce(shot)
}
points(x1, y1, pch = 19, col = &quot;blue&quot;, cex = 1.8)
</pre>
<p>The post <a rel="nofollow" href="https://r.prevos.net/euler-problem-144/" target="_blank">Laser Beams and Elliptical Billiards: Euler Problem 144</a> appeared first on <a rel="nofollow" href="https://r.prevos.net/" target="_blank">The Devil is in the Data</a>.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://r.prevos.net/euler-problem-144/"> The Devil is in the Data</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/Zrz9JJ74lKg" height="1"></a></div>