Matlab/Photons/Monte Carlo

[Stiganator]

I'm writing a matlab program for a photon simulation.
I'm supposed to yield R=.3616, T=.3565, but I end up with R~0.4 and T~0.28.

If anyone has any input it would be appreciated.

Does anyone know how this method would compare to using for loops? Running this with 10000 photons takes ~6 secs per trial, but I can't imagine using for loops would be any faster.









==================================================
format loose;
whitebg('w');
clc;
clear all; close all;

Results=zeros(20,4);

for(j=1:20)
%INPUTS
%================================================
%NumTrials = 5;%input('Number of trials= ');
NumPhotons = 1000;%input('Number of Photons= ');
d = 0.02;%input('Tissue Thickness (in cm)= ');
%deltaz = .025;%input('Step Size (in cm)= ');
musuba = 10;%input('Absorption Coefficient (cm^-1)= ');
musubs = 90;
g = 0;
%n1 = 1;
%n2 = 1;


%if(mod(NumTrials,3)==0)
% PlotRows = (NumTrials/3);
%else
% PlotRows = (((NumTrials/3)-(mod(NumTrials,3)/3))+1);
%end
%Fig1=figure();


%INITIALIZE
tic;
A=0;
T=0;
R=0;
Z=zeros(1,NumPhotons);
Cz=ones(1,NumPhotons);


while(NumPhotons>0)
ds=-log(rand(1,NumPhotons))/(musuba+musubs);
Z=Z+(Cz.*ds);



ind=find(rand(1,(length(ds)))<(musuba*ds));%Absorbed
A=A+length(ind);
NumPhotons=NumPhotons-length(ind);
for(k=length(ind):-1:1)
Z(ind(k))=[];
Cz(ind(k))=[];
ds(ind(k))=[];
end
ind=find(Z<0);%Reflected out
R=R+length(ind);
NumPhotons=NumPhotons-length(ind);
for(k=length(ind):-1:1)
Z(ind(k))=[];
Cz(ind(k))=[];
ds(ind(k))=[];
end

ind=find(Z>d);%Transmitted
T=T+length(ind);
NumPhotons=NumPhotons-length(ind);
for(k=length(ind):-1:1)
Z(ind(k))=[];
Cz(ind(k))=[];
ds(ind(k))=[];

end



%====================
n_rn = length(Cz);
if (g==0)
mu = 2*rand(1,n_rn) - 1;
else
mu = (1+g.^2 - ( (1-g.^2)./ (1-g+ 2*g.*rand(1,n_rn)) ).^2) ./ (2*g);
end
theta=acos(mu);
phi = 2*pi*rand(1,n_rn);
CzTemp=Cz;
i = find( abs(Cz) > 0.99999 );
CzTemp(i) = Cz(i).*cos(theta(i))./abs(Cz(i));

i = find( abs(Cz) <= 0.99999 );
CzTemp(i) = -sin(theta(i)) .* cos(phi(i)) .* sqrt(1-Cz(i).^2)+ Cz(i).*cos(theta(i));

Cz = CzTemp;
%========================
end
Time=toc;
Results(j,1)=R;
Results(j,2)=T;
Results(j,3)=A;
Results(j,4)=Time;
end
mean(Results,1))/1000
mean(Results,2))/1000
mean(Results,3))/1000
sum(Results,4))

 

[CycloWizard]

The way you've done it is definitely faster than looping. As for the disparity in results, in Monte Carlo simulations, the numbers you're getting aren't that far from the expected values. If you want to improve the accuracy, the rule of thumb is to simply run more simulations (i.e. increase NumPhotons) to as high a value as you can within your computational time constraint. As the number of samples becomes infinite, your answers will approach your expected answers as long as your code is right. 1000-10000 samples is really pretty low. I would probably try 10^6 or so and see if the results are any better. The time it takes will only be 600 seconds=10 minutes, so just let it run once and see if the answers are "better." If that doesn't work, I can take a look at your code tomorrow and see if I can find something wrong.