load('raytrace_longrun_data_N=2e8.octavebinary')
dsize = size(ray_density_ss);
Nx = dsize(1)
Ny = dsize(2)
% Work out x and y coordinates for ray_density_ss (might be an off-by-one error here):
x = ((1:Nx)/Nx - 0.5) * (2+L);
y = ((1:Ny)/Ny - 0.5) * 2;


figure(1)
hold off
imagesc(x,y,ray_density_ss')
hold on
axis('equal')
xlabel('x')
ylabel('y')

% extract density subarray
margin=0.1
xrange = find((L/2-margin < x) & (x < L/2+margin));
yrange = find((-1+margin < y) & (y < 1-margin));
density_subarray = ray_density_ss(xrange,yrange);

figure(2)
hold off
imagesc( x(xrange), y(yrange), density_subarray' )
hold on
axis('auto')
xlabel('x')
ylabel('y')

