Python——链式反应

1、题目

      放射性物质的链式反应是一个随机过程，借助M-C方法模拟研究链式反应，计算其的倍增系数k以及临界质量M。

2、算法及分析

     设轴块为长方体a×a×b，发生裂变的轴核位于轴块内随机点（x0,y0,z0）处。随机点坐标的值域为

      1.给定轴块质量M、轴块边长之比s=ab，旧裂变个数N0，密度为1。则确定：

      2.产生九个0~1的随机数

      3.旧裂变位置

      4.旧裂变放出的两个中子的方向

      5.中子的飞行距离

      6.可能发生新裂变的位置

      7.检查上述位置是否在轴块内，如果在N的值加1

      8.重复2~7，执行N0，计算k

3、程序 

#链式反应
import numpy as np
import random as rd
import math as ma
from matplotlib import pyplot as plt
from scipy import signal
def f(m,s1,N0):
    N=0
    a=(m*s1)**(1/3);b=(m*(s1**(-2)))**(1/3)
    for p in range(N0):
        r=[]
        for i in range(9):
            r.append(rd.random())
        x0=a*(r[0]-0.5);y0=a*(r[1]-0.5);z0=b*(r[2]-0.5)
        dr1=2*ma.pi*r[3];theta1=ma.acos(2*r[4]-1);dr2=2*ma.pi*r[5];theta2=ma.acos(2*r[6]-1)
        d1=r[7];d2=r[8]
        x1=x0+d1*ma.sin(theta1)*ma.cos(dr1)
        y1=y0+d1*ma.sin(theta1)*ma.cos(dr1)
        z1=z0+d1*ma.cos(theta1)
        x2=x0+d2*ma.sin(theta2)*ma.cos(dr2)
        y2=y0+d2*ma.sin(theta2)*ma.cos(dr2)
        z2=z0+d2*ma.cos(theta2)
        if -0.5*a<=x1<=0.5*a and -0.5*b<=y1<=0.5*b and -0.5*a<=z1<=0.5*a:
            N=N+1
        if -0.5*a<=x2<=0.5*a and -0.5*b<=y2<=0.5*b and -0.5*a<=z2<=0.5*a:
            N=N+1
    k=N/N0
    return k
m=[];s1=[];c1=[];p=50;K=50
for jo in range(1,p):
    for ji in range(1,K):
        m.append(jo*0.1)
        s1.append(ji*0.1)
        t1=m[ji-1+(jo-1)*(p-1)];t2=s1[ji-1]
        c1.append(f(t1,t2,10**4))
m=np.array([m])
m=m.reshape(-1,p-1)
s1=np.array([s1])
s1=s1.reshape(-1,p-1)
c1=np.array([c1])
c1=c1.reshape(-1,p-1)
df=0
if type(len(c1)/2)!=type(1):
    df=len(c1)-4
if type(len(c1)/2)==type(1):
    df=len(c1)-3
c1=signal.savgol_filter(c1,df,4)
d3=plt.axes(projection='3d')
v1=np.ones(c1.shape)
d3.plot_surface(m,s1,v1,cmap='rainbow')
d3.plot_surface(m,s1,c1,cmap='rainbow')
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
d3.set_zlabel('倍增系数k')
d3.set_xlabel('质量M/kg')
d3.set_ylabel('轴边长之比s')
d3.set_title('链式反应（经过4阶光滑处理）')
plt.show()	
#平衡态
import random as rd
from matplotlib import pyplot as plt
def f(M):
    t1=[];N=100000;B=0
    for i in range(N):
        t1.append(1)
    for j in range(M):
        r=rd.random()
        k=int(r*N)
        t1[k]=-t1[k]
        B=B-t1[k]
    F=B/N
    return F
x1=[];x2=[]
for op in range(1,100):
    x1.append(f(op*10**4))
    x2.append(op*10**4)
plt.plot(x2,x1,c='k',label='N=100000')
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.xlabel('时间t(s)')
plt.ylabel('粒子数比n/N')
plt.title('平衡态')
plt.legend()
plt.show()

4、计算结果 

FIG1. 链式反应

FIG2. 气体平衡态