"""
   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 #survate 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.
    ebd = av*A - ao*A**(2./3) - aa*((A-2*Z)**2)/A - ac*Z**2/A**(1./3) + 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}'.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)