Create Graphs for Bethe-Weizsäcker binding energy per nukleon

bethe_weizsaecker_binding.py

# SPDX-License-Identifier: Unlicense

import matplotlib.pyplot as plt
import numpy as np


# Graph Limits
Zmin = 1
Zmax = 200
Nmin = 1
Nmax = 200

# Bethe-Weizsäcker values in MeV
# from https://de.wikipedia.org/wiki/Bethe-Weizs%C3%A4cker-Formel (2020-05-25)
a_V = 15.67
a_O = 17.23
a_C = 0.714
a_S = 93.15
a_P = 11.2

# Bethe-Weizsäcker formula
# from https://de.wikipedia.org/wiki/Bethe-Weizs%C3%A4cker-Formel (2020-05-25)
def E(Z, N):
    A = Z + N
    if Z % 2 == 0:
        if N % 2 == 0:
            sgn = 1
        else:
            sgn = 0
    else:
        if N % 2 == 0:
            sgn = 0
        else:
            sgn = -1
    out = a_V * A - a_O * A**(2/3) - a_C * Z * (Z-1) * A**(-1/3) - a_S * (N-Z)**2 / (4*A) + sgn * a_P * A**(-1/2)
    if out < 0:
        out = np.nan
    else:
        out /= A
    return out

# Arrays
Z_array = range(Zmin, Zmax+1)
N_array = range(Nmin, Nmax+1)
E_array = [[None for N in N_array] for P in Z_array]

# Computation
for P in Z_array:
    for N in N_array:
        E_array[P-Zmin][N-Nmin] = E(P, N)

# Graph
plt.figure()
plt.title('binding energy per nucleon')
plt.xlim(Nmin, Nmax)
plt.xlabel('N')
plt.ylim(Zmin, Zmax)
plt.ylabel('Z')
graph = plt.pcolormesh(N_array, Z_array, E_array)
plt.colorbar(graph, label='binding energy in MeV')
graph.axes.set_aspect('equal')
plt.grid(True)
plt.tight_layout()
plt.show()