In [7]:

    

import numpy as np
import matplotlib.pyplot  as plt


    

In [8]:

    

# Assign indices for each of the four reactions
i_2Hdn  = 0 #  2H(d,n)3He "D-D"
i_2Hdp  = 1 #  2H(d,p)3H  "D-D"
i_3Hdn  = 2 #  3H(d,n)4He "D-T"
i_3Hedp = 3 # 3He(d,p)4He "D-He3"
reactions = (i_2Hdn, i_2Hdp, i_3Hdn, i_3Hedp)

# Assign numerial values for all the B-H parameters. Notation 
# ("A","B","B_G") corresponds exactly as found in the original
# B-H reference for clarity.

# Allocate the Bosch-Hale coefficients for each fusion reaction
B_G_2Hdn  = 31.3970
B_G_2Hdp  = 31.3970
B_G_3Hdn  = 34.3827
B_G_3Hedp = 68.7508
B_G = (B_G_2Hdn, B_G_2Hdp, B_G_3Hdn, B_G_3Hedp)

#                     A_1       A_2        A_3        A_4         A_5
A_2Hdn  = np.array([5.3701e4, 3.3027e2, -1.2706e-1, 2.9327e-5, -2.5151e-9 ])
A_2Hdp  = np.array([5.5576e4, 2.1054e2, -3.2638e-2, 1.4987e-6,  1.8181e-10])
A_3Hdn  = np.array([6.9270e4, 7.4540e8,  2.0500e6,  5.2000e4,   0.0000e0  ])
A_3Hedp = np.array([5.7501e6, 2.5226e3,  4.5566e1,  0.0000e0,   0.0000e0  ])
A = (A_2Hdn, A_2Hdp, A_3Hdn, A_3Hedp)

#                      B_1         B_2        B_3        B_4
B_2Hdn  = np.array([ 0.000    ,  0.000    , 0.000    , 0.000   ])
B_2Hdp  = np.array([ 0.000    ,  0.000    , 0.000    , 0.000   ])
B_3Hdn  = np.array([ 6.380e1  , -9.950e-1 , 6.981e-5 , 1.728e-4])
B_3Hedp = np.array([-3.1995e-3, -8.5530e-6, 5.9014e-8, 0.000   ])
B = (B_2Hdn, B_2Hdp, B_3Hdn, B_3Hedp)

# Assign the valid energy ranges for each reaction [E_min, E_max] [keV]
En_2Hdn = np.array([0.5, 4900]) # 
En_2Hdp = np.array([0.5, 5000])
En_3Hdn = np.array([0.5, 550])
En_3Hedp = np.array([0.3, 900])
En = (En_2Hdn, En_2Hdp, En_3Hdn, En_3Hedp)


    

In [9]:

    

# Compute the B-H cross sections. The cross section will be arrays
# corresponding to cross section as functions of energy

# List to hold the cross-section arrays for each reaction
sigma = []

# Iterate over each of the four reactions
for r in reactions:
    E = np.logspace(0, 3, 500) # [keV]
    S_n = A[r][0] + E*(A[r][1] + E * (A[r][2] + E * (A[r][3] + E * A[r][4])))
    S_d = 1 + E * (B[r][0] + E * (B[r][1] + E * (B[r][2] + E * B[r][3])))
    S = S_n / S_d
    cs = S / (E * exp(B_G[r] / sqrt(E))) # [millibarn]

    # Transform units from [millibarn] to [m^2]
    millibarn_2_barn = 1e-3
    barn_2_m2 = 1e-28
    cs = cs * millibarn_2_barn * barn_2_m2 # [m2]

    # Append this reaction cross section to container
    sigma.append(cs) # [m2]


    

In [10]:

    

# Draw all four cross sections on the same log-log plot over
# reasonably useful energy ranges

fig, ax = plt.subplots()

plt.hold(True)
ax.loglog(E, sigma[0], '-bo', linewidth=2, markersize=14, markevery=30, label='$\mathrm{^2\!H(d,n)^3\!He}$')
ax.loglog(E, sigma[1], '-k^', linewidth=2, markersize=12, markevery=30, label='$\mathrm{^2\!H(d,p)^3\!H}$')
ax.loglog(E, sigma[2], '-rs', linewidth=2, markersize=11, markevery=30, label='$\mathrm{^3\!H(d,n)^4\!He}$')
ax.loglog(E, sigma[3], '-g*', linewidth=2, markersize=15, markevery=30, label= '$\mathrm{^3\!He(d,p)^4\!He}$')
ax.hold(False)

ax.legend(loc='lower right', fontsize=24, numpoints=1)

font = {'family' : 'sans-serif',
        'style'  : 'normal',
        'weight' : 'normal',
        'size'   : 28}
matplotlib.rc('font', **font)

ax.grid(axis='both')
plt.title('Bosch-Hale Nuclear Reaction Cross Sections', fontsize=24)
plt.xlabel('Energy [keV]')
plt.ylabel('Cross section [ m$^2$]')

ax.xaxis.set_tick_params(which='major', width=1, length=10)
ax.xaxis.set_tick_params(which='minor', width=1, length=5)
ax.yaxis.set_tick_params(which='major', width=1, length=10)
ax.yaxis.set_tick_params(which='minor', width=1, length=5)

pylab.xlim([1,1000])
pylab.ylim([1e-35,1e-27])

gcf().subplots_adjust(top=0.92, bottom=0.2, left=0.2, right=0.95)

fig.set_size_inches(10,8)

plt.show()


    

In [11]:

    

# Save the figure as a vector-graphic EPS file
fig.savefig('BH_CrossSections.eps')