Bouncing Photon (Part 2)

First we'll add two type synnonyms to make things simpler:

1
2

type Vector = (Double, Double)
type Time = Double

Now, we define a function that calculates the new position based on the velocity and the travel time:

1
2

travelT :: Vector -> Vector -> Time -> Vector
travelT (rx, ry) (vx, vy) t = (rx + t*vx, ry + t*vy)

Now, we need to calculate when and where a photon collides on the mirror given it's current position and velocity. So, suppose this collission takes time $t$ and we are currently at $(r_x,r_y)$ with position $(v_x, v_y)$ . Then the new position $(r_x + t v_x, r_y + t v_y)$ satisfies $(r_x + t v_x)^2 + (r_y + t v_y)^2 = 1 \quad \quad t > 0$

1
2
3
4
5
6
7
8

nextCol :: Vector -> Vector -> (Vector, Time)
nextCol (rx, ry) (vx, vy) = ((x,y), t)
    where
        t = let (t1, t2) =  ((- (rx*vx + ry*vy) - sqrt (vx^2 - ry^2 * vx ^ 2 + 2 * rx * ry * vx * vy + vy^2 - rx^2 * vy^2))/(vx^2 + vy^2), (- (rx*vx + ry*vy) + sqrt (vx^2 - ry^2 * vx ^ 2 + 2 * rx * ry * vx * vy + vy^2 - rx^2 * vy^2))/(vx^2 + vy^2)) in
            if t1 > 0 && t2 > 0   then min t1 t2 else
            if t1 > 0            then t1
            else             t2
        (x,y) = travelT (rx, ry) (vx, vy) t

Now, given the velocity vector of a photon, we would like to find the new velocity vector after it bounces. To do that, we extract the component of the velocity in the direction of the normal to the surface at that point and reverse it. And it so happens that the position vector on the unit circle is an unit vector in the opposite direction to that of the normal.

1
2
3
4
5

reflectV :: Vector -> Vector -> Vector
reflectV (ux, uy) (x, y) = (ux - 2*cx, uy - 2*cy) --assuming (x, y) is a unit vector, which it is
    where
        comp = ux * x + uy * y
        (cx, cy) = (comp*x, comp*y)

Now, we will add a function that takes the current position, velocity and time and updates it to that until the next reflection.

1
2
3
4
5
6

update :: (Vector, Vector, Time) -> (Vector, Vector, Time)
update (r, v, t) = (r', v', t')
    where
        (r', t'')   = nextCol r v
        t'          = t + t'' --time is additive, unlike position
        v'          = reflectV v r'

And now, we are ready, we update the state repeatedly until the photon escapes. When it does, we read off the time.

1
2
3
4
5

ans = head $ dropWhile f (iterate update initial)
    where
        initial = ((-0.6,0.3),(0.5,(sqrt 3)/2),0)
        f ((x,y), _, _) = not $ let tanTheta = (y/x) in 
            0 < x && -0.0874887 < tanTheta && tanTheta < 0.0874887 --is in the right half of the plane and is within 5 degrees