We simulated a 1 x 1 x 11 cm3 parallelepiped, with a 1MeV neutron cource directed along 0 0 1 positioned at 0 0 -1. This geometry as cut by 10 planes from 0 to 9 cm, on each plane the neutron spectra was accumulated through the tally:
F2:n 10 11 12 13 14 15 16 17 18 19
e2 1e-10 99log 10
This tally was extracted from mcnp.inpm by using TallyRead.py script. This generated the file mOut/f2 which contains the data that we are going to plot.
In [1]:
# let us import the libraries to read and plot the tally output
import numpy as np # data manipulation library
import matplotlib.pyplot as plt # data plotting library
# things to make the plot prettier...
import matplotlib as mpl
mpl.rcParams['axes.titlesize'] = 22
mpl.rcParams['axes.labelsize'] = 22
mpl.rcParams['lines.linewidth'] = 3
mpl.rcParams['lines.markersize'] = 10
mpl.rcParams['xtick.labelsize'] = 22
mpl.rcParams['ytick.labelsize'] = 22
mpl.rcParams['legend.fontsize'] = 24
In [2]:
# import data into a matrix
matrix = np.loadtxt("mOut/f2")
# crete an array where each column is the result of the tally on one surface
# we have 10 surfaces (0,9). In matrix the first column is the value of the surface
# so we generate an array with 10 empty fields
surfaces = [[] for i in range(10)]
# lets read the matrix into the surface
for i in range(10):
surfaces[i] = matrix[matrix[:,0]==i]
In [3]:
# create the plotting window with pubblication standards
plt.subplots(1,1,figsize=(12, 9), dpi=300)
# load neutron spectra into matlotlib
for surface in surfaces:
# import data with errors
plt.errorbar(surface[:,3],surface[:,4],yerr=surface[:,5], label="%i cm"%(surface[0,0]))
# set scale of the axis
plt.xscale("log")
plt.yscale("log")
# set y range
plt.xlim((1e-11,10))
# set x range
# plt.ylim((5e5,3e8))
plt.xlabel("Neutron Energy (MeV)")
plt.ylabel("Neutron Flux (cm$^{-2}$ s$^{-1}$)")
plt.legend(loc="center left")
plt.show()
In [ ]: