function [pt,d]=project1(r)

hold off;
% Initialization
p0=[0.6;0.1];
v=[1;0];
T=10;
r=1/3;

% Computation 
[ptmp,vtmp,t]=nextpoint(p0,v,T,r);

hold on;
axis equal;
iter=0;
while t<T,
    [pt,vt,dt]=nextpoint(ptmp,vtmp,T-t,r);
    t=t+dt;
    ptmp=pt;
    vtmp=vt;
    iter=iter+1;
end
% iter is the number of reflections
iter
% d is the distance we are looking for.
% only the first four digits are right. About 12 digits are lost along the
% computation. This is a very ill-conditioned problem. To get more digits, 
% we will have to use a software (e.g., Maple or Mathematica) 
% which can provide greater accuracy.
d=norm(pt);


function [pt,vt,t]=nextpoint(p0,v,T,r)
% starting from p0 with the velocity v, find the point pt where the photon
% hits the circle and new velocity vt after the reflection.
% march along the straight line until it intersects with one of the
% circles.
% Warning: This method may miss those mirrors on which the photon will hit
% upon almost tangentially. But it works fine for this particular initial
% setting.
dt=0.1;
ptmp=p0+dt*v;
p0c=round(ptmp);
rtmp=norm(ptmp-p0c);
t=dt;
while (rtmp > r) && (t < T),
    t=t+dt;
    ptmp=ptmp+dt*v;
    p0c=round(ptmp);
    rtmp=norm(ptmp-p0c);
end

% once we have found the circle where the photon is going to hit upon, we can
% find the point where the photon hits the circle by solving a quadratic
% equation about the hitting time.
if rtmp<r
    pd=ptmp-dt*v-p0c;
    a=norm(v)^2;
    b=2*pd'*v;
    c=norm(pd)^2-r^2;
    dis=sqrt(b^2-4*a*c);
    % b is always negative, so bring the numerator down to denominator to
    % avoid cancellation error.
    ddt=2*c/(dis-b);
    t=t-dt+ddt;
end

if (t>T)
    pt=p0+T*v;
    vt=v;
    t=T;
else
    pt=p0+t*v;
    vt=nextvelocity(v,pt,p0c,r);
end

% now plot the path and the circle
h=0.01;
tp=0:h:t;
x1=p0(1)+v(1)*tp;
y1=p0(2)+v(2)*tp;

theta=0:h:2*pi;
x2=p0c(1)+r*cos(theta);
y2=p0c(2)+r*sin(theta);

plot(x1,y1,'r-',x2,y2,'b-')

function vt=nextvelocity(v,p,p0c,r)
% given the incident velocity, the reflecting point, and the center of the
% circle, compute the reflected velocity by Householder reflector.
a=p-p0c;
vt=v-2*(a'*v)/(a'*a)*a;