%% rectangular billiard table 
teps = 10^(-10);
% set up wall parameters wW,wC where real(wW(i)'.z) = wC(i)
a = 1; b= 1;
I = sqrt(-1);
wW = [1,I,-1,-I];
wC = [a,b,a,b];

%% 
N = length(wW); 
tX =[];
tY =[];
for j = 1:N
    k = mod(j+N,N)+1;
    A = [real(wW(j)), imag(wW(j)); real(wW(k)), imag(wW(k))];
    B = [wC(j); wC(k)];
    X = A\B;
    tX = [tX,X(1)];
    tY = [tY,X(2)];
    
end
tX = [tX,tX(1)];
tY = [tY,tY(1)]; 

%% graph pool table
figure(1)
hold off
table = plot(tX,tY,'r-');
set(table,'LineWidth',3);
axis equal
hold on

%% set up billiard shot
% set up initial point
z0 = 0.8+0.3*I;
w0 = 1+2*I;
w0 = w0/abs(w0);
ball = plot(real(z0),imag(z0),'bo');
set(ball,'LineWidth',3);


%% iterate bounces
numits = 20;
z = z0;
w = w0;

for k=1:numits 
    times = (wC'-real(wW'*z))./real(wW'*w);
    mt = min(times(find(times > teps)));
    wallnum = min(find(times==mt))
    nextz = z+mt*w;
    w1 = wW(wallnum); 
    nextw = -(w1)^2/w;
    nextw = nextw/abs(nextw);
    disp 'press a key' 
    pause
    plot([real(z) real(nextz)],[imag(z) imag(nextz)],'k-');
    z = nextz;
    w = nextw;
    

end;

 %% 
hold off