# -*- coding: utf-8 -*-
"""
Created on Fri Feb 12 14:59:36 2016

@author: Complex_Batac
"""

import matplotlib.pyplot as plt

def ex2_5(Z):
    # Definition of variables
    a_1 = 15.67
    a_2 = 17.23
    a_3 = 0.75
    a_4 = 93.2

    # Creates empty containers for all of the binding energies per
    # nucleon (B divided by A) and all of the corresponding mass numbers
    bepn = []
    mass_num = []
    # Solves for binding energy per nucleon for each A from A = Z to A = 3Z
    for i in xrange(Z, 3 * Z + 1):
        A_i = i
        if A_i % 2 != 0:
            a_5 = 0.0
        elif A_i % 2 == 0 and Z % 2 == 0:
            a_5 = 12.0
        elif A_i % 2 == 0 and Z % 2 != 0:
            a_5 = -12.0
        B_i = (a_1 * A_i) - (a_2 * A_i**(2.0/3.0)) - \
              (a_3 * Z**2.0 / A_i**(1.0/3.0)) - \
              (a_4 * (A_i - 2.0 * Z)**2.0 / A_i) + \
              (a_5 / (A_i**(1.0/2.0)))
        bepn.append(B_i/A_i)
        mass_num.append(A_i)
    # Prints out the A and B of the most stable nucleus
    print "MOST STABLE NUCLEUS OF ELEMENT WITH Z = " + str(Z)
    print "Mass number: " + str(mass_num[bepn.index(max(bepn))])
    print "Binding energy per nucleon: " + str(max(bepn))

    # Creates empty containers for all atomic numbers, all
    # of the corresponding maximum binding energies per nucleon
    # and all of the corresponding mass numbers
    atomic_num = []
    max_bepn = []
    mass_num_ = []
    # Iterates from Z = 1 to Z = 100
    for m in xrange(1, 101):
        atomic_num.append(m)
        bepn_m = []
        mass_num_m = []
        # Iterates from A = Z to A = 3Z
        for n in xrange(m, 3 * m + 1):
            A_m_n = n
            if A_m_n % 2 != 0:
                a_5 = 0.0
            elif A_m_n % 2 == 0 and m % 2 == 0:
                a_5= 12.0
            elif A_m_n % 2 == 0 and m % 2 != 0:
                a_5 = -12.0
            B_m_n = (a_1 * A_m_n) - (a_2 * A_m_n**(2.0/3.0)) - \
                    (a_3 * m**2.0 / A_m_n**(1.0/3.0)) - \
                    (a_4 * (A_m_n - 2.0 * m)**2.0 / A_m_n) + \
                    (a_5 / (A_m_n**(1.0/2.0)))
            bepn_m.append(B_m_n/A_m_n)
            mass_num_m.append(A_m_n)
        # Exception handling for when bepn_m is an empty sequence
        try:
            max_bepn.append(max(bepn_m))
        except:
            max_bepn.append(0)
        mass_num_.append(mass_num_m[bepn_m.index(max(bepn_m))])
    # Prints the most stable value of A for each Z
    for o in xrange(100):
        print (max_bepn[o], mass_num_[o])
    # Prints the value of Z at which the highest
    # binding energy per nucleon occurs
    print "\n HIGHEST BINDING ENERGY PER NUCLEON AT Z = " + \
          str(atomic_num[max_bepn.index(max(max_bepn))])

    # Graph the binding energy per nucleon for the most stable
    # nucleus of each element as a function of atomic number
    plt.plot(atomic_num, max_bepn)
    plt.xlabel("Atomic number (Z)")
    plt.ylabel("Maximum binding energy per nucleon")
    plt.show()
    plt.close()

ex2_5(40)