查看: 2116|回复: 2|关注: 0

[已答复] 用蒙特卡罗法模拟光散射问题

[复制链接]

新手

5 麦片

财富积分


050


1

主题

1

帖子

0

最佳答案
发表于 2016-5-9 16:27:06 | 显示全部楼层 |阅读模式
我用蒙特卡罗法模拟光在沙尘环境中的多次散射问题,编写了如下的程序,程序可以运行,但是感觉运行结果不正确。希望写过类似程序的高手可以指点一下。
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<L)  % ?? 判断光子是否在介质中
    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<R
        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

新手

5 麦片

财富积分


050


0

主题

2

帖子

0

最佳答案
发表于 2018-11-19 14:32:12 | 显示全部楼层
我想问,你写好了吗?能正常运行了吗

新手

5 麦片

财富积分


050


0

主题

1

帖子

0

最佳答案
发表于 2019-3-15 14:11:37 | 显示全部楼层
你运行结果是啥
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

快速回复 返回顶部 返回列表