A photon moving at speed 1 in the plane starts at at heading due east. Around every integer lattice point in the plane, a circular mirror of radius has been erected. How far from the origin is the photon at ?
Solution:
You have to be careful to use enough digits here, because there are 13 reflections, and each reflection from a curved surface amplifies the error. We needed 11 additional digits.
The state of the photon is where is the position, the velocity, the distance travelled. We always have .
Starting in this state, you hit a mirror centred at if
(assuming that mirror is in the forward direction, and you don't hit anything else first).
If so, let be the point where it hits. This is where
. Then the next state is
where .
> | plots[display]([plottools[circle]([0,0],1), plots[arrow]([2,-2],[-1,2],shape=arrow,colour=blue), plots[arrow]([1,0],[1,2],shape=arrow,colour=blue), plots[arrow]([1,0],[2,0],shape=arrow,colour=red) ],axes=none,scaling=constrained); |
Global variables:
x, y, dx, dy, s
for current state;
if going mainly in
direction,
iy = 0
and
ix
is integer to go in
direction to next mirror;
if going mainly in direction, ix = 0 and iy is integer to go in direction to next mirror.
> | move:= proc() # return true if you hit the mirror or end of path # and update state in this case, # false if you don't (leaving state unchanged) local mx, my, t, px, py,h,ds; global x,y,dx,dy,s,ix,iy; if iy = 0 then # going in x direction mx:= round(x+ix-.001); my:= round(y+dy/dx*ix); ds:= sqrt(1+(dy/dx)^2)*abs(mx-x); else # going in y direction my:= round(y+iy); mx:= round(x+dx/dy*iy); ds:= sqrt(1+(dx/dy)^2)*abs(my-y); fi; if abs((mx-x)*dy - (my-y)*dx) < 1/3 then t:= min(fsolve((x+t*dx-mx)^2 + (y+t*dy-my)^2=1/9,t)); px:= x+t*dx; py:= y+t*dy; ds:= sqrt((px-x)^2+(py-y)^2); if s + ds >= 10 then x:= x + (10-s)*dx; y:= y + (10-s)*dy; s:= 10; return true; else h:= 18*((mx-px)*dx+(my-py)*dy); s:= s + ds; dx:= dx - h*(mx-px); dy:= dy - h*(my-py); # dx^2 + dy^2 should stay constant, but make sure h:= sqrt(dx^2+dy^2); dx:= dx/h; dy:= dy/h; x:= px; y:= py; return true; fi # no hit, but maybe we'd go too far elif s + ds >= 10 then x:= x + (10-s)*dx; y:= y + (10-s)*dy; s:= 10; return true; fi; false end; |
Initial conditions.
> | x,y,dx,dy,s:= 0.5, 0.1, 1, 0, 0; path:= [x,y]; ix:= 1; iy:= 0; _Envsignum0:= 0; Digits:= 50; |
> | while s < 10 do if move() then # reassess directions if abs(dx) >= abs(dy) then ix:= signum(dx); iy:= 0 else ix:= 0; iy:= signum(dy); fi; path:= path,[x,y]; else # go to next mirror in same direction ix:= ix+signum(ix); iy:= iy+signum(iy); fi od: |
The final state:
> | x,y,dx,dy,s; |
The answer:
> | ans:=sqrt(x^2+y^2); |
> | with(plots): display({plot([path]),seq(seq(plottools[circle]([i,j],1/3),i=-1..1),j=-1..2)},scaling=constrained,axes=frame); |
Warning, the name changecoords has been redefined
Just for comparison, here it is with 40 digits.
> | x,y,dx,dy,s:= 0.5, 0.1, 1, 0, 0: path:= [x,y]: ix:= 1: iy:= 0: _Envsignum0:= 0: Digits := 40: while s < 10 do if move() then # reassess directions if abs(dx) >= abs(dy) then ix:= signum(dx); iy:= 0 else ix:= 0; iy:= signum(dy); fi; path:= path,[x,y]; else # go to next mirror in same direction ix:= ix+signum(ix); iy:= iy+signum(iy); fi od: ans40:= sqrt(x^2+y^2); |
> | ans40 - ans; |
And here's what we would have got with a starting point different by in the y direction.
> | path0:= path: x,y,dx,dy,s:= 0.5, 0.1+1e-10, 1, 0, 0: path:= [x,y]: ix:= 1: iy:= 0: _Envsignum0:= 0: Digits := 40: while s < 10 do if move() then # reassess directions if abs(dx) >= abs(dy) then ix:= signum(dx); iy:= 0 else ix:= 0; iy:= signum(dy); fi; path:= path,[x,y]; else # go to next mirror in same direction ix:= ix+signum(ix); iy:= iy+signum(iy); fi od: |
> | display({plot([path],colour=blue),plot([path0],colour=red),seq(seq(plottools[circle]([i,j],1/3),i=-1..1),j=-1..2)},scaling=constrained,axes=frame); |
> | seq(path[i]-path0[i],i=1..nops([path])); |
> | plots[logplot]([seq([i,sqrt(%[i][1]^2+%[i][2]^2)],i=1..nops([path]))],style=point,symbol=circle); |
Back to Solutions to the SIAM $100 Challenge