"""
   This is a simple python script that reads the file binding_energies.txt and plots binding energies of nuclei in various fomats.

   It also computes the binding energy according to von Weizsaecker's droplet model 

   author: Bob Wimmer, CAU Kiel, Germany
   email:  wimmer@physik.uni-kiel.de
   date:   April 17, 2015
   
 Copyright (C) 2015  Bob Wimmer

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You can find a copy of the GNU General Public License
    at <http://www.gnu.org/licenses/>.
"""
import sys
from numpy import *
import matplotlib.pyplot as plt

def droplet(A,Z):
    #constants for droplet model in MeV
    av = 15.84 #volume term
    ao = 18.33 #survace term
    ac = 0.714 #coulomb term
    ap = 11.2  #pair term
    aa = 23.2  #asymmetry term
    def dp(A,Z):
#        if (A%2 == 0) and (Z%2 == 0): 
#            return 1.
#        elif (A*Z%2 == 1):
#            return 0.
#        else: return -1.
        return 1 - (A-Z)%2 - Z%2
    ebd = av*A - ao*A**(2./3) - aa*((A-2*Z)**2)/A - ac*Z**2/A**(1./3) + ap*dp(A,Z)/sqrt(A)
#    print A, Z, dp(A,Z)
    return ebd

def plot_ebtpn_vs_A(ebtpn,A,Z):
    plt.plot(A,ebtpn,'b+')
    plt.xlabel('atomic mass number (A) []')
    plt.ylabel('binding energy per nucleon [MeV/nuc]')
    plt.show()

def plot_ebt_vs_ebd(ebt,ebd):
    plt.plot(ebd,ebt,'b+')
    plt.xlabel('droplet model binding energy [MeV]')
    plt.ylabel('measured binding energy [MeV]')
    plt.show()

def plot_ebtpn_vs_ebdpn(ebt,ebd,A):
    plt.plot(ebd/A,ebt/A,'b+')
    plt.xlabel('droplet model binding energy/nuc [MeV/nuc]')
    plt.ylabel('measured binding energy/nuc [MeV/nuc]')
    plt.show()

def plot_deb_vs_N(ebt,ebd,A,Z):
    plt.plot(A-Z,ebd-ebt,'b+')
    plt.xlabel('neutron number []')
    plt.ylabel('E_{droplet} - E_{meas} [MeV]')
    plt.show()


if  __name__ == '__main__':
    print('usage: binding_energies.py A Z')
    print('outputs binding energy for nucleus with mass number A and nuclear charge Z')
    print('or:    binding_energies.py')
    print('reads file and plot comparisons')
    print(len(sys.argv))
    if len(sys.argv) == 3: 
        A = float(sys.argv[1])
        Z = float(sys.argv[2])
        print('Droplet model binding energy for nucleus ({0:d},{1:d}) is: {2:6.3f} MeV'.format(int(A),int(Z),droplet(A,Z)))
    else:
        A, Z, El, M, ebt, ebtpn = loadtxt('binding_energies.txt',dtype={'names':('A', 'Z', 'El', 'M', 'ebt', 'ebtpn'),'formats':(float,float,'|S5',float,float,float)},skiprows=9,unpack=True)
        print(A[10], Z[10], El[10])
        plot_ebtpn_vs_A(ebtpn,A,Z)
        ebd = []
        for i in arange(0,len(A)):
            ebd.append(droplet(A[i],Z[i]))
        plot_ebt_vs_ebd(ebt,ebd)
        plot_ebtpn_vs_ebdpn(ebt,ebd,A)
        plot_deb_vs_N(ebt,ebd,A,Z)
        print('correlation between measured binding energy and droplet model is:' )
        print(corrcoef(ebt,ebd))