The following solution was supplied by Brian Medley
Hi Marijke, Since you expressed curiosity about my photon program, it is attached here. I am NOT proud of the programming; don't ask why all the main code appears twice - that's just the way it turned out. In the program: (dx,dy) is a direction vector for the photon path (Oi,Oj) is the centre of the last mirror (Ni,Nj) is the centre of the (possibly) next mirror L is the value of a parameter in the vector equation of the photon's current path that generates the postion of closest approach of that (projected) path to (Ni,Nj) (Cx,Cy) and (Px,Py) give the projected position of closest approach relative to (Ni,Nj) and the origin, repectively (Ax,Ay) and (Newax, Neway) are the coordinates of impacts with the mirrors Dot is the value of the scalar product of vectors used to find the resolved part of (dx,dy) parallel to the normal to the mirror - it's used in finding the rebound direction.
2 point 49
20 Dist=0
50 screen 823
60 for I=-2 to 2
61 for J=-2 to 2
62 circle (300+111*I,300-111*J),111/3,1
63 next J
64 next I
100 Dx=1:Dy=0
110 Ax=0.5:Ay=0.1
120 Oi=0:Oj=0
130 Ni=1:Nj=0
140 L=(Dx*(Ni-Ax)+Dy*(Nj-Ay))/(Dx*Dx+Dy*Dy)
150 Cx=Ax-Ni+L*Dx
160 Cy=Ay-Nj+L*Dy
161 Px=Ni+Cx
162 Py=Nj+Cy
170 D=Cx*Cx+Cy*Cy
180 if D<1/9 then 'print "hit";Ni,Nj
190 Backup=sqrt(1/9-D)
200 D=sqrt(Dx*Dx+Dy*Dy):'print "if ";D;"<>1 then we are in trouble"
210 Newax=Px-Dx*Backup/D
220 Neway=Py-Dy*Backup/D
225 line (300+Ax*111,300-Ay*111)-(300+Newax*111,300-Neway*111),4
230 'print "contact at";Newax,Neway
240 Rx=(Newax-Ni):Ry=(Neway-Nj)
241 'print Rx,Cx,Ry,Cy
250 Dot=Dx*Rx+Dy*Ry:'print "if";Dot;"is positive we're in trouble"
260 Newdx=Dx-Rx*Dot*18
270 Newdy=Dy-Ry*Dot*18
280 'print "zooming off in direction";Newdx;Newdy
290 Dirlength=sqrt(Newdx*Newdx+Newdy*Newdy)
291 'print "if";Dirlength;"<>1 were in trouble"
292 Newdx=Newdx/Dirlength
293 Newdy=Newdy/Dirlength
300 Dist=sqrt((Newax-Ax)^2+(Neway-Ay)^2)
310 'print "distance so far=";Dist
400 Dx=Newdx:Ax=Newax
410 Dy=Newdy:Ay=Neway
412 circle (300+111*Ax,300-111*Ay),3,7
415 Oi=Ni:Oj=Nj
420 Xd=sgn(Dx):Yd=sgn(Dy)
430 K=1
440 repeat
450 Ni=Oi+(K+1)*Xd:Nj=Oj-Yd
460 repeat
465 Ni=Ni-Xd:Nj=Nj+Yd
467 'print "let's try";Ni,Nj
470 L=(Dx*(Ni-Ax)+Dy*(Nj-Ay))/(Dx*Dx+Dy*Dy)
480 Cx=Ax-Ni+L*Dx
490 Cy=Ay-Nj+L*Dy
510 Px=Ni+Cx
520 Py=Nj+Cy
530 D=Cx*Cx+Cy*Cy
535 'print "approach distance";sqrt(D)
536 'input A
540 'if D<1/9 then print "hit";Ni,Nj else if D<0.12 then print "shave at";D
550 until D<1/9 or Ni=Oi
560 K=K+1
570 until D<1/9
600 Backup=sqrt(1/9-D)
610 D=sqrt(Dx*Dx+Dy*Dy):'print "if ";D;"<>1 then we are in trouble"
620 Newax=Px-Dx*Backup/D
630 Neway=Py-Dy*Backup/D
635 line (300+Ax*111,300-Ay*111)-(300+Newax*111,300-Neway*111),4
640 'print "contact at";Newax,Neway
650 Rx=(Newax-Ni):Ry=(Neway-Nj)
660 'print Rx,Cx,Ry,Cy
670 Dot=Dx*Rx+Dy*Ry:'print "if";Dot;"is positive we're in trouble"
680 Newdx=Dx-Rx*Dot*18
690 Newdy=Dy-Ry*Dot*18
700 'print "zooming off in direction";Newdx;Newdy
800 Dirlength=sqrt(Newdx*Newdx+Newdy*Newdy)
810 'print "if";Dirlength;"<>1 were in trouble"
820 Newdx=Newdx/Dirlength
830 Newdy=Newdy/Dirlength
835 ExtraDist=sqrt((Newax-Ax)^2+(Neway-Ay)^2)
840 Dist=Dist+ExtraDist
850 'print "distance so far=";Dist
855 if Dist>10 then goto 1000
860 'input A
870 goto 400
1000 Overshoot=Dist-10
1010 Prop=Overshoot/ExtraDist
1020 Finalx=Ax+(1-Prop)*(Newax-Ax)
1030 Finaly=Ay+(1-Prop)*(Neway-Ay)
1035 circle (300+111*Finalx,300-111*Finaly),3,6
1040 'print "at t=10, co-ords are";Finalx,Finaly
1050 print "final dist from O is";sqrt(Finalx^2+Finaly^2)