1/17/19 - I've had this peach of a problem in the back of my mind for a couple of years now, here's an attempt to get some of the math laid out.

### Bucket problem

Jacob is camping in the mountains with his friends and, like a responsible citizen, fills up a bucket with water before he starts the fire for that night. When he goes to the water spout to fill up his big orange Home Depot bucket (which has a glossy reflective interior), he notices that his headlamp makes a very strange reflection at the bottom of his bucket! Perplexed, he decides to forget about it, saying "oh well I probably need to know like abstract linear calculus to figure out why this happens." Which was wrong of him.

a piece of paper laid down on the bottom of a metal cylindrical rice cooker (contrast slider --------->+ to make it look cooler)

### Summary of approach

First, I use a computer simulation, programmed in MATLAB, to model the reflection of light inside the bucket.

Next, I introduce the idea of the catacaustic as the mathematical basis for the shapes formed on the bucket.

Last, I solve for the catacaustic under various configurations of the light source and bucket.

### Theoretical set-up

As with any problem, we need to make some simplifying assumptions! Suppose the light source is a height $h$ above the base of the bucket (radius $R$) and the bucket's walls extend arbitrarily high. Furthermore, the light source is assumed to be very small, scattering light evenly in all directions. Looking directly down into the bucket, the light source is offset from the center-axis by a distance $d.$

### Computer simulation

I wrote a simple piece of code to simulate the reflection of light inside the cylinder, plotting where the light hits. I used MATLAB, and here is the pseudocode
  		
loop for # of iterations:
set initial light position (always same)
set initial light direction (randomized)
loop until broken:
store current position
advance position using light direction and large step size
check if the light is outside the bucket walls, if so
reset position to the stored previous position
loop until broken:
store current position
advance position using light direction and small step size (more precision)
check if the light is outside the bucket walls, if so
store previous position in original variable
recalculate light direction (explained below)
break
check if the light is below the bucket floor, if so
reset position to the stored previous position
loop until broken:
store current position
advance position using light direction and small step size
if the light is below the bucket floor, if so
store previous position as final (i.e. where the light hits the bucket floor)
break
break parent loop


The direction of the light can be randomized using spherical coordinates: the direction can be described using two angles $\theta$ and $\alpha,$ the horizontal angular displacement from the $x$ axis and the angle of depression, respectively. So it would look something like

theta = rand*2*pi
alpha = rand*pi/2
v = [cos(theta)*cos(alpha), sin(theta)*cos(alpha), -sin(alpha)]

Changing the light direction based on how it reflects off the circular interior looks tricky but is quite easy. Compressing the light position $\vec{x}$ and light direction vectors $\vec{v}$ into their solely their horizontal components, you can confirm that the new light direction is $$\vec{v} - \frac{2\vec{v}\cdot \vec{x}}{\text{norm}(x)}\vec{x}.$$ Here are the results of the program (click to enlarge, arrow keys to compare):

Click these pretty pictures!

(I couldn't get it as a vectorized image because MATLAB crashed trying to save the figure, had to screenshot it.) This heatmap scatter plot is the result of the program iterated 1 million times (for images 2&3 only a fraction of these were plotted) with the initial parameters R = 1, d = 0.8, h = 2.

Hopefully these three different visualizations of the data can give you some idea of the complex and interesting shapes, each of different intensity, formed in this extraordinarily simple configuration.

### Mathematics

To understand the geometry of the shapes produced on the bottom of the bucket, it is first necessary to reduce the problem into a two-dimensional one. How do we determine where the light ray is stopped by the bottom of the bucket? Imagine it is simply travelling and reflecting in the plane of the circular cross section of the bucket. Then after a time $t$ that is dependent on the height $h$ and angle of depression $\alpha,$ it will be stopped, and we record the point where it stops as a point of light on the bottom of the bucket.

So the shapes formed in the bucket are closely related to the possible configurations of reflected light rays within the circle! Given the set of all ray configurations (if you traced a ray and all its reflections within the circle), the "density" of lines in a certain area will increase the probability that the light ray is "stopped" within that area*. If a point on the circle has a high likelihood of being the location where a light ray is "stopped," then it will form a shape on the bottom of the bucket. So the problem of the ray "density" is essentially* equivalent to the problem of the shapes on the bucket.

*This is not exactly true, the probability that a light ray is "stopped" only after a very long time is exceedingly small because $\alpha$ would have to very near zero when in fact the light rays are distributed uniformly over $\alpha$. In other words, the "probabilities" we are speaking of are tied to the distribution of $\alpha$, which we will ignore here. This does not affect our analysis of the shapes formed on the bottom of the bucket, but it would affect the intensities of the shapes formed on the bottom of the bucket.

(Generated in Geogebra with $d=0.8R,$ you can very clearly see that the shapes formed match those that form on the bottom of the bucket. The first image is the set of all possible rays following one reflection inside the circle, the second image is the set of all possible rays following exactly two reflections inside the circle, and the third image is the set of all possible rays following exactly three reflections inside the circle. You could continue on, but the successive images are not as visible on the bucket because it is not likely that a ray will have the space to reflect so many times.)

The very loose idea of the ray "density" can be formalized with the idea of envelopes. For a set of curves (lines in this case), the envelope is a curve that is tangent to all the curves and formed by the points of tangency:

(Image from wikipedia) a simple example of an envelope for a family of curves.

In the three images presented above, the envelopes are very clearly visible. Mathworld (link below) gives the mathematical definition as: for a one-parameter family of curves given by $U(x,y,c) = 0,$ the envelope is the set of points $(x,y)$ such that there exists a number $c$ where $\dfrac{\partial U}{\partial c} = 0$ and $U(x,y,c) = 0.$

We can use the simple example of the gif above. The set of curves is given by $$y + \dfrac{k}{1-k}x - k= 0$$ for $k\in(0,1)$ and $x\in(0,1).$ Simultaneously solving $$\dfrac{\partial}{\partial k}\left(y+\dfrac{k}{1-k}x-k\right) = 0\qquad\text{ and }\qquad y + \dfrac{k}{1-k}x - k = 0,$$ we find $y=(1-\sqrt{x})^2,$ which is the red curve.

We are looking for a way to express the set of all possible ray reflections in a similar way as the range of a function of $x,y,$ and some other parameter. We'd like to be able to find an algebraic way of expressing all possible ray reflections.

Considering the angles involved in the problem, polar coordinates may seem like a good idea, but it didn't work out so well when I tried it. Here's a method using the simple idea behind reflection and some degree of wishful thinking.

The tricky part of the problem is algebraically expressing a reflection. The reflected ray is simple the reflection of the original ray over the line connecting the center and the point of reflection. We think: this would be pretty easy if the line connecting the center and the point of reflection were the $x$-axis, then we just negate the function, yeah? Once that goes through your mind, the solution is pretty obvious and it just becomes algebra.

(Another way to do this is with vectors, using the dot product. It could even be easier, but the result is the same nonetheless.)

TLDR: switch coordinate frame to make reflection easy, then switch back.

Suppose the ray begins at the point $(0,-d)$ and it is directed towards a point on the circle given by $(R\cos\theta,R\sin\theta),$ then the line joining these two is $$y = -d + \dfrac{R\sin\theta+d}{R\cos\theta}x.$$ We use the rotation formulas: \begin{align*} x &= x'\cos\theta - y'\sin\theta\\ y &= x'\sin\theta + y'\cos\theta. \end{align*} So now the ray can be reexpressed as \begin{align*} x'\sin\theta + y'\cos\theta &= -d + \left(\sin\theta + \dfrac{d}{R}\right)(x'-y'\tan\theta)\\ y'\left(\cos\theta + \theta\theta\sin\theta + \dfrac{d\tan\theta}{R}\right) &= \dfrac{dx'}{R} - d\\ y'\left(\dfrac{R+d\sin\theta}{R\cos\theta}\right) &= \dfrac{dx'}{R} = d\\ y' &= \left(\dfrac{d\cos\theta}{R+d\sin\theta}\right)x' - \dfrac{dR\cos\theta}{R+d\sin\theta}. \end{align*} Now the reflected ray is easily found to be $$y' = \left(\dfrac{-d\cos\theta}{R+d\sin\theta}\right)x'+\dfrac{dR\cos\theta}{R+d\sin\theta}.$$ Using the rotation formulas: \begin{align*} x' &= x\cos\theta + y\sin\theta\\ y' &= -x\sin\theta + y\cos\theta \end{align*} Substituting: \begin{align*} -x\sin\theta + y\cos\theta &= \left(\dfrac{-d\cos\theta}{R+d\sin\theta}\right)(x\cos\theta + y\sin\theta) + \dfrac{dR\cos\theta}{R+d\sin\theta}\\ (x\sin\theta - y\cos\theta)(R+d\sin\theta) &= d\cos^2\theta x + d\sin\theta\cos\theta y - dR\cos\theta\\ (R\sin\theta + d(\sin^2\theta-\cos^2\theta))x + dR\cos\theta &= (R\cos\theta + 2d\sin\theta\cos\theta)y\\ y&= \dfrac{R\sin\theta - d\cos 2\theta}{R\cos\theta + d\sin 2\theta}x+\dfrac{dR\cos\theta}{R\cos\theta + d\sin 2\theta} \end{align*} as our equation for the reflected ray. This can be confirmed:

Plot in Mathematica (excuse the half blue half orange circle!)

So the set of all possible reflected rays can be given by $$U(x,y,\theta)=(R\cos\theta + d\sin2\theta)y + (-R\sin\theta + d\cos 2\theta)x - d R \cos\theta = 0.$$ Differentiating with respect to $\theta$ and simultaneously solving the two resulting equations is unimaginably bad. So we employ a clever trick: reparametrizing $U$ by making the substitution $t = e^{i\theta}.$ Now the curve family becomes a polynomial instead of a sum of several trigonometric functions. With this new substitution we calculate \begin{align*} \cos\theta &= \dfrac{1}{2}(t+t^{-1})\\ \sin\theta &= \dfrac{i}{2}(t^{-1}-t)\\ \sin 2\theta &= \dfrac{i}{2}(t^{-2}-t^2)\\ \cos 2\theta &= \dfrac{1}{2}(t^2+t^{-2}). \end{align*} Substituting this, we find that we can reexpress $U$ as $$U(x,y,t) = (dx-idy)t^4 + (Ry-dR+iRx)t^3+(Ry-dR-iRx)t+(dx+idy)=0.$$ Both $U=0$ and $\dfrac{\partial U}{\partial t} = 0$ for some $t$ to solve for the envelope. Instead of directly solving for $t$ using some kind of quartic equation, we notice that this is only true if the polynomial $U(t)$ has a double root. Hence, $U(t)$ must have a discriminant of zero. Using the discriminant formula for a quartic polynomial, we can find that $$-4 d^6 R^6 - 60 d^6 R^4 x^2 - 12 d^4 R^6 x^2 - 192 d^6 R^2 x^4 + 312 d^4 R^4 x^4 - 12 d^2 R^6 x^4 + 256 d^6 x^6 - 192 d^4 R^2 x^6 - 60 d^2 R^4 x^6 - 4 R^6 x^6 + 24 d^5 R^6 y - 192 d^5 R^4 x^2 y + 48 d^3 R^6 x^2 y + 384 d^5 R^2 x^4 y - 192 d^3 R^4 x^4 y + 24 d R^6 x^4 y + 48 d^6 R^4 y^2 - 60 d^4 R^6 y^2 - 384 d^6 R^2 x^2 y^2 + 600 d^4 R^4 x^2 y^2 - 72 d^2 R^6 x^2 y^2 + 768 d^6 x^4 y^2 - 576 d^4 R^2 x^4 y^2 - 72 d^2 R^4 x^4 y^2 - 12 R^6 x^4 y^2 - 192 d^5 R^4 y^3 + 80 d^3 R^6 y^3 + 768 d^5 R^2 x^2 y^3 - 384 d^3 R^4 x^2 y^3 + 48 d R^6 x^2 y^3 - 192 d^6 R^2 y^4 + 288 d^4 R^4 y^4 - 60 d^2 R^6 y^4 + 768 d^6 x^2 y^4 - 576 d^4 R^2 x^2 y^4 + 36 d^2 R^4 x^2 y^4 - 12 R^6 x^2 y^4 + 384 d^5 R^2 y^5 - 192 d^3 R^4 y^5 + 24 d R^6 y^5 + 256 d^6 y^6 - 192 d^4 R^2 y^6 + 48 d^2 R^4 y^6 - 4 R^6 y^6 =0$$ is the formula for the envelope. (Also known as the catacaustic.) Here is it graphed in Mathematica:

Parameters: $R=1,\quad d=\{0.2,0.4,0.6,0.8,1.0\}$ from left to right

And this matches really well with our simulation for $d=0.8$ above.

### Next steps?

I'd like to be able to calculate the catacaustic for the set of twice-reflected light rays, or thrice-reflected light rays. I think I can do it --- need some more time?

### Sources used:

• MATLAB was used for the program and plotting, Inkscape for drawings and diagrams, Geogebra used where cited
• Mathematica was used for the heavy algebra towards the end. The notebook can be found here.
• the MATLAB densityscatter() function by Jose Manuel Amigo was used to create the "heatmap"
• My mathematical attempt was initially too complex (i.e. polar coordinates) and overlooked a much better approach at the problem. So a lot, the notion of catacaustics in particular, is inspired from these pages:
• The mathematical tricks employed when calculating the envelope were from this Wikipedia page.