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