clear?all;clc;
%本程序不適用于各向異性因子為0的條件
N=10^4;?????%光子數
%起始坐標
x=zeros（1N）;
y=zeros（1N）;
z=zeros（1N）;
%上下邊界
up=10;down=0;
sup=zeros（1N）;
sdown=zeros（1N）;
%各向異性因子
g=0.75;
%介質的吸收系數和散射系數
ua=0.15;???%吸收系數
us=10;??????%散射系數
ut=ua+us;??%相互作用系數
%散射角th與方位角psi
th=pi/2;
psi=0;
%光子能量權重
w=1;
res_w=zeros（1N）;???????%最終光子能量結果保存空間
%傳播方向余弦
ux=zeros（1N）;
uy=zeros（1N）;
uz=ones（1N）;
flg=ux*ux‘+uy*uy‘+uz*uz‘;???%標志，用于判斷monte?carlo模擬是否完成的標準
%檢測器
p=2;r=0.5;
%%%%%%%%%%%%%%%
count=1;
while?flg>1e-10
????count=count+1;
????%光程
????srand=rand（1N）;????????%光程隨機數
????s=-log（srand）/ut;
????%光子能量吸收
????w=w*us/ut;
????%散射角和方位角
????thrand=rand（1N）;???????%散射角隨機數，注：方位角與散射角的隨機量是否為同一個？
????psirand=thrand;?????????%方位角隨機數
????costh=（1+g^2-（（1-g^2）./（1-g+2.*g.*thrand））.^2）/（2*g）;
????sinth=（1-costh.^2）.^0.5;
????psi=2*pi*psirand;
%?????%可加入（頂、底）邊界判斷條件，若判斷溢出，需將對應位置的uxuyuz置零，并將對應位置的能量權重保存至res_w中
%?????????????k=find（uz<0）;kk=find（uz>0）;
%?????????if?~isempty（k）
%?????????????sdown（k）=（z（count-1k）-down）./uz（k）;
%?????????????out=find（s%?????????end
%?????????if?~isempty（kk）
%?????????????sup（kk）=（up-z（count-1kk））./uz（kk）;
%?????????????out=find（s>sup）;%ux（outj）=0;uy（outj）=0;uz（outj）=0;res_w（outj）=w;
%?????????end
????%運動方向與組織表面法向接近
????i=find（uz.*uz>0.99999）;
????if?~isempty（i）
????????ux（i）=sinth（i）.*cos（psi（i））;
????????uy（i）=sinth（i）.*sin（psi（i））;
????????uz（i）=sign（uz（i））.*costh（i）;
????end
????%運動方向與組織表面法向不接近
????j=setdiff（1:Ni）;
????if?~isempty（j）
????????tempx（j）=（sinth（j）.*（ux（j）.*uz（j）.*cos（psi（j））-uy（j）.*sin（psi（j）））./（1-uz（j）.^2）.^0.5）+ux（j）.*costh（j）;
????????tempy（j）=（sinth（j）.*（uy（j）.*uz（j）.*cos（psi（j））+ux（j）.*sin（psi（j）））./（1-uz（j）.^2）.^0.5）+uy（j）.*costh（j）;
????????tempz（j）=（-sinth（j）.*cos（psi（j））.*（1-