数学建模-核反应堆屏蔽层设计问题

function sim_zhongzi

n=input('中子个数：');%中子个数
xx=[];
yy=[];

N=10;%每个中子最多碰撞次数

d=2;
D=3*d;
H=10*10;

c=zeros(1,3);

x=zeros(1,n);
y=H*rand(1,n)/10;

%中子运动
for j=1:N
    if isempty(x)
        break;
    end
    R=-d*log(rand(1,length(x)));
    seta=2*pi*rand(1,length(x));
    
    
    x=x+(R.*cos(seta));
    y=y+(R.*sin(seta));
    
    %c(2)穿透
    t=find(x>D);
    c(2)=c(2)+length(t);
    
    if length(t)>0
        xx=[xx,x(t)];
        yy=[yy,y(t)];
        
        x(t)=[];
        y(t)=[];
        
    end
    
    %c(1)返回反应堆数量
    t=find(x<0);
    c(1)=c(1)+length(t);
    
    if length(t)>0
        xx=[xx,x(t)];
        yy=[yy,y(t)];
        
        x(t)=[];
        y(t)=[];
        
    end
    
    if j==N
        c(3)=c(3)+length(x);
        
    end
    
end

xx=[xx,x];
yy=[yy,y];

check=sum(c)-n
check=length(xx)-n

bili=c/n

plot(xx,yy,'r')
hold on
line([0,0],[-H,H])
hold on
line([D,D],[-H,H])

text(-4*d,H-5*d,'返回')

text(-4*d,H-9*d,sprintf('(%6.2f%%)',c(1)/n*100))
text(D/2-d,H-5*d,'吸收')
text(D/2-2*d,H-9*d,sprintf('(%6.2f%%)',c(3)/n*100))
text(D+2*d,H-5*d,'穿透')
text(D+d,H-9*d,sprintf('(%6.2f%%)',c(2)/n*100))

hold off
end