rng('shuffle'); %use current time as seed
%Define constants
mu_a=0.1; %absorption coefficient [1/cm]
mu_s=100; %scattering coefficient [1/cm]
mu_t=mu_a+mu_s; %total interaction attenuation coefficient [1/cm]
g=0.9;%expected value of cos theta
N=1000;% number of photons
z1=201;%number of grid elements in z
r1=201;%number of grid elements in r
A=zeros(r1,z1); %create energy deposition array
dz=0.005;%grid line seperation in z
dr=0.005;%grid line speration in r
ni=1;%refractive index of air
nt=1.37;%refractive index of tissue
m=1000; %1 chance in m of photon packet surviving
Rsp=(ni-nt)^2/(ni+nt)^2;%Specular reflectance
for j=1:N
%Initialize
x=0;
y=0;
z=0;
mu_x=0;
mu_y=0;
mu_z=1;
W=1;
%inital loss of W due to specular reflectance
W=W-Rsp;
s=0;
dead=0;%photon is initally not dead
while dead==0;
if s==0;
%Launch photon
s=-log(rand)/mu_t;
end
%Find distance to boundary db
%Check if photon will hit boundary
%Compute distance to boundary
if mu_z==0;
db=inf;
elseif mu_z<0;
db=(0-z)/mu_z;
elseif mu_z>0;
db=inf; %(If layer is finite change condition)
end
%Decide if step size is greater than db
if db<=s
%Move photon to boundary
s=s-db;
x=x+mu_x*db;
y=y+mu_y*db;
z=z+mu_z*db;
ai=acos(abs(mu_z));
ac=asin(ni/nt);
if ai>ac;
r_ai=1;
else
at=asin((ni*sin(ai))/nt);
r_ai=1/2*(((sin(ai-at)^2)/(sin(ai+at)^2))+((tan(ai-at)^2)/(tan(ai+at)^2)));
end
%Compute probabilty of a photon packet being internally reflected
%or transmitted
if rand<=r_ai
mu_z=-1*mu_z;
else
%Transmits;
dead=1;
end
else
%For photon that does not hit boundary
x=x+(s*mu_x);
y=y+(s*mu_y);
z=z+(s*mu_z);
s=0;
%Find new direction cosines
theta=acos((1/(2*g))*(1+(g^2)-((1-g^2)/(1-g+2*g*rand))^2));
psi=2*pi*rand;
if abs(mu_z)>0.99999
mu_x2=sin(theta)*cos(psi);
mu_y2=sin(theta)*sin(psi);
mu_z2=sign(mu_z)*cos(theta);
else
mu_x2=(sin(theta)*(mu_x*mu_z*cos(psi)-mu_y*sin(psi))/sqrt(1-mu_z^2))+mu_x*cos(theta);
mu_y2=(sin(theta)*(mu_y*mu_z*cos(psi)+mu_x*sin(psi))/sqrt(1-mu_z^2))+mu_y*cos(theta);
mu_z2=-sin(theta)*cos(psi)*(sqrt(1-mu_z^2))+mu_z*cos(theta);
end
mu_x=mu_x2;
mu_y=mu_y2;
mu_z=mu_z2;
%Find indices of current position
r=sqrt(x^2+y^2);
iz=z/dz;%index of z
ir=r/dr;%index of r
izup=ceil(iz);%upper z grid
izdown=floor(iz);%lower z grid
irup=ceil(ir);%upper r grid
irdown=floor(ir);%lower r grid
%compute distances to grids
dW=(mu_a/mu_t)*W;
%Check if index goes beyond boundary
if irup+1>r1||irdown+1>r1||izup+1>z1||izdown+1>z1
%Skip;
elseif z>2;
dead=1;
else
%Index is within boundary
%Score weight
dist1=abs(z-izup*dz)/dz;
dist2=abs(z-izdown*dz)/dz;
dist3=abs(r-irup*dr)/dr;
dist4=abs(r-irdown*dr)/dr;
A(irup+1,izup+1)= A(irup+1,izup+1)+dist1*dist3*dW;
A(irup+1,izdown+1)=A(irup+1,izdown+1)+dist2*dist3*dW;
A(irdown+1,izup+1)=A(irdown+1,izup+1)+dist4*dist1*dW;
A(irdown+1,izdown+1)=A(irdown+1,izdown+1)+dist4*dist2*dW;
%Reduce weight
W=W-dW;
end
%Roulette
if W<=0.0001
if rand <= 1/m
W=m*W;
else
W=0;
dead=1;
end
end
end
end
end
figure
fluence=(sum(A,2)/(dz*N*mu_a));
fluence(201,:)=[];
fluence(1,:)=[];
semilogy(.005:.005:.9950,fluence)
ylim([0.5,10]);
xlabel('z [cm]')
ylabel('Fluence[-]')
title('100,000 photon packages')
% saveas(gcf,'1000photons-1.png')
