mat model ver 2.1

import math
import pylab
from matplotlib import mlab


t = 2700.0 + 273.0
T = 2973.0
p = 9678.1 #кг/м^3
Cp = 619.0
small_lambda = 3.64 #теплопроводность
big_lambda = 10196.0
m = 0.61 * 0.000001
beta = 8.555 * 0.00001 #температурный коэффициент линейного расширения
energy = [0.01, 0.1, 0, 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])

#Текущая ячейка - 1
#Две строки, три столбца.
pylab.subplot(2, 2, 1)
pylab.plot(energy, Rc_nat)
pylab.xscale('log')
pylab.yscale('log')
pylab.title("R capture for nature uran")

pylab.subplot(2, 2, 2)
pylab.plot(energy, Rc_o)
pylab.xscale('log')
pylab.yscale('log')
pylab.title("R capture for nature UO2")

pylab.subplot(2, 2, 3)
pylab.plot(energy, Rf_nat)
pylab.xscale('log')
pylab.yscale('log')
pylab.title("R fission for nature uran")

print(*Rf_o)
pylab.subplot(2, 2, 4)
pylab.plot(energy, Rf_o)
pylab.xscale('log')
pylab.yscale('log')
pylab.title("R fission for UO2")

pylab.show()