核反应堆屏蔽层设计说明

./2011高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号〔由赛区组委会评阅前进行编号：赛区评阅记录〔可供赛区评阅时使用：评阅人评分备注全国统一编号〔由赛区组委会送交全国前编号：全国评阅编号〔由全国组委会评阅前进行编号：.核反应堆屏蔽层设计摘要核反应堆,又称为原子反应堆或反应堆,是装配了核燃料以实现大规模可控制裂变链式反应的装置。为防护中子、γ射线和热辐射,必须在反应堆和大多数辅助设备周围设置屏蔽层。通常设计最外层屏蔽时应将辐射减到人类允许剂量水平以下,常称为生物屏蔽。核反应堆屏蔽层是用一定厚度的铅把反应堆四周包围起来,用以阻档或减弱反应堆发出的各种射线。在各种射线中,中子对人体伤害极大,因此,屏蔽设计,主要是了解中子穿透屏蔽的百分比<或概率>,这对反应堆的安全运行是至关重要的。问题一在假定屏蔽层厚度已知的情况下,要求求解中子穿透屏蔽层的百分比。采用传统的的蒙特卡罗法对问题进行分析,利用软件计算出当屏蔽层时的时的概率。首先假定屏蔽层是理想的均匀铅质平板；中子在屏蔽层相继两次碰撞之间游动的距离服从指数分布；两次碰撞之间中子的平均游动距离为；中子在轴上的位置表示为,中子运动方向与轴夹角为,即用描述第次碰撞后的中子的运动状态；中子碰撞后的弹射角服从的均匀分布；从而得到第次碰撞后在屏蔽层的位置。由假设可能得到式中为两次碰撞之间中子的平均游动距离；是区间上均匀分布的随机数。然后采用逆转法求解平均分布和指数分布的随机数。问题二用蒙特卡罗模拟要求穿透率为时,需试验次。为节省时间,减少计算,在问题一上做改进。采用蒙特卡罗模拟收敛法,先用分析法求出穿透率为时屏蔽层厚度,则三倍此厚度即可达到穿透率不大于的要求。再用计算机收索法进行验证,得出相同结果且符合实际要求。最后根据上述估计为日本核危机善后提出建议。关键词：蒙特卡罗模拟分析法计算机收索法一、问题重述1.1问题背景核反应堆,又称为原子反应堆或反应堆,是装配了核燃料以实现大规模可控制裂变链式反应的装置。为防护中子、γ射线和热辐射,必须在反应堆和大多数辅助设备周围设置屏蔽层。有些屏蔽材料俘获中子后放射出γ射线,因此在中子屏蔽外要有一层γ射线屏蔽。通常设计最外层屏蔽时应将辐射减到人类允许剂量水平以下,常称为生物屏蔽。核反应堆屏蔽层是用一定厚度的铅把反应堆四周包围起来,用以阻档或减弱反应堆发出的各种射线。在各种射线中,中子对人体伤害极大,因此,屏蔽设计,主要是了解中子穿透屏蔽的百分比<或概率>,这对反应堆的安全运行是至关重要的。1.2问题提出⑴假设屏蔽层,在大数定理的意义下,中子穿透屏蔽层的百分比多少？⑵在实际应用中,要求中子穿诱屏蔽层的概率极小,一般数量级为,即穿入屏蔽层的中子若为几百万个,也只能有几个中子穿过屏蔽层。问题是多厚的屏蔽层才能使它被穿的概率小于?⑶根据上述估计,并查阅相关文献,尝试为日本福岛核泄漏事件的核危机善后工作提出约2000字的建议。二、模型假设⑴假定屏蔽层是理想的均匀铅质平板；⑵假定中子在屏蔽层相继两次碰撞之间游动的距离服从指数分布；⑶两次碰撞之间中子的平均游动距离为；⑷中子在轴上的位置表示为,中子运动方向与轴夹角为；⑸中子碰撞后的弹射角服从的均匀分布；符号说明：精确度：单层屏蔽厚度：试验次数：两次碰撞间中子游动距离：穿透率：单层穿透率：随机数：为中子在轴上的位置：中子运动的方向和轴的夹角模型的建立与求解4.1问题一问题分析首先考虑一个中子进入屏蔽层后运动的物理过程：假定屏蔽层是理想的均匀铅质平板,中子以初速度和方向角射入屏蔽层,如图1所示,运动一段距离后,在处与铅核碰撞之后,中子获得新的速度及方向,再运动一段距离后,与铅核第二次碰撞,并获得新的状态等等,经若干次碰撞后,发生以下情况之一则终止运动过程：<1>弹回反应堆；<2>穿诱屏蔽层；<3>第i次碰撞后,中子被屏蔽层吸收。图1蒙特卡罗法作为一种计算方法,是由S.M.乌拉姆和J.·诺伊曼在20世纪40年代中叶为研制核武器的需要而首先提出来的。在此之前,该方法的基本思想实际上早已被统计学家所采用了。例如,早在17世纪,人们就知道了依频数来决定概率的方法。采用传统的的蒙特卡罗法对问题进行分析,利用软件计算出当屏蔽层时的概率。使用蒙特卡罗方法进行分子模拟计算是按照以下步骤进行的：⑴使用随机数发生器产生一个随机的分子构型。⑵对此分子构型的其中粒子坐标做无规则的改变,产生一个新的分子构型。⑶计算新的分子构型的能量。⑷比较新的分子构型于改变前的分子构型的能量变化,判断是否接受该构型。⑸若新的分子构型能量低于原分子构型的能量,则接受新的构型,使用这个构型重复再做下一次迭代。若新的分子构型能量高于原分子构型的能量,则计算玻尔兹曼常数,同时产生一个随机数。⑹若这个随机数大于所计算出的玻尔兹曼因子,则放弃这个构型,重新计算。⑺若这个随机数小于所计算出的玻尔兹曼因子,则接受这个构型,使用这个构型重复再做下一次迭代。⑻如此进行迭代计算,直至最后搜索出低于所给能量条件的分子构型结束。模型的建立假定屏蔽层平行板得厚度为,其中为两次碰撞之间中子的平均游动距离；每次碰撞后中子因损失一部分能量而速度下降,假设第十次碰撞后,中子的速度下降到某一很小的数值而终止运动〔被吸收。由于对穿透屏蔽层的中子感兴趣,故用描述第次碰撞后的中子的运动状态,其中为中子在轴上的位置,中子运动的方向和轴的夹角。中子在再屏蔽层相继两次碰撞之间游离的距离服从指数分布,中子碰撞后的弹射角服从上的均匀分布。从而得到第次碰撞后在屏蔽层的位置⑴式中,是中子第次碰撞后的弹射角度,是中子从第次碰撞到第次碰撞时所游动的距离。由假设可能得到：式中为两次碰撞之间中子的平均游动距离；是区间上均匀分布的随机数。式⑴表明了中子运动的位置和方向是随机的。生成连续随机变数的方法很多。选择哪些算法取决于我们所希望生成变数的的分布,需要考虑到随机变量的准确性、计算和存储效率以及算法的复杂性等因素。最常用的两种算法是逆转法和接受-排除法。下面采用逆转法求解平均分布和指数分布的随机数：设指数分布的密度函数为：分布函数为：生成随机数,设求解得。由此可得重新整理后得：取方程两边的自然对数,可得最后求解可得为了简化计算用代替。由于是随机数,所以也是随机数,因此指数分布的过程生成器将是：分析在区间上均匀分布的随机变量,该分部的由下列函数给出：用法从这个随机变量生成观测值,分布函数求得：要使用从均匀分布生成观测值,首先生成随机数,然后设以求解。由此可得：求解可得：它是均匀分布的生成器。可以用下述公式对进行模拟：式中：是均匀分布随机数。中子在屏蔽层里随机游动,类似布朗运动,第次碰撞后,中子的位置有三种情况：⑴,中子返回反应堆；⑵,中子穿出屏蔽层；⑶若,则中子在屏蔽层继续运动；若,则中子被屏蔽层吸收。中子的三种运动模式如图2。针对中子的三种运动模式,利用程序模拟求解,模拟结果如下表1：中子数〔个返回〔%吸收〔%穿透〔%10082.00%2.00%12.00%100081.60%5.80%12.60%300081.036.73%12.23%500081.28%6.30%12.42%1000080.98%6.47%11.40%1500081.69%6.30%11.97%2000081.50%5.83%12.67%10000081.75%6.18%12.07%中子扩散图像如下：中子穿透屏蔽层的百分比约为12%。4.2问题二模型的建立从表1可知,当要求精确度时,试验次数。在实际应用中,要求穿透概率很小,比如,即是说,穿入屏蔽层的中子数若为几百万个,只允许几个穿透屏蔽层。此时,要达到千分之一的精度,则试验次数为,要获得更高的精度,则问题一的模拟法已行不通。现在,讨论加快模拟收敛的方法。把均匀平行板分为厚度相同的层,只取一层进行模拟。设中子在一层中被吸收和弹回的概率之和为,则穿过一层的概率是,因此穿透层的概率是。由于中子穿过一层的平均速度有所下降,因而总的穿透概率。对方法,试验次数可保证的精度好于,经层后有,若,就可获得穿透概率。这样处理后,不必作高达的试验,只需作次的试验就可达到的精度。这一改进比直接的模拟方法大大加快了收敛速度,从而大大减少了模拟时间。模型求解分析法求解：由于模型模拟过程中设计较多中子,又每个中子碰撞最多达次,因此模拟万个中子,则实际运行最多要处理万次碰撞。通过循环语句逐个处理则费时,因此采用加快模拟收敛的方法,可以利用的功能,编写成值文件,提高模拟效率。模拟得当时模拟的结果如下表所示：中子数/个返回/%穿透/%吸收/%100080.413.16.550008112.786.221000081.612.116.29表格2：模拟结果图表1：模拟结果图示模拟结果表明取屏蔽层厚不合适,中子穿透屏蔽层的百分比超过了,在实际应用中要求穿透屏蔽的概率很小,数量级为,即进入屏蔽层的中子若为百万个以上,只允许几个穿过屏蔽层。设计屏蔽层的厚度,可把屏蔽层看成层厚度均为的平行板,由于中子穿过一层后的平均速度有所下降,若记中子穿过厚度屏蔽层的概率为,则穿过整个屏蔽的概率有：利用前面的模拟结果：时,,令,即：即：因此取蔽层的厚度时,可使中子穿透屏蔽层的概率。经过三层屏蔽后,则中子穿透概率,整个屏蔽层的厚度为。计算机搜索法验证：逐步增大屏蔽层厚度,如：屏蔽层厚度开始取,若则增加厚度为,直到穿透概率为止。经过三层屏蔽后,则可使中子穿透屏蔽层的概率。运用模拟结果如下：厚度返回/%穿透/%吸收/%82.2211.656.1381.982.4915.5381.760.817.4482.81017.19表格3：模拟结果图表2：模拟结果图示验证结果与分析法相符合,因此,经过三层屏蔽后,整个屏蔽层的厚度为,则中子穿透概率。建议报告2011年3月11日日本本州岛附近海域发生强烈地震,强震导致福岛第一核电站发生爆炸引发核危机,从这一刻开始就引发了核泄漏。福岛第一核电站的一号、三号和二号反应堆分别在当地时间日时分、日时分和日晨时分发生氢气爆炸。福岛发生的核泄漏事件死亡人、失踪人,辐射泄漏最初一星期向大气释放的辐射物质活度高达至万万亿贝克勒尔,对海洋领域造成的损失不可估计。福岛核电站从选址到设计都有严格质量监控,对抗震也有很高的技术要求,福岛核电站发生核泄漏的可能原因是日本遭受的特震超过了当初的抗震标准。另一方面核反应堆在停堆后仍需要对堆芯进行冷却,因为核燃料有自衰变余热,若长时间得不到冷却,会使得堆芯温度上升,导致核燃料棒熔化。一般停堆后,冷却系统泵所需的电力需要从外部输入,同时备有多台应急发电机供电。但日本强震后,福岛核电厂外电网瘫痪,应急柴油发电机短暂运行后也陷入中断状态。为防止核反应堆大量放射线从堆芯穿透出来,一般反应堆设有三到四道屏障。第一道是把核燃料放在瓷芯块中,使得大部分裂变产物和气体产物保存在芯块中；第二道是合金制造的包壳,把瓷芯块密封在其中,使得其在高温下不与水发生反应；第三道是耐压力容器,避免放射性物质泄漏到反应堆厂房；第四道是反应堆安全壳,一般是厚近厘米的钢筋混凝土墙,里面衬有防辐射金属材料,防止放射性物质进入环境。安全壳是最重要的一道屏障,以防前面的措施失效,它能保证周围居民遭受辐射在安全围。核反应堆,又称为原子反应堆或反应堆,是装配了核燃料以实现大规模可控制裂变链式反应的装置。制裂变链式反应的装置。反应堆主要有活性区,反射层,外压力壳和屏蔽层组成。活性区又由核燃料,慢化剂,冷却剂和控制棒等组成；反射层是为减少堆芯的中子泄漏、降低临界质量和临界尺寸、提高并展平堆芯中子注量率而布置在堆芯四周的材料或物体。应具有大的中子散射截面和小的中子吸收截面；屏蔽层是为防护中子、射线和热辐射,必须在反应堆和大多数辅助设备周围设置屏蔽保护层。其中慢化剂、控制棒、冷却剂、屏蔽层是核心组件。慢化剂就是用来将快中子能量减少,使之慢化成为中子或中能中子的物质,常采用石墨、重水、轻水；控制棒在反应堆中起补偿和调节中子反应性以及紧急停堆的作用；冷却剂由主循环泵驱动,在一回路中循环,从堆心带走热量并传给二回路中的工质,使蒸汽发生器产生高温高压蒸汽,以驱动汽轮发电机发电,是既在堆心中工作又在堆外工作的一种反应堆成分,这就要求冷却剂必需在高温和高中子通量场中工作是稳定的。屏蔽层是为防护中子、射线和热辐射,必须在反应堆和大多数辅助设备周围设置屏蔽保护层,其设计要力求造价便宜并节省空间。目前日本福岛核电站的放射性物质控制在安全壳,仍未造成真正意义上的核扩散,但一旦安全壳密闭失效将会造成更为严重的后果。根据问题一、问题二的求解可知核反应堆的屏蔽层设计的厚度严重影响了中子穿透屏蔽层的百分比。当屏蔽层的厚度〔为两次碰撞间中子的平均游动距离时,中子穿透屏蔽层的百分比为,远远大于实际应用中要求穿透屏蔽层的概率。在实际应用中要求穿透屏蔽层的概率的数量级为,即进入屏蔽层的中子若为几百万个以上,只允许几个穿透屏蔽层。于是用蒙特卡罗方法把屏蔽层看成层厚度均为的平行板,得到当屏蔽层的厚度时,可是穿透屏蔽层的概率达到实际应用要求。可见增加屏蔽层的厚度可以在很大程度上减小放射性物质的泄露。一旦发生核泄漏,应及时用人工方式使反应堆冷却,通过注入海水是堆芯降温；同时应采取措施防止冷却水的流失,封闭核污染区,保证人生安全和环境清洁,及时疏散周边人员。日本政府应加大对1、2、3、4号机组的监测与处理工作,以反应堆安全壳的破坏造成更大规模的核泄漏和更严重的损失。年发生的史上最大的切尔诺贝利核电站核泄漏事故导致事故后前个月有人死亡,之后年有万人死亡,万人遭受各种程度的辐射疾病折磨,方圆公里地区的万多民众被迫疏散；年日本福岛发生的核泄漏事件死亡人、失踪人,辐射泄漏最初一星期向大气释放的辐射物质活度高达至万万亿贝克勒尔,对海洋领域造成的损失不可估计。核泄漏造成的严重损失以及对人类生命财产的致命威胁,为我们敲响了要和平利用核能源与预防核泄漏的警钟。如何有效地预防核泄漏以及在核泄漏发生后应采取怎样的措施使核泄漏造成的威胁和伤害达到最低,这就需要我们对核反应堆、核泄漏发生的原因和核辐射有深入的了解。通过福岛核泄露的原因分析,需要加强核反应堆的防震功能来应对地震对核反应堆的破坏；同时要设有多台应急发电机,以防地震造成的电网瘫痪而是反应堆堆芯无法冷却；加强反应堆停堆系统、应急冷却系统和安全壳等安全措施的发展是其在发生核泄漏时能发挥应有的作用,从而保证核能的安全开发和利用。六、模型评价6.1模型评价优点问题一,蒙特卡罗法作为一种计算方法,是由S.M.乌拉姆和J.·诺伊曼在20世纪40年代中叶为研制核武器的需要而首先提出来的。在此之前,该方法的基本思想实际上早已被统计学家所采用了。例如,早在17世纪,人们就知道了依频数来决定概率的方法。属于比较成熟的分析方法,所得结果比较可靠。问题二,有些系统构造虽不复杂,却属部机理不明的"黑箱"系统,我们无法用机理分析建立模型,系统模拟是很最重要的一种方法。每模拟一次中子的运动就相当于对中子运动过程做一次"试验"或"观察",按理论来讲必须模拟无限次,才能求出中子穿出屏蔽层的概率。在实际中常进行较大次数的模拟,用所得频率作为所求概率的近似解。设计模拟模型是为了"观察"系统并收集必要信息,这和观察一个实现系统非常相似,模拟系统实现系统的仿真。6.1.2缺点模拟模型本质是一种求解问题的试验方法,是对思维结果的一种验证。由于有随机数每次运行结果不同,与实际有偏差。方法比较单一,运行的精确度虽然已经很高,但是对于实际情况来说仍然需要提高,特别是在当今倡导人权的社会中,需要给予核设施附近的工作人员更好的保护,因此需要将精确度进一步提高。七、参考文献[1]阮晓青,周义仓,数学建模引论,：高等教育,20XX。[2]袁新生,邵大宏,郁时炼,：科学,20XX。[3]吕金虎,陆君安,士华,：大学,20XX。[4]乃文,：清华大学,20XX。八、附录附录一：值文件functionsim_zhongzi2%模拟核反应堆屏蔽层中中子的运动%为了提高程序运行速度,这里没有逐个动态显示中子的运动,%所以一开始就产生所有中子的初始状态；当然,也可以很容易修改为动态显示程序。%n=input<'中子个数：'>;%中子个数xx=[];%存储第10次碰撞后中子的横坐标yy=[];%存储第10次碰撞后中子的纵坐标N=10;%每个中子最多碰撞次数d=2;D=3*d;H=10*10;c=zeros<1,3>;%c<1>返回反应堆数量,c<2>穿透,c<3>吸收%产生中子x=zeros<1,n>;y=H*rand<1,n>/10;%中子运动forj=1:N,%考虑至多10次碰撞ifisempty<x>,break;endR=-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>;iflength<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>;iflength<t>>0,xx=[xx,x<t>];yy=[yy,y<t>];x<t>=[];y<t>=[];endifj==N,c<3>=c<3>+length<x>;endendxx=[xx,x];yy=[yy,y];check=sum<c>-n%right-0check2=length<xx>-n%right-0bili=c/n%return%将此行注释或删除后才能执行下列代码%重绘所有中子,当n增大时,会花较长时间plot<xx,yy,'r.'>holdonline<[0,0],[-H,H]>holdonline<[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>>holdoff附录二：程序subplot<2,2,1>n=10000;xx=[];%存储第10次碰撞后中子的横坐标yy=[];%存储第10次碰撞后中子的纵坐标N=10;%每个中子最多碰撞次数d=2;D=3*d;H=10*10;c=zeros<1,3>;%c<1>返回反应堆数量,c<2>穿透,c<3>吸收%产生中子x=zeros<1,n>;y=H*rand<1,n>/10;%中子运动forj=1:N,%考虑至多10次碰撞ifisempty<x>,break;endR=-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>;iflength<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>;iflength<t>>0,xx=[xx,x<t>];yy=[yy,y<t>];x<t>=[];y<t>=[];endifj==N,c<3>=c<3>+length<x>;endendxx=[xx,x];yy=[yy,y];check=sum<c>-n%right-0check2=length<xx>-n%right-0bili=c/n%return%将此行注释或删除后才能执行下列代码%重绘所有中子,当n增大时,会花较长时间plot<xx,yy,'r.'>holdonline<[0,0],[-H,H]>holdonline<[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>>holdoffsubplot<2,2,2>n=10000;xx=[];%存储第10次碰撞后中子的横坐标yy=[];%存储第10次碰撞后中子的纵坐标N=10;%每个中子最多碰撞次数d=2;D=6*d;H=10*10;c=zeros<1,3>;%c<1>返回反应堆数量,c<2>穿透,c<3>吸收%产生中子x=zeros<1,n>;y=H*rand<1,n>/10;%中子运动forj=1:N,%考虑至多10次碰撞ifisempty<x>,break;endR=-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>;iflength<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>;iflength<t>>0,xx=[xx,x<t>];yy=[yy,y<t>];x<t>=[];y<t>=[];endifj==N,c<3>=c<3>+length<x>;endendxx=[xx,x];yy=[yy,y];check=sum<c>-n%right-0check2=length<xx>-n%right-0bili=c/n%return%将此行注释或删除后才能执行下列代码%重绘所有中子,当n增大时,会花较长时间plot<xx,yy,'r.'>holdonline<[0,0],[-H,H]>holdonline<[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>>holdoffsubplot<2,2,3>n=10000;xx=[];%存储第10次碰撞后中子的横坐标yy=[];%存储第10次碰撞后中子的纵坐标N=10;%每个中子最多碰撞次数d=2;D=8*d;H=10*10;c=zeros<1,3>;%c<1>返回反应堆数量,c<2>穿透,c<3>吸收%产生中子x=zeros<1,n>;y=H*rand<1,n>/10;%中子运动forj=1:N,%考虑至多10次碰撞ifisempty<x>,break;endR=-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>;iflength<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>;iflength<t>>0,xx=[xx,x<t>];yy=[yy,y<t>];x<t>=[];y<t>=[];endifj==N,c<3>=c<3>+length<x>;endendxx=[xx,x];yy=[yy,y];check=sum<c>-n%right-0check2=length<xx>-n%right-0bili=c/n%return%将此行注释或删除后才能执行下列代码%重绘所有中子,当n增大时,会花较长时间plot<xx,yy,'r.'>holdonline<[0,0],[-H,H]>holdonline<[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>>