mat model ver 3.0

import math
import pylab
from matplotlib import mlab


p = 9678.1 #кг/м^3
m = 0.61 * 0.000001
energy = [0.01, 0.1, 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000]
sigma_t_o = [4, 4, 4, 4, 5, 1, 6, 6, 6, 6]
sigma_c_o = [3*10**(-4), 10**(-4), 3*10**(-5), 10**(-5), 5*10**(-6), 6*10**(-6),
             10**(-5), 3*10**(-5), 10**(-4), 10**(-4)]
sigma_f_238 = [3*10**(-5), 10**(-5), 3*10**(-6), 7*10**(-7), 10**(-6), 3*10**(-4),
               2*10**(-4), 5*10**(-5), 3*10**(-2), 2]
sigma_c_238 = [4, 1,  5*10**(-1), 9*10**(-1), 3*10, 4, 9*10**(-1), 2*10**(-1), 1*10**(-1), 2*10**(-4)]
sigma_f_235 = [10**3, 300, 80, 70, 20, 9, 3, 2, 1, 2]
sigma_c_235 = [200, 40, 10, 70, 10, 4, 1, 0.6, 0.1, 10**(-3)]
for i in range(10):
    sigma_t_o[i] = sigma_t_o[i] * 10**(-28)   #metres^2
    sigma_c_o[i] = sigma_c_o[i] * 10**(-28)
    sigma_c_235[i] = sigma_c_235[i] * 10**(-28)
    sigma_c_238[i] = sigma_c_238[i] * 10**(-28)
    sigma_f_235[i] = sigma_f_235[i] * 10**(-28)
    sigma_f_238[i] = sigma_f_238[i] * 10**(-28)

Na = 6.022*10**23
Mr = ((0.007*235 + 0.993*238)*0.001)
Mr_UO2 = 0.27
N_238_nat = Na*19050/Mr
N_235_nat = Na*18700/Mr
N_o = Na*p/Mr_UO2
Rf_nat = []
Rc_nat = []
Rf_o = []
Rc_o = []
fi1 = 5.0*10**13/0.0001 # 1/(metr^2*sec)
fi2 = 2.5*10**13/0.0001
for i in range(10):
    if energy[i] < 1:
        fi = fi1
    else:
        fi = fi2
    Rf_nat.append(fi * (N_235_nat * sigma_f_235[i] + N_238_nat * sigma_f_238[i]))

for i in range(10):
    if energy[i] < 1:
        fi = fi1
    else:
        fi = fi2
    Rc_nat.append(fi * (N_235_nat * sigma_c_235[i] + N_238_nat * sigma_c_238[i]))

for i in range(10):
    if energy[i] < 1:
        fi = fi1
    else:
        fi = fi2
    Rf_o.append(fi * (32.4 * N_o * sigma_f_238[i] + 0.9 * N_o * sigma_f_235[i]))

for i in range(10):
    if energy[i] < 1:
        fi = fi1
    else:
        fi = fi2
    Rc_o.append(fi * (32.4 * N_o * sigma_c_238[i]+ 0.9 * N_o * sigma_c_235[i])
                + 33.3 * N_o * sigma_c_o[i])


pylab.subplot(2, 2, 1) #Текущая ячейка - 1; Две строки, три столбца
pylab.plot(energy, Rc_nat)
pylab.grid(linestyle='--', linewidth=0.5)
pylab.xlabel('Neutron Energy (eV)')
pylab.ylabel('R capture')
pylab.xscale('log')
pylab.yscale('log')
pylab.title("R capture for nature uran")

pylab.subplot(2, 2, 2)
pylab.plot(energy, Rc_o)
pylab.grid(linestyle='--', linewidth=0.5)
pylab.xlabel('Neutron Energy (eV)')
pylab.ylabel('R capture')
pylab.xscale('log')
pylab.yscale('log')
pylab.title("R capture for UO2")

pylab.subplot(2, 2, 3)
pylab.plot(energy, Rf_nat)
pylab.grid(linestyle='--', linewidth=0.5)
pylab.xlabel('Neutron Energy (eV)')
pylab.ylabel('R fission')
pylab.xscale('log')
pylab.yscale('log')
pylab.title("R fission for nature uran")

pylab.subplot(2, 2, 4)
pylab.plot(energy, Rf_o)
pylab.grid(linestyle='--', linewidth=0.5)
pylab.xlabel('Neutron Energy (eV)')
pylab.ylabel('R fission')
pylab.xscale('log')
pylab.yscale('log')
pylab.title("R fission for UO2")

pylab.show()