% Plot probability density from a ray tracing.
function density_image = density(xs, ys, L, delta)

N = numel(xs);

% Set up density density_image coordinates
Nx = 1+ceil((2+L)/delta);
Ny = 1+ceil(2/delta);
xcoord = (0:(Nx-1))*delta;
ycoord = (0:(Ny-1))*delta;
density_image = zeros(Nx,Ny);

for i=1:N-1 % ith ray connects vertex i to vertex i+1.
   xlength = xs(i+1)-xs(i);
   ylength = ys(i+1)-ys(i);
   length = sqrt( xlength.^2 + ylength.^2 );
   Npoints = ceil(3*length/delta);
   s = ((1:Npoints) - 0.5)/Npoints; % fractional arclength coordinate along ray
   xarr = xs(i) + s*xlength; % array of Npoints x-coordinates, equally spaced along ray
   yarr = ys(i) + s*ylength; % array of Npoints y-coordinates, equally spaced along ray
   xsub = ceil( (xarr + (1+L/2))/delta ); % array of x indices in density_image array
   ysub = ceil( (yarr + 1)/delta ); % array of y indices in density_image array
   density_image1Dsubs = Nx*(ysub-1) + xsub;
   density_image(density_image1Dsubs) = density_image(density_image1Dsubs) + 1;
end

hold off
imagesc(density_image')
axis equal








