我用蒙特卡罗法模拟光在沙尘环境中的多次散射问题，编写了如下的程序，程序可以运行，但是感觉运行结果不正确。希望写过类似程序的高手可以指点一下。

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