% Plot probability density from a ray tracing.
% Attempt to improve efficiency of density.m.
function density_image = density2(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);

max_length = 1.25*(2+L); % Biggest chord we can draw, with some margin
Npoints_max = max_length / delta; % Biggest size s array we should need, with some margin
s = ((1:Npoints_max) - 0.5); % dimensionless arclength array
xarr = zeros(1,Npoints_max);
yarr = zeros(1,Npoints_max);
xsub = zeros(1,Npoints_max);
ysub = zeros(1,Npoints_max);
density_image1Dsubs = zeros(1,Npoints_max);
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.*xlength + ylength.*ylength );
   Npoints = ceil(length./delta);
   xarr(1:Npoints) = xs(i) + s(1:Npoints).*(xlength/Npoints); % array of Npoints x-coordinates, equally spaced along ray
   yarr(1:Npoints) = ys(i) + s(1:Npoints).*(ylength/Npoints); % array of Npoints y-coordinates, equally spaced along ray
   xsub(1:Npoints) = ceil( (xarr(1:Npoints) + (1+L/2))./delta ); % array of x indices in density_image array
   ysub(1:Npoints) = ceil( (yarr(1:Npoints) + 1)./delta ); % array of y indices in density_image array
   density_image1Dsubs(1:Npoints) = Nx*(ysub(1:Npoints)-1) + xsub(1:Npoints);
   density_image(density_image1Dsubs(1:Npoints)) ++;
end

hold off
imagesc(density_image')
axis equal








