% This is a PROTOTYPE Monte-Carlo simulation code
% for photon transport through a homogeneous slab.
% The scattering statistics in particular are wrong.

clear all
N = 20000;              % number of photon groups to track
depth = 5;              % thickness of slab
plane_loc = 5.5;        % location of imaging plane
mu = 1;                 % mean distance before scattering
alpha = 0.1;            % probability of getting absorbed
sigma = 0.5;            % variance of distance before scattering
xpos = zeros(N,3);      % initial position of photons (all at the origin)
dir = [zeros(N,1) zeros(N,1) ones(N,1)];        % initial direction of photons (up)
intensity = ones(N,1);  % initial proportion of photons in each group
hit_list_x = [];        % The hit list keeps track of where photons hit imaging plane
hit_list_y = [];
hit_intensity = [];

iter = 1;
while (iter < 1000) & (length(xpos) > 0),   % length(xpos) = 0 means no photons are left

  [m,n] = size(xpos);
  dir = dir + 0.1*randn(m,n);                                    % compute new scattering directions
  dir_len = max(sum(dir.^2,2).^0.5, eps*ones(m,1));              % length of each new direction vector
  dir = [dir(:,1)./dir_len dir(:,2)./dir_len dir(:,3)./dir_len]; % all directions have length one
  sdist = sqrt(sigma)*randn(m,1) + mu;                           % new scattering distances
  step = [dir(:,1).*sdist dir(:,2).*sdist dir(:,3).*sdist];

  xpos_new = xpos + step;                                        % new position    

  active_index = find(xpos_new(:,3) <= depth);                   % Identify which photons are still in
  escaped_index = find(xpos_new(:,3) > depth);                   %   the medium and which ones escaped.
  if length(escaped_index) > 0                                   % For photons that escaped, figure out
    escaped = xpos_new(escaped_index,:);                         %   where they hit imaging plane and
    esc_dir = dir(escaped_index,:);                              %   tally up the hits.
    t = (plane_loc - escaped(:,3))./esc_dir(:,3);
    hit_loc = escaped(:,1:2) + [t.*esc_dir(:,1) t.*esc_dir(:,2)];
    hit_list_x = [hit_list_x; hit_loc(:,1)];
    hit_list_y = [hit_list_y; hit_loc(:,2)];
    hit_intensity = [hit_intensity; intensity(escaped_index)];
  end
  xpos = xpos_new(active_index,:);                               % Only photons that didn't escape make
  dir = dir(active_index,:);                                     %   it to the next iteration.
  intensity = (1-alpha)*intensity(active_index);                 % Assume alpha percent of photons absorbed.
  iter = iter + 1;
end
plot(hit_list_x,hit_list_y,'b+'); axis('equal'); drawnow; pause

xmax = 4; ymax = 4;                                            % Make a 2d histogram
nbins = 64;                                                 
hit_list_x = max(min(hit_list_x,xmax),-xmax);
hit_list_y = max(min(hit_list_y,ymax),-ymax);
xlist = round((nbins-1)*(hit_list_x+xmax)/(2*xmax))+1;
ylist = round((nbins-1)*(hit_list_y+ymax)/(2*ymax))+1;
iimage = zeros(nbins,nbins);
for j = 1:length(hit_intensity)                                % There must be a way to do this without a loop!
  iimage(xlist(j),ylist(j)) = iimage(xlist(j),ylist(j))+hit_intensity(j);
end
imagesc(iimage); axis('image')