% What is the Lyapunov exponent for Stadium billiards?

N = 1000
L = 5
x0 = L*(rand-0.5)
y0 = L*(rand-0.5)
phi0 = 2*pi*rand
epsilon = eps*10


[x1,y1] = raytrace(N, x0, y0, phi0, L);

x0 += randn*epsilon;
y0 += randn*epsilon;
phi0 += randn*epsilon;

[x2,y2] = raytrace(N, x0, y0, phi0, L);

% Calculate the curve showing divergence of trajectory (x2,y2) from (x1,y1):
distances = sqrt( (x2-x1).^2 + (y2-y1).^2 );
steps = ( 1:length(x1) ) -1;

% Find where things start to come apart.
for k=1:N
   if (distances(k) > 0.1)
      break % get out of the loop
   end
end

figure(1)
hold off
plot(x0,y0,'ro')
hold on
plot([x0,x1(1:k)],[y0,y1(1:k)],'b-')
plot([x0,x2(1:k)],[y0,y2(1:k)],'g--')
axis('equal')
xlabel('x')
ylabel('y')
legend('Starting Point','Trajectory 1','Trajectory 2')
print('lyapunov_trajectories.pdf','-dpdfwrite')


% Where shall I end the separation plot?
k2 = min(k+20,N)

figure(2)
hold off
plot(steps(1:k2),log(distances(1:k2)),'b')
hold on
plot(steps(1:k2), log(distances(1)) + steps(1:k2), 'r--')
xlabel('step')
ylabel('log_e distance')
legend('Trajectory separation','Lyapunov exponent = 1')
print('lyapunov_separation.pdf','-dpdfwrite')