我用蒙特卡罗法模拟光在沙尘环境中的多次散射问题，编写了如下的程序，程序可以运行，但是感觉运行结果不正确。希望写过类似程序的高手可以指点一下。
fai_T=45*pi/180;%发散角
sita_T=10*pi/180;%发射仰角 45
fai_R=45*pi/180;%接收角
sita_R=10*pi/180;%接收仰角 45
Ksray=0.1;%瑞利散射系数 0.145
Ksmie=1.5;%米氏散射系数 0.261a
Ka=0.6;%吸收函数 0.039
Ks=Ksray+Ksmie;%消光系数
K=Ka+Ks; %吸收散射之和
w=1;%权重
A=1;%存活率0.01
g=0.7;%不对称因子0.2
f=0.5;%散射因子
Ar=1.8*10^(-2);%接收孔径
r=4.8;
gama=0.017;%非对称因子 大气粒子尺寸分布？
L=10;%传输距离
%for i=1:1:1000
%%发射端初始位置%%
w=1;
x=0;y=0;z=0;
ux=0;uy=0;uz=1;
while 1
ksi=unifrnd(0,1);
s=-log(ksi)/K;
x=x+ux*s;y=y+uy*s;z=z+uz*s;
if (z>0)&&(z
w=Ks/K.*w; %改变权重
if w>0.0001%权重是否太小
ksi=unifrnd(0,1);
fai=2*pi*ksi;
temp=(1-g.^2)/(1-g+2*g*ksi);
cos_sita=(1+g.^2-temp.^2)/(2*g);
sin_sita=sqrt(1-cos_sita.^2);
ux1=ux;uy1=uy;uz1=uz;
if abs(uz1)>0.99999
ux=sin_sita*cos(fai);
uy=sin_sita*sin(fai);
uz=uz1/abs(uz1)*cos(fai);
else
temp1=sqrt(1-uz1.^2);
ux=sin_sita*(ux1*uz1*cos(fai)-uy1*sin(fai))/temp1+ux1*cos_sita;
uy=sin_sita*(uy1*uz1*cos(fai)+ux1*sin(fai))/temp1+uy1*cos_sita;
uz=-sin_sita*cos(fai)*temp1+uz1*cos_sita;
end %确定下一个光子的方向
else
ksi=unifrnd(0,1);
if ksi<=1/10
w=10*w;
ksi=unifrnd(0,1);
fai=2*pi*ksi;
temp=(1-g.^2)/(1-g+2*g*ksi);
cos_sita=(1+g.^2-temp.^2)/(2*g);
sin_sita=sqrt(1-cos_sita.^2);
ux1=ux;uy1=uy;uz1=uz;
if abs(uz1)>0.99999
ux=sin_sita*cos(fai);
uy=sin_sita*sin(fai);
uz=uz1/abs(uz1)*cos(fai);
else
temp1=sqrt(1-uz1.^2);
ux=sin_sita*(ux1*uz1*cos(fai)-uy1*sin(fai))/temp1+ux1*cos_sita;
uy=sin_sita*(uy1*uz1*cos(fai)+ux1*sin(fai))/temp1+uy1*cos_sita;
uz=-sin_sita*cos(fai)*temp1+uz1*cos_sita;
end
%确定下一个光子的方向
else
break%此光子结束，追踪下一个光子
end
end
else
sita_i=1/cos(uz);
n_i=1;n_t=2;%ni nt 分别是介质和空气中光子的折射率
sita_t=asin(n_i*sin(sita_i)/n_t); %菲涅尔定律
R=0.5.*((sin(sita_i-sita_t)/sin(sita_i+sita_t)).^2+(tan(sita_i-sita_t)/tan(sita_i+sita_t)).^2);
ksi=unifrnd(0,1);
if ksi
w=w*R;
if z<0
z=-z;
else
z=2*L-z;
end
uz=-uz;
else
w=w*(1-R);
break
end
end
end
%receive(i)=w;
%end
