MODULE mass_stuff
USE gem_kind
IMPLICIT NONE
PRIVATE

PUBLIC :: getmass, prepare_mass, prepare_mass_tf, massfr, droplet_mass
PUBLIC :: get_pairing


INTEGER, PARAMETER, PRIVATE :: n_mass = 4200
INTEGER, DIMENSION(0:100), PRIVATE :: chz
REAL (KIND=r4), DIMENSION(n_mass), PRIVATE :: frmass,exmass,beta2

INTEGER, DIMENSION(0:100), PARAMETER, PUBLIC :: cha1 = (/ &
1,  1,  3,  5,  6,  8,  9, 11, 13, 14, 16, 18, 19, 22, 22,  &
24, 25, 27, 29, 30, 32, 34, 35, 37, 38, 40, 42, 44, 46, 48, &
51, 53, 55, 57, 59, 61, 63, 66, 68, 70, 72, 74, 77, 79, 81, &
83, 86, 88, 90, 92, 94, 97, 99, 101,103,106,108,110,113,115,&
118,120,123,125,128,130,133,136,138,141,143,146,149,151,154,&
156,159,162,165,167,170,173,175,179,183,184,190,193,196,199,&
203,206,206,207,212,222,226,233,230,231,232/)

INTEGER, DIMENSION(0:100), PARAMETER, PUBLIC :: cha2 = (/  &
1,  3,  7, 10, 14, 17, 19, 21, 24, 27, 31, 34, 37, 39, 41,  &
44, 48, 51, 54, 57, 60, 64, 67, 70, 73, 76, 79, 82, 85, 88, &
91, 94, 97,100,103,106,109,111,114,117,120,123,125,128,131, &
133,136,139,141,144,146,149,151,154,156,159,161,164,166,168,&
171,173,175,177,180,182,184,186,188,191,193,195,197,199,201,&
203,205,207,209,211,213,210,214,215,218,219,226,229,230,232,&
236,238,240,241,246,247,251,251,255,256,258/)

CONTAINS

SUBROUTINE getmass(iZ,iA,mass_exp,mass_ld,beta_two)
! 
! FUNCTIONAL DESCRIPTION:       
! 
! This subroutine looks up the mass excess (MeV) of a nucleus (Z,A).
! 
! DUMMY ARGUMENTS:
! 
! input  Iz: INTEGER :: - proton number of nucleus
!        Ia: INTEGER :: - nucleon number of nucleus
!        type: logical =1 mass defect from Yakawa+exponential model
!                           =0 experimental value
!
! output mass: REAL (KIND=r4) ::   mass_exp = mass defect with pairing and d
!                                  shell correction(MeV)
!                                  mass_ld = mass defect from liquid drop
!                                  term only
!      only experimental mass for Z < 11 are included in table at
!      present
! 

INTEGER, INTENT(IN) :: iZ,iA
REAL (KIND=r4), INTENT(OUT), OPTIONAL :: mass_ld,mass_exp,beta_two


INTEGER :: ind
REAL (KIND=r4) :: m_ld,m_exp,beta22

beta22 = 0.
    
IF (iZ > 100) THEN
   CALL massfr(REAL(iZ),REAL(iA),m_ld)
   m_exp = m_ld
ELSE IF (iA >= cha1(iZ) .AND. iA <= cha2(iZ)) THEN
   ind = chz(iz) - cha1(iz) + iA
   if (ind > n_mass)THEN
       WRITE (UNIT=*, FMT=*) "increase n_mass in mass.f90"
       STOP
   END IF 

   m_ld = frmass(ind)
   beta22 = beta2(ind)

   IF (iZ <= 100) THEN
       m_Exp = exmass(ind)
   ELSE
       WRITE (UNIT=*, FMT=*) "No exp mass for Z,A =",iZ,iA
       m_Exp = frmass(ind)
   END IF
ELSE
   CALL massfr(REAL(iZ),REAL(iA),m_ld)
   m_exp = m_ld
END IF

IF (PRESENT(mass_ld)) mass_ld = m_ld
IF (PRESENT(mass_exp)) mass_exp = m_exp
IF (PRESENT(beta_two)) beta_two = beta22

RETURN
END SUBROUTINE getmass



SUBROUTINE prepare_mass ()
INTEGER :: izz,I_READ_ERROR,ind
INTEGER :: Inn,iaa
REAL (KIND=r4) :: x1,x2,x3,x4,x5,x6,x7,x8,x9,m1,m2,m3,m4,m5,m6
REAL (KIND=r4) :: pairing

! the moller table starts at Z=8, to get
! lighter elements, use the thomas-fermi table
CALL prepare_mass_tf ()

chz(0) = 1
DO izz=1,100
    chz(izz) = chz(izz-1) + cha2(izz-1) - cha1(Izz-1) + 1
END DO


OPEN (                &
 UNIT=56,             &
 FILE = gemini_home//"mass.tbl",   &
 STATUS = "OLD",      &
 ACTION= "READ")
WRITE (*,*) "opened mass.tbl"

DO 
    READ (UNIT=56, FMT="(3I5,9F10.3,6F10.2)", IOSTAT=I_READ_ERROR)&
   IZZ,Inn,IAA,X1,X2,X3,X4,X5,X6,X7,X8,X9,M1,M2,M3,M4,M5,M6
    
    IF (I_READ_ERROR /= 0) EXIT
    IF (izz > 100) EXIT
    IF (iAA >= cha1(iZZ) .AND. iAA <= cha2(iZZ)) THEN
        ind = chz(izz) - cha1(izz) + iAA

        IF (ind > n_mass)THEN
             WRITE (UNIT=*, FMT=*) "increase n_mass in mass.f90"
             STOP
        END IF 

        CALL get_pairing(izz,iaa,pairing)
             

        frmass(ind) = m2 - m1 - pairing
           
        exmass(ind) = m2

        beta2(ind) = x6

!       IF (m3 /= 0) THEN
!            exmass(ind) = m3
!       ELSE
!            exmass(ind) = m2
!       END IF

    END IF 
END DO
CLOSE (UNIT=56)
END SUBROUTINE prepare_mass

SUBROUTINE prepare_mass_tf ()
INTEGER :: i, i_READ_Error,ind
INTEGER :: Inn,izz,iaa
REAL (KIND=r4) :: x1,x2,x3,x4,x5,m1,m2
REAL (KIND=r4) :: shell,pairing
LOGICAL :: be

INQUIRE (FILE = gemini_home//"mass_tf.gdat", &
 EXIST = be)

!open file and read plm2 numbers
if (be) THEN
   OPEN (             &
   UNIT=2,            &
   FILE = gemini_home//"mass_tf.gdat", &
   STATUS = "OLD",    &
   ACTION = "READ",    &
   FORM = "UNFORMATTED")   

   READ (UNIT=2) chz,frmass,exmass
   CLOSE (UNIT=2)
   WRITE (*,*) "opened mass_tf.gdat"
   RETURN
END IF 


chz(0) = 1
DO izz=1,100
    chz(izz) = chz(izz-1) + cha2(izz-1) - cha1(Izz-1) + 1
END DO


OPEN (                    &
 UNIT=56,                 &
 FILE = gemini_home//"mass_tf.tbl",    &
 STATUS = "OLD",          &
 ACTION="READ")
WRITE (*,*) "opened mass_tf.tbl"
! skip introductory remarks in file
DO i=1,66
    READ (UNIT=56, FMT=*)
END DO

DO 
    READ (UNIT=56,                                                &
    FMT="(I3,I4,I4,F5.1,F5.1,F8.2,F6.2,F5.2,F6.2,F8.2,F7.2,F7.2)" &
    ,IOSTAT=I_READ_ERROR)                                         &
    IZZ,Inn,IAA,x1,x2,x3,shell,pairing,x4,x5,M1,m2              

    IF (I_READ_ERROR /= 0) EXIT
    IF (izz > 100) EXIT
    IF (iAA >= cha1(iZZ) .AND. iAA <= cha2(iZZ)) THEN
       ind = chz(izz) - cha1(izz) + iAA

       IF (ind > n_mass)THEN
           WRITE (UNIT=*, FMT=*) "increase n_mass in mass.f90"
           STOP
       END IF 

       CALL get_pairing(izz,iaa,pairing)
 

       frmass(ind) = m1 - shell - pairing

!      if you want to backshift level densities by the
!      congruence then uncomment the following line
!      and the liune in the GET_PAIRING subroutine      
!      frmass(ind) = m1 - shell - pairing - x4

        exmass(ind) = m1

!       IF (m2 /= 0) THEN
!           exmass(ind) = m2
!       ELSE
!           exmass(ind) = m1
!       END IF
    END IF 
END DO
CLOSE (UNIT=56)

!having got all the masses, put them into a file
OPEN (                 &
UNIT=2,                &
FILE = gemini_home//"mass_tf.gdat", &
STATUS = "NEW",        &
ACTION = "WRITE",      &
FORM = "UNFORMATTED")
 
WRITE (UNIT=2) chz,frmass,exmass
CLOSE (UNIT=2)
WRITE (*,*) "created mass_tf.gdat"
END SUBROUTINE prepare_mass_tf

PURE SUBROUTINE massfr(Z,AA,mass)
       
! 
! FUNCTIONAL DESCRIPTION:       
! 
! Calculates macroscopic finite range model masses of spherical
! nucleus using formula of Krappe, Nix, and Sierk 
! Reference- (Phys Rev C20, 992 (1979))
! modified to use the parameters of Moller + Nix Nucl. Phys. A361(1981)
! 117. Pairing correction term for odd-odd nuclei
! is included, as this is the most appropriate ground state for hot nuclei
! where shell and pairing effects have washed out.
! 
! DUMMY ARGUMENTS:
! 
! input       Z : REAL (KIND=r4) :: number of protons
!             N : REAL (KIND=r4) :: number of neutrons
!     
! output  mass : REAL (KIND=r4) :: mass defect in MeV
! 
! 
REAL (KIND=r4), INTENT(IN) :: Z,AA
REAL (KIND=r4), INTENT(OUT) :: mass

REAL (KIND=r4) :: N,A13,A23,Fiss,Enz,Evol,x,fact,Es,Ec,Ewig
REAL (KIND=r4) :: ad,Ecd,akf,Ecpf,Ea0,Eca,deltau,deltal,Ep,I

REAL (KIND=r4), PARAMETER :: mn = 8.071431
REAL (KIND=r4), PARAMETER :: mh = 7.289034
REAL (KIND=r4), PARAMETER :: e2 = 1.4399764
REAL (KIND=r4), PARAMETER :: b = 0.99
REAL (KIND=r4), PARAMETER :: w = 36.
REAL (KIND=r4), PARAMETER :: ael = 1.433e-5
REAL (KIND=r4), PARAMETER :: r0 = 1.16
REAL (KIND=r4), PARAMETER :: a = 0.68
REAL (KIND=r4), PARAMETER :: as = 21.13
REAL (KIND=r4), PARAMETER :: ks = 2.3
REAL (KIND=r4), PARAMETER :: av = 15.9937
REAL (KIND=r4), PARAMETER :: kv = 1.927
REAL (KIND=r4), PARAMETER :: rp = 0.8
REAL (KIND=r4), PARAMETER :: c0 = 4.4
REAL (KIND=r4), PARAMETER :: ca = 0.212

N = AA - Z
A13 = AA**(1./3.)
A23 = AA**(2./3.)

! relative neutron excess
I = (N-Z)/AA
fiss = Z**2./A13

! neutron-proton terms
Enz = mn*N + mh*Z

! Volume energy 
Evol = -av*(1.-kv*I**2.)*AA

! Surface energy
x=a/(r0*A13)
fact=1.-3.*x**2. &
+ (1.+1./x)*(2.+3.*x+3.*x**2.)*exp(-2./x)
Es = as*(1.-ks*I**2.)*A23*fact

! Coulomb energy
fact = fiss-0.76361*Z**1.3333/A13
Ec = 0.6*e2/r0*fact

! Wigner term
Ewig = w*abs(I)-ael*Z**2.39

! correction to Coulomb energy for diffuse surface
! see Davies & Nix Phys. Rev. C14 (1976) 1977
ad=0.7071*b
x = ad/(r0*A13)
fact= 1 - 1.875*x+2.625*x**3. &
-.75*exp(-2./x)*(1.+4.5*x+7.*x**2.+3.5*x**3.)
Ecd = -3.*Z**2.*e2*ad**2./(r0*a13)**3.*fact

! correction to coulomb energy due to proton form factor
akf=1./r0*(7.06858*z/AA)**.33333
x=(rp*akf)**2.
Ecpf=-0.125*rp**2.*e2/r0**3. &
*(3.0208-0.113541667*x+0.0012624*x**2.)*Z**2./AA

!A0 term
Ea0=c0

!Charge asymmetry term
Eca=ca*(Z**2.-N**2.)/AA

! pairing term for odd-odd nuclei
deltau = 12./sqrt(AA)
deltal = 20./AA
Ep =deltau - 0.5 * deltal

! add all terms
mass = Enz+Evol+Es+Ec+Ewig+Ecd+Ecpf+Ea0+Eca+Ep

RETURN
END SUBROUTINE massfr

PURE SUBROUTINE droplet_mass(Z,A,mass)
 
!+
! 
! gives the Droplet model masses from the work of W.D. Myers
! "Droplet model of atomic nuclei" (IFI, Plenum, New York, 1977)
!-
REAL (KIND=r4), INTENT(IN) :: Z,A
REAL (KIND=r4), INTENT(OUT) :: mass

REAL (KIND=r4) :: N,delta,epsilons,A3,A23,A43
REAL (KIND=r4) :: Ev,Es,Ec

REAL (KIND=r4), PARAMETER :: a1=15.96
REAL (KIND=r4), PARAMETER :: a2=20.69
REAL (KIND=r4), PARAMETER :: J=36.8
REAL (KIND=r4), PARAMETER :: c1=0.73219
REAL (KIND=r4), PARAMETER :: Q=17.
REAL (KIND=r4), PARAMETER :: K=240.
REAL (KIND=r4), PARAMETER :: L=100.
REAL (KIND=r4), PARAMETER :: c2=0.00016302
REAL (KIND=r4), PARAMETER :: c3=1.28846
REAL (KIND=r4), PARAMETER :: c4=.55911
REAL (KIND=r4), PARAMETER :: c5=.00049274

N = A - Z
A3 = A**0.33333333
A23 = A**0.6666666
A43 = A**1.3333333
delta = ((N-Z)/A+3./16.*c1/Q*Z/A23)/(1.+9./4.*J/Q/A3)
epsilons = (-2*a2/A3+L*delta**2.+c1*Z**2./A43)/K
Ev = (-a1+j*delta**2.-0.5*K*epsilons**2.)*A
Es = (a2+9./4.*J**2./Q*delta**2.)*A23 
Ec =  c1*Z**2./A3 - c2*Z**2.*A3 - c3*Z**2./A &
 - c4/2**0.33333*Z - c5*Z**2.
mass = Ev + Es + Ec

! Wigner term
IF (NINT(n) == NINT(Z) .AND. MODULO(INT(N),2) == 1) THEN
    mass = mass + 30./A
ELSE
    mass = mass + 30.* ABS(N-Z)/A
END IF

! odd-even
! IF (MODULO(INT(Z),2) == 0 .AND. MODULO(INT(N),2) == 0) THEN
!   mass = mass - 12./A**0.5 +10/A
!ELSE IF (MODULO(INT(Z),2) == 1 .AND. MODULO(INT(N),2) == 1) THEN
!   mass = mass  + 12./A**0.5 - 10/A
!ELSE
!   mass = mass + 10/A
!END IF

mass = mass + N*8.07169 + Z*7.28922

! remove atomic binding
! mass = mass - 0.00001433*Z**2.39
RETURN
END SUBROUTINE droplet_mass


PURE SUBROUTINE get_pairing(iz,ia,pairing)
  
![subroutine_header_comments]
INTEGER, INTENT(IN) :: iz,ia
REAL (KIND=r4), INTENT(OUT) :: pairing

INTEGER :: inn,oez,oen

inn = ia - iz
oez = MODULO(iz,2)
oen = MODULO(inn,2)
if (inn == iz .AND. oez == 1) THEN
   pairing = 4.8/REAL(inn)**0.33333333 &
   + 4.8/REAL(iz)**0.33333333 &
   - 6.6/REAL(ia)**0.6666666 + 30/REAL(Ia)
ELSE IF (oez == 1 .AND. oen == 1) THEN
   pairing = 4.8/REAL(inn)**0.33333333 &
   + 4.8/REAL(iz)**0.33333333         &
   - 6.6/REAL(ia)**0.6666666
ELSE IF (oez == 1 .AND. oen == 0) THEN
   pairing = 4.8/REAL(iz)**0.33333333
ELSE IF (oez == 0 .AND. oen == 1) THEN
   pairing = 4.8/REAL(inn)**0.33333333
ELSE
   pairing = 0                                       
END IF 

! zero pairing for odd-odd
! pairing = pairing - 4.8/REAL(inn)**0.33333333 &
! - 4.8/REAL(iz)**0.33333333                   & 
! + 6.6/REAL(ia)**0.6666666

! congruence energy
! if you want to backshift by the congruence energy then
! uncomment the follow line
! pairing = pairing  - 10.*EXP(-4.2*REAL(ABS(iz-inn))/REAL(ia))
  
RETURN 
END SUBROUTINE get_pairing

END MODULE mass_stuff