%Code from Leathers et al., 2004, NRL report
% PSICAM Analysis
% Matlab code for computing the probability of photon survival from wall to wall within a PSICAM



a = 0.2; %absorption coeff
b = 50; %scattering coeff
model = 1; % 1=isotropic, 2=HG, 3=Petzold
g = 0.92; %scattering asymmetry factor (required for HG only)
r = 0.05; % radius of sphere (m)
N=1000; %number of photons traced
minweight=0.0001; % min photon weight before discarding it
maxns=10; %maximum number of scatters for each photon
c = a+b; % attenuation coefficient
w0 = b/c; %single-scattering albedo
count=0;
for i=1:N, % trace photons one at a time
weight = 1; % new photon, generated at position (0,0,-r): bottom of sphere
ns = 0;
y = 0;
z = -r;
x = 0;
val=(i-0.5)/N; %step uniformly over interval (0,1)
mu=sqrt(1.0-val); % diffuse light directions
gamma = mu; % alpha, beta, gamma = direction cosines
beta = 0; %move in x-z plane
alpha = sqrt(1.0-mu^2); %alpha^2+beta^2+gamma^2=1
while (ns <= maxns) & (weight > minweight) % keep tracing photon path
%how far does it go?
l=-(1/c)*log(1-rand); %total path length
z = z+gamma*l; %new position
x = x+alpha*l;
y = y+beta*l;
if ((x^2+y^2+z^2) > r^2); %does it reach the wall?
count=count + weight;
break %if the photon hits the wall, break out of while loop and don't scatter it.
end
% scatter photons
weight=weight*w0; % fraction scattered (versus absorbed) given by w0
ns = ns+1; %number of times the photon has been scattered
gammas = scat(model,g); % compute polar scattering angle
phis = 2*pi*(rand);
[alpha,beta,gamma] = chgdir(alpha,beta,gamma,gammas,phis);
end
end
F=count/N;
disp(sprintf('The probability of photon survival from wall to wall is %g',F))