function [xs,ys] = raytrace(N, x0, y0, phi0, L)
% N = number of intersections to find
% x0, y0 = initial position
% phi0 = initial direction (ccw from +x)
% L = length of straight side (the circular ends have unit radius).
%
% Variables for current working position:
y1 = y0;
x1 = x0;
phi = phi0;
t0 = cputime;

xs = zeros(1,N);
ys = zeros(1,N);
for i=1:N 
   % (0) Validate initial position (inside stadium?)
   
   
   % (1) Find next vertex
   % Approach: follow the parametric ray paths for positive t:
   %  x = x1 + t*cos(phi)
   %  y = y1 + t*sin(phi)
   % Intersection with top or bottom side (y = +/- 1):
   sinphi = sin(phi);
   cosphi = cos(phi);
   ttop = (1-y1)/sinphi;
   tbot = (-1-y1)/sinphi;
   tside = max([ttop, tbot]);
   xside = x1 + tside*cosphi;
   if (xside < -L/2)
      % Intersection with left circular end
      quadratic_left = [1, L*cosphi + 2*x1*cosphi + 2*y1*sinphi, (L^2)/4 - 1 + L*x1 + x1^2 + y1^2];
      tleft = roots(quadratic_left);
      t = max(tleft); % should be 2 real roots; the most positive one is wanted.
      x1 = x1 + t*cosphi;
      y1 = y1 + t*sinphi;
      phi_normal = atan2(y1, x1 + L/2) + pi;
   elseif (xside > L/2)
      % Intersection with right circular end
      quadratic_right = [1, -L*cosphi + 2*x1*cosphi + 2*y1*sinphi, (L^2)/4 - 1 - L*x1 + x1^2 + y1^2];
      tright = roots(quadratic_right);
      t = max(tright); % should be 2 real roots; the most positive one is wanted.
      x1 = x1 + t*cosphi;
      y1 = y1 + t*sinphi;
      phi_normal = atan2(y1, x1 - L/2) + pi;
   else % The remaining possibility is just a side intersection.
      yside = y1 + tside*sinphi;
      x1 = xside;
      y1 = yside;
      phi_normal = pi * (1 + sign(y1)/2);
   end   
   
   % (2) Determine new value of phi. For neatness, keep in the range 0..2pi.
   phi = mod(2 * phi_normal - phi - pi, 2*pi);
   % (3) Save current vertex
   xs(i) = x1;
   ys(i) = y1;

   % This part is for monitoring long runs. Typically spits out an update every 250 sec or so.
   if (i/1000000 == round(i/1000000))
      current_time = ctime(time)
      completed_iteration = i
      out_of = N
      cpu_time_sec = cputime - t0
      x1
      y1
   end
end