sim_norm_circle.m

用MATLAB编写的蒙特卡罗算法,用于检测光子成象.非常不错,大牛编写!!!!开发环境:MATLAB

clear;
%光源坐标
light=[0,0,0];%坐标原点在圆柱的中心
hight=18;%高度的一半  
R=16; %半径

%h=5;%光源长度 
%step=0.49383; %步长

n1=1.37;%液体折射率
n2=1;%外界折射率
rs=3;%探测器半径
hs=[11+rs,11-rs,rs,-rs,-11+rs,-11-rs];%探测器高度范围
%hw=20;%探测器角宽度
%=======================
N=10^3;%光子数
%======================
%临界入射角的余弦值
costh=(1-(n2/n1)^2)^0.5;
%介质的吸收系数和散射系数
ua=0.025;
us=10;
uss=2;
a=uss/(uss+ua);
step=-1/(uss+ua);%步长因子

%光子传输的终止条件中的权重阈值
threshold=10^(-3);
%光子传输的终止条件中自行设置的常数m
m=10;

up=0;%上界面溢出数
down=0;%下界面溢出数
beside=0;%侧面溢出数
disapear=0;
sen(1:24)=0;%探测器接受数

en=[-1,1];
A=zeros(370,370); %用于计算能量分布的矩阵
xi=0;        %用于表示能量矩阵的行
xj=0;        %用于表示能量矩阵的列

tic
clk1=clock;
T=471.4180/100000;

%==================N个光子计算开始======================
for n=1:N
%－－－－－－－－－估计计算剩余时间－－－－－－－－－－
if n==300
    clk2=clock;%确定计算1000次的时间
    T=(clk2-clk1)/300;%计算1次的时间
    T=T(end)+T(end-1)*60;
end
day=double(int64((N-n)*T/3600/24));
hour=double(int64((N-n)*T/3600))-day*24;
minute=double(int64((N-n)*T/60-(hour+day*24)*60));
second=((N-n)*T-(hour+day*24)*3600-minute*60);
fprintf('===========完成百分比%2.5g===========\n\n\n\n\n',n/N*100);%完成百分比
if n<300
    fprintf('请等待完成时间的估算......');
end
if n>300
   fprintf('===========完成百分比%2.5g===========\n===========估计剩余时间 %g 天 %g 小时 %g 分 %2.2g 秒===========\n\n\n\n\n',n/N*100,day,hour,minute,second );%剩余时间
end
if n==301
    fprintf('是否要继续计算，若要请按任意键，若否请按Ctrl+Break\n\n\n\n\n\n\n');
    pause
end
%－－－－－－－－－－－－－－－－－－－－－－－－－－－

%－－－－－－－－确定出射位置，出射方向，初始条件－－－－－－－－
e=rand(1);       %出射位置符合均匀抽样，这里将光源视作点光源，如果视作3D形状，则x，y，z方向都要抽样;
e1=rand(1);e2=rand(1);  %光子出射方向的抽样
fai=2*pi*e1;
costheta=2*e2-1;
sintheta=(1-costheta^2)^0.5;
ux0=sintheta*cos(fai);
uy0=sintheta*sin(fai);
uz0=costheta;

w=1;       %初始权数值
xyz=light;direct=[ux0,uy0,uz0];
info=[xyz(1:3),direct(1:3)];
%－－－－－－－－－计算单个光子过程开始－－－－－－－－－－－－－－－－－－
tou=0;%循环记号，还没透射
jishu=0;
%remain=0;
infotemp=[];
while ( tou==0&jishu<=1000000 )
    jishu=jishu+1;
    for k=1:6
        infotemp(k)=info(k);
    end
    info=iterat(info,step);

%－－－－－－－－能量衰减处理－－－－－－－－－－－－－－－－－－－
w=w*a;
if (w<=threshold)
    %remain=1;
    %有一次机会
    if (rand(1)<1/m)%利用随机数和常数的乘积，判断最后一次传输的光子权重如何
        w=w*m;
    else
        w=0;tou=1;%disapear=disapear+1;%消失
    end
end
			
%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

%计算一个光子散射一次后能量对整个能量矩阵的影响，划分成0.1*0.1的小网格
%if(info(3)>en(1)&&info(3)<en(2))  %en
%       xi=int16((10*info(2)));   %取整
%		xi=xi+160;          %去掉负值的影响
%		xj=int16((10*info(1)));
%		xj=xj+160;
%       if(xi>0&xi<=320&xj>0&xj<=320)	
%          A(xi,xj)=A(xi,xj)+(1/a-1)*w*(1-remain*(m-1)/m);
%      end
% end

%remain=0;

%－－－－－－－－－边界处理－－－－－－－－－－－－－－－－－－－－
if (info(3)<(-hight)) %down
    tempp=abs(info(6));%夹角的余弦值
    if (tempp>costh)
        tou=1;
        down=down+1;%透射
    else
        info(6)=abs(info(6));%反射
        info(3)=hight*hight/info(3);
    end
else
    if (info(3)>hight) %up
        tempp=abs(info(6));%夹角的余弦值
        if (tempp>costh)
            tou=1;
            up=up+1;%透射
        else
            info(6)=-abs(info(6));
            info(3)=hight*hight/info(3);%反射
        end
    else
        if (sqrt(info(1)*info(1)+info(2)*info(2))>R) %sqrt(Norm(info,2))>R-step/10)
            tempp=sqrt((info(1))^2+(info(2))^2+(info(3))^2);
			temp1=info(1)/tempp;
			temp2=info(2)/tempp;
			temp3=info(3)/tempp;
			tempp=temp1*info(4)+temp2*info(5)+temp3*info(6);%夹角的余弦值
            if (tempp>costh)
                tou=1;
                beside=beside+1;%透射
            %jiao=atan2(info(2),info(1))/pi*180;
			    temp4=0;j=0;
                while (temp4==0&j<=7)
                    j=j+1;
                    %if(jiao>((j-1)*45-180)&jiao<((j-1)*45+hw-180))
                    jiao=((j-1)*45-180)/180*pi;
                    for k=1:3
                        if (info(3)<hs(2*k-1)&info(3)>hs(2*k))  %hs 
                            temp6=sqrt((info(1))^2+(info(2))^2+(info(3)-hs(2*k-1)/2-hs(2*k)/2)^2);
							temp6=(cos(jiao)*info(1)+sin(jiao)*info(2))/temp6;
                            if (temp6>(R/sqrt(R*R+rs*rs)))
                                temp4=1; %探测器接收到了
                                sen(j-8+8*k)=sen(j-8+8*k)+1;
                                break;
                            end
                        end
                    end
                end
           else
            temp4=info(4);temp5=info(5);%反射
           
		    info(4)=(temp2*temp2-temp1*temp1)*temp4-2*temp1*temp2*temp5;
		    info(5)=(temp1*temp1-temp2*temp2)*temp5-2*temp1*temp2*temp4;

			tempp=info(1)*info(1)+info(2)*info(2);
			info(1)=R*R/tempp*info(1);
			info(2)=R*R/tempp*info(2);
            end
        end
    end
end

end
end


toc
N
[up,down,beside]
sen
%surf(1:370,1:370,A);
axis equal