蒙特卡洛源程序

clear;
clc;
N=10e3;

prob_of_survival = 0.75;       % b/c
c=0.303;
g=0.924;
receiver_z=30;

minpower=1e-3;          % minimum power value for photons before they are terminated by rouletting

beamWaist=0.001;%波束宽度
randVals = rand（N1）;
diverg=0.00125;

radius = beamWaist*sqrt（-log（1-randVals））;     % 半径采样值
invF = -diverg/beamWaist;


phiAng = （2*pi）.* rand（N1）;   % 方位角初始  phi0 
divAng = -invF.*radius;        % 散射角theta0


cosPhi = cos（phiAng）;
sinPhi = sin（phiAng）;


rec_loc = zeros（N4）;     % location of the received photon on the zy rec. plane
total_rec_packets = 0;              % Total number of packets to cross detector
total_rec_power = 0;                % Total power to cross detector
rec_dist = zeros（N1）;    % total distance photon traveld at point of reception
totaldist = zeros（N1）;   % record total distance traveled by each photon
weight = ones（N1）;  % set weights to one
max_uz=0.999;



%初始传播矢量
ux=sin（divAng）.*cos（phiAng）;
uy=sin（divAng）.*sin（phiAng）;
uz=cos（divAng）;


x = radius.*cosPhi;
y = radius.*sinPhi;
z= zeros（N1）;        
%hist3（[x y][100 100]）
photonre=N;
num_photons=N;  
photon = zeros（num_photons8）;
photon（:8） = ones（num_photons1）;  % set all active

tic;

while photonre>0
    rand1=rand（N1）;
    phi =rand（N1）.*2.*pi;   % 方位角   
for i=1:N
    
    if （photon（i8） == 1）%未被接收
         dn=（-1/c）*（log（rand1（i）））;%散射前后光子位置的几何距离  步长
         theta=acos（（1/g*2）*（（1+g^2）-（（1-g^2）/（1-g+（2*g.*rand（N1）））.^2）））;   %散射角
   
         xstep=dn*ux（i）;
         ystep=dn*uy（i）;
         zstep=dn*uz（i）;
            if（（z+zstep）>receiver_z）
                   if（uz~=0）
                    % z距离接收面距离
                    z_dist_rec_intersection = receiver_z -z;                   
                    y_dist_rec_intersection = z_dist_rec_intersection*uy（i）/uz（i）;                   
                    x_dist_rec_intersection = z_dist_rec_intersection*ux（i）/uz（i）;
                   else
                       disp（zstep）;
                   end
             % euclidian distance to the reciever plane
             dist_to_rec = z_dist_rec_intersection / uz（i）;    % z / mu_z
             rec_loc（i1） = x（i） + x_dist_rec_intersection;   % x-axis location of reception
             rec_loc（i2） = y（i） + y_dist_rec_intersection;    % y-axis location of reception
             rec_loc（i3） = z（i）;                             % incident angle （mu_x）
             rec_loc（i4） = ux（i）;                             % for statistics should be uniform （mu_y）
             rec_loc（i5） = uy（i）;                             % for statistics should be uniform （mu_z）
             total_rec_packets = total_rec_packets + 1;
             total_rec_power = total_rec_power +weight（i）;        % total power at the receiver plane