%Code from Leathers et al., 2004, NRL report
% 3-D Diffusion From Point Source
% Monte Carlo Matlab script for computing how light spreads from a point
% source.


a = 1; % absorption coefficient (1/m)
c = 6; % attenuation coefficient (1/m)
w0 = (c-a)/c; % single-scattering albedo
P0 = 1; % source intensity (W)
N = 10e4; % number of photons to trace
ns = 30; % number of scatters to trace
Nd = 20; % number of detectors
dd = 0.02; % detector spacing (m)
A = zeros(1,Nd); % initalize detectors (all photons)
for i=1:N
p=[0,0,0]; % initial position (x,y,z)
muz = 1-2*rand; % isotropic source
phi = 2*pi*rand; % isotropic source
mux = sqrt(1-muz^2)*cos(phi); % initial x-direction cosine
muy = sqrt(1-muz^2)*sin(phi); % intitial y-direction cosine
r=0; % intial distance from source
w = P0; % initial photon weight
for j=1:ns
s = -log(rand)/c; % geometric path length
p = p + s*[mux,muy,muz]; % mo
r= norm(p); % radial position
index = ceil(r/dd); % find detector number
if (index<=Nd), % is photon within detector grid?
A(index) = A(index) + (1-w0)*w; % count absorption
end
w = w*w0;
mu= 1 - 2*rand; % isotopic scattering
phis = 2*pi*rand; % isotopic scattering
[mux,muy,muz]=chgdir(mux,muy,muz,mu,phis); % new direction
end
end
r = [1:Nd]*dd; % radii for detectors
A = A/N; % normalize detector response
A = A./(dd*4*pi*r.^2); % convert response to W/m^3
ht = 3*a*c*P0*exp(-r*sqrt(3*a*c))./(4*pi*r); % theory
semilogy(r,ht,r,A,'*')
grid on
xlabel('r'), ylabel('W/m^3')
legend('theory','MC (total)')