monte calo模拟

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-