!
! FUNCTIONAL DESCRIPTION:
!
! Calculates the decay-widths in MeV for light particle emission Z<3 from the
! Hauser-Feshbach formalism. The summation over both spin and orbital angular
! momentum is performed.
!
! DUMMY ARGUMENTS:
!
! input: mode (INTEGER ::)
!        = 1 FOR NEUTRON EMISSION
!        = 2 FOR PROTON
!        = 3     DEUTERON
!        = 4     TRITON
!        = 5     HELIUM-3
!        = 6     ALPHA PARTICLE
!        = 7     Lithium-6
!        = 8     Lithium-7
!
!      Z0 (REAL (KIND=4) ::) proton number of decaying nucleus
!      A0 (REAL (KIND=4) ::) nucleon number of decaying nucleus
!      Ex0 (REAL (KIND=4) ::) excitation energy (MeV) of decaying nucleus
!      J0 (REAL (KIND=4) ::) spin of decaying nucleus
!
! IMPLICIT INPUTS:
!
!    none
!
! IMPLICIT OUTPUTS:
!
! I_Channels - number of open light particle channels
! gamma(1:I_channels) - decay widths for each light particle channel
! Za(1:_I_channels) - proton number of light fragment
! Aa(1:I_channels)  - nucleon number of light fragment
! Ex_residue(1:I_channels) - excitation energy of daughter nucleus

MODULE evap
USE gem_kind
USE mass_stuff
USE functs
USE yrast_mod
USE stat_mod
IMPLICIT NONE
PRIVATE 

PUBLIC :: evaporation_simple,prepare_evap


INTEGER, PRIVATE :: mode_gs,Iz2,Ia2,S0,s2,s2_best,jj
INTEGER, PRIVATE :: i
INTEGER, PRIVATE :: l,lms,l_min,s2_min


REAL (KIND=r4), PRIVATE :: Z2,A2,Z23,A23,const,Ec,pair_energy
REAL (KIND=r4), PRIVATE :: separation_c,r2,level_density
REAL (KIND=r4), PRIVATE :: mass2,mass2ld,shell_corr,ba,entropy
REAL (KIND=r4), PRIVATE :: E_rot,E_sad,J2,U2,ef0,aden
REAL (KIND=r4), PRIVATE :: ran,temp,ef,m_inertia_orbit,m_inertia_2
REAL (KIND=r4), PRIVATE :: fract,prob,prob2,ek,ranges
REAL (KIND=r4), PRIVATE :: m_inertia_saddle

INTEGER, PARAMETER, PRIVATE :: N_mode=15
REAL (KIND=r4), DIMENSION(n_mode), PARAMETER, PRIVATE :: Z1 = &
(/0.0,1.0,1.0,1.0,2.0,2.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0/)

REAL (KIND=r4), DIMENSION(N_mode), PARAMETER, PRIVATE :: A1 = &
(/1.0,1.0,2.0,3.0,3.0,4.0,6.0,6.0,6.0,6.0,6.0,6.0,7.0,7.0,7.0/)

REAL (KIND=r4), DIMENSION(N_mode), PARAMETER, PRIVATE :: S = &
(/0.5,0.5,1.0,0.5,0.5,0.0,1.0,3.0,0.0,2.0,2.0,1.0,1.5,0.5,3.5/)

REAL (KIND=r4), DIMENSION(N_mode), PARAMETER, PRIVATE :: DEF = &
(/8.07,7.28,13.13,14.949,14.93,2.424, &
14.0854,16.368,17.648,18.40,19.46,19.73,& !6Li
14.9067,15.384,19.537/)   !7Li

INTEGER, DIMENSION(N_mode), PARAMETER, PRIVATE :: level_no = &
(/0,0,0,0,0,0,0,1,2,3,4,5,0,1,2/)

! separation for Coulomb barrier
REAL (KIND=r4), DIMENSION(5,2:8), PARAMETER, PRIVATE :: sep_c = RESHAPE(&
(/ 4.509, 1.466, -0.5413, -3.218,   0.0,&   !proton
   6.723, 1.205, -0.5938, -0.2433,  0.0,&   !deuteron
   6.319, 1.346, -0.5901, -0.8494,  0.0,&   !triton
   5.812, 1.542, -0.8719,  0.7053,  0.0,&   !3He
   5.935, 1.308, -0.5129, -0.4928,  0.0,&   !alpha particle
   6.425, 1.543, -0.8504,  0.3892,  0.0,&   !6Li
   6.245, 1.536, -0.8019,  1.5555,  0.0/),& !7Li
(/5,7/))

!angular momentum dependence of Coulomb barrier, i.e. (l(l+1)) dependence
REAL (KIND=r4), DIMENSION(6,8), PARAMETER, PRIVATE :: vc0l = RESHAPE(&
(/-0.0111, -0.02063 , -0.8885, 1.508, 0.03129 , -0.02062 , & !neutron
  -0.0111, -0.02063 , -0.8885, 1.508, 0.03129 , -0.02062 , & !proton
  -0.0001, -469.3   , -2.555,  5.238, 0.1851  , -0.6560 ,  & !deuteron
  -0.0011, -24.61   , -3.773,  6.793, 0.2540  , -0.5061 ,  & !triton
  -0.0001, -376.1   , -3.886,  6.426, 0.1821  , -0.4884 ,  & !3He
  -0.0061, -0.6657  , -3.841,  6.913, 0.2497  , -0.2416 ,  & !alpha
   0.0021, -99.80   , -3.291,  5.537,-0.3208  ,  1.616721, & !6Li
   0.0021, -89.63   , -5.523,  8.092,-0.0324  ,  1.35290/),& !7Li
(/6,8/))

CONTAINS
SUBROUTINE evaporation_simple(mode,Z0,A0,Ex0,J0)

INTEGER, INTENT(IN) :: mode
REAL (KIND=r4), INTENT(IN) :: Z0,A0,EX0,J0

IF (mode > n_mode) RETURN

z2 = z0 - z1(mode)
a2 = a0 - A1(mode)
IF (z1(mode) > z2) RETURN
IF (A1(mode) > A2) RETURN
IF (z2 >= A2) RETURN
A23 = A2**0.3333333333
Z23 = Z2**0.3333333333
!-------------------------------------------------------------------------

SELECT CASE(mode)
CASE (-2:6)
   mode_gs = mode
CASE (7:12)
   mode_gs = 7
CASE (13:15)
   mode_gs = 8
END SELECT


CALL get_pairing(NINT(z2),NINT(A2),pair_energy)
if (z0 < Z_shell) pair_energy = 0 


IF (mode == 1) THEN
   Ec = 0.0
ELSE
! find radius for calculation of Coulomb barrier
   separation_c = sep_c(1,mode_gs) + sep_c(2,mode_gs)*A23&
  + sep_c(3,mode_gs)*Z23 + sep_c(4,mode_gs)/Z2 
   Ec = Z1(mode)*Z2*1.44/separation_c
END IF

!find radius for calculation of l dependence term of the coulomb barrier
const = 1.0/(vc0l(2,mode_gs)*(1.0-EXP(-vc0l(1,mode_gs)*A2)) &
   + vc0l(3,mode_gs)*Z23 + vc0l(4,mode_gs)*A23 &
   + vc0l(5,mode_gs)*z23**2.0 + vc0l(6,mode_gs)*A23**2.0 )


M_inertia_Orbit = 20.424/const
R2 = r0*A2**0.3333333
M_inertia_2 = 0.4*A2*R2**2
M_inertia_saddle = M_inertia_2 + M_inertia_orbit


! calculate the mass defect of the residual system
iz2 = z2
ia2 = a2
CALL getmass(iZ2,iA2,Mass2,mass2ld)
IF (para%mass_option /= 1 .AND. z0 > Z_shell) THEN
    shell_corr = mass2 - mass2ld - pair_energy
ELSE
    shell_corr = 0.0
END IF

! calculate the binding energy
ba = Mass2 + def(mode) - Mass0


IF (A0 < 8) THEN
  E_rot = J0**2.0*20.8993/M_inertia_saddle
  E_sad = Ba + E_rot + Ec
  J2 = J0
ELSE
  J2 = J0*M_inertia_2/M_inertia_saddle
  CALL yrast(Z2,A2,J2,0.0,0.0,E_rot)
  E_rot = E_rot + (J0-J2)**2.0*20.8993/M_inertia_orbit
  E_sad = Ba + E_rot + Ec
END IF
s2_best = INT(J2) + 1

U2 = Ex0 - E_sad
IF (U2 < 0) RETURN
jj = 0
s0 = INT(J0)


CALL get_little_a(A2,Z2,j2,U2,bs,shell_corr,E_rot,aden,entropy,&
pair_energy,1.0)

ef0 = entropy

IF (U2 > 40) THEN
    CALL light_particle_quick (mode,Ex0)
ELSE
    CALL light_particle (mode,Ex0)
END IF

IF (jj == 0) RETURN
IF (storage(jj)%gam < 0) RETURN
I_channels = I_channels + 1
Za(I_channels) = NINT(Z1(mode))
Aa(I_channels) = NINT(A1(mode))
gamma(I_channels) = storage(jj)%gam


dau(I_channels)%Ex_residue =  Ex0 - Ba - Ec  - REAL(l*(l+1))*const - ek 
dau(I_channels)%J_residue = s2
dau(I_channels)%s_particle = s(mode)
dau(I_channels)%level = level_no(Mode)
IF (para%I_angle) THEN
   dau(I_channels)%ll = l
   dau(I_channels)%e_kinetic = ek + ec + REAL(l*(l+1))*const
   dau(I_channels)%l_p_s = l
   dau(I_channels)%s_particle = 0
END IF
RETURN

END SUBROUTINE evaporation_simple

SUBROUTINE light_particle (mode,Ex0)

INTEGER, INTENT(IN) :: mode
REAL (KIND=r4), INTENT(IN) :: EX0
REAL (KIND=r4) :: eff


! loop over lms = l - s2
! l = orbital angular momentum of the emitted particle
! s2 = spin of residual
! lms is confined to limits -J0 to J0
lms = -s0
eff = ef0
s2_min = s0 + 1

DO
   IF (lms > s0) EXIT
   IF ((ef0-eff) > 4.0 .AND. s2_min < s2_best) EXIT
!  find minimum values of J2 and l which are consistent
!  with this value of lms
   l_min = INT(REAL(s0+lms+1)/2.0)
   s2_min = l_min - lms
   s2 = s2_min
   ef = ef0
! loop over spin of residual nucleus
   DO
       IF (ef0-ef > 4) EXIT
       l = s2 + lms      !orbital angular momentum
       CALL yrast(Z2,A2,REAL(s2),0.0,0.0,E_rot) ! rotational energy of residue
       U2 = Ex0 - Ba - Ec - E_rot - REAL(l*(l+1))*const
       if (u2 > 0)CALL get_little_a(A2,Z2,REAL(s2),U2,bs,shell_corr,E_rot,&
          aden,entropy,pair_energy,1.0)

       IF (U2 > 0.0) THEN
           ef = entropy
           temp = (U2/aden)**0.5
           IF (s2 == s2_min) eff = ef
           IF (ee-ef < 65.0 .AND. u2 > 0.0) THEN

               CALL get_level_density(U2,REAL(s2),aden,&
               M_inertia_2,entropy,level_density)

               jj = jj + 1
               IF (jj > n_light) THEN
                   WRITE (UNIT=*, FMT=*) "increase N_light in stat_mod.f90"
                   STOP 
               END IF
     
               storage(jj)%gam = level_density*temp*(2.0*s(mode)+1.0)
               storage(jj)%energy = REAL(l)
               storage(jj)%spin = s2
               IF (jj > 1) storage(jj)%gam = &
               storage(jj)%gam + storage(jj-1)%gam
           END IF
       ELSE
           ef = ef0 - 5.0
           IF (s2 == s2_min) eff = ef0 - 5.0
       END IF
       s2 = s2 + 1
   END DO
   lms = lms + 1
END DO
IF (jj == 0) RETURN
IF (storage(jj)%gam < 0.0) RETURN

i = 1

CALL RANDOM_NUMBER(prob)

DO 
   fract = storage(i)%gam/storage(jj)%gam
   IF (fract > prob) EXIT
   i = i + 1
END DO

l = NINT(storage(i)%energy)
s2 = storage(i)%spin
CALL yrast(Z2,A2,REAL(s2),0.0,0.0,E_rot)   ! rotational energy

U2 = Ex0 - Ba - Ec - E_rot - REAL(l*(l+1))*const

CALL get_little_a(A2,Z2,REAL(s2),U2,bs,shell_corr,E_rot,aden,entropy,&
pair_energy,1.0)
temp = (U2/aden)**0.5


ef = entropy
prob2 = 1.0
prob = 0.0
Ranges = MIN(4.0*temp,U2)

DO
   IF (prob2 < prob) EXIT
   CALL RANDOM_NUMBER(ran)
   ek = ranges*RAN
   prob = EXP(2.0*(aden*(U2-ek))**0.5 - ef)&
 *(1.0+U2**2.0)/(1.0+(U2-ek)**2.0)
   CALL RANDOM_NUMBER(prob2)
END DO

RETURN
END SUBROUTINE light_particle

SUBROUTINE light_particle_quick (mode,Ex0)

INTEGER, INTENT(IN) :: mode
REAL (KIND=r4), INTENT(IN) :: EX0


REAL (KIND=r4) :: t
REAL (KIND=r4) :: fact1,fact2,fact3,bit,x_min,integrate,ss1,ss2
REAL (KIND=r4) :: M_inertia_star,bb(200),s2_max,ss,eg

M_inertia_star = 1.0/(1.0/M_inertia_orbit+1.0/M_inertia_2)



! loop over lms = l - s2
! l = orbital angular momentum of the emitted particle
! s2 = spin of residual
! lms is confined to limits -J0 to J0
lms = -s0
ef = ef0
s2_min = s0 + 1

DO
  IF (lms > s0) EXIT
  IF ((ef0-ef) > 4.0 .AND. s2_min < s2_best) EXIT
!  find minimum vales of J2 and l which are consistent
!  with this value of lms
   l_min = INT(REAL(s0+lms+1)/2.0)
   s2_min = l_min - lms
   CALL yrast(Z2,A2,REAL(s2_min),0.0,0.0,E_rot) ! rotational energy of residue
   U2 = Ex0 - Ba - Ec - E_rot - REAL(l_min*(l_min+1))*const
   IF (U2 > 0.01) THEN

      CALL get_little_a(A2,Z2,REAL(s2_min),U2,bs,shell_corr,e_rot,aden,entropy,&
      pair_energy,1.0)
      ef = entropy
      IF (ee-ef < 65.0 .AND. u2 > 0.0) THEN
          temp = (U2/aden)**0.5
          fact1 = (M_inertia_2*(REAL(l_min)+0.5)      &
         +M_inertia_orbit*(REAL(s2_min)+0.5))**2.0    &
         /M_inertia_2/M_inertia_orbit/M_inertia_saddle&
         /temp*20.424
          fact2 = (M_inertia_star*temp/20.424)**0.5
          bit = (REAL(lms)+0.5)*M_inertia_star       &
         /M_inertia_orbit + M_inertia_star/2.0/M_inertia_2
          x_min = (REAL(s2_min)-0.5+bit)/fact2

!         now we are going to approximate the integral
!         from x_min to infinity of EXP(-x**2.0) using
!         the approximation for the normal probability
!         from "Handbook of Mathematical functions",
!         by Abramowitz and Stegan, Page932, accurate
!         to 7.5e-8
!
          t = 1.0/(1.0+0.2316419*1.4142135*x_min)
          fact3 = t*0.31938 - t**2.0*0.356563782   &
          + t**3.0*1.781477937 - t**4.0*1.821255978 &
          + t**5.0*1.330274429
          fact3 = fact3*0.3989423*1.77245385

          integrate = EXP(fact1-x_min**2.0)*(fact2**2.0 &
         + fact2*(1.0-2.0*bit)*fact3)

!         check to see if the integration along s2 gives a
!         value greater than that for s2_min
          IF (integrate/(2.0*REAL(s2_min)+1.0) < 1.0) &
          integrate = 2.0*REAL(s2_min)+1.0
          jj = jj + 1

          IF (jj > n_light) THEN
              WRITE (UNIT=*, FMT=*) "increase N_light"
              STOP 
          END IF
    

          CALL get_level_density(U2,0.0,aden,M_inertia_2,&
          entropy,level_density)

          storage(jj)%gam = temp*(2.0*s(mode)+1.0)*integrate*level_density

          bb(jj) = bit
          storage(jj)%spin = lms
          IF (jj > 1) storage(jj)%gam = storage(jj)%gam + storage(jj-1)%gam
      END IF
   ELSE
      ef = 0.0
   END IF
   lms = lms + 1
END DO
IF (jj == 0) RETURN
IF (storage(jj)%gam < 0) RETURN

fract = 0
i = 0

!find the values of (l - s2) form the stored particle widths
CALL RANDOM_NUMBER(prob)
DO
  IF  (fract > prob) EXIT
  i = i + 1
  fract = storage(i)%gam/storage(jj)%gam
END DO
lms = storage(i)%spin

! now find the range over which we will choose s2 the residue's spin
l_min = INT(REAL(s0+lms+1)/2.0)
s2_min = l_min - lms
CALL yrast(Z2,A2,REAL(s2_min),0.0,0.0,E_rot)   ! rotational energy
U2 = Ex0 - Ba - Ec - E_rot - REAL(l_min*(l_min+1))*const

CALL get_little_a(A2,Z2,REAL(s2_min),U2,bs,shell_corr,E_rot,aden,entropy,pair_energy,1.0)

temp = (U2/aden)**0.5
fact2 = (M_inertia_star*temp/20.454)**0.5
s2_max = ((8.0*fact2**2.0+(2.0*bb(i)-1)**2.0)**0.5 &
       -(2.0*bb(i)+1.0))/4.0
s2_max = MAX(REAL(s2_min)-.5,s2_max)
ef = ((s2_max+bb(i))/fact2)**2.0
eg = 2.0*s2_max+1.0
! find range over which to pick s2
ss = s2_max
prob = 1.0
DO
  IF (prob < 0.05) EXIT
  ss = ss + fact2/2.0
  prob = (2.0*ss+1.0)/eg*EXP(ef-((ss+bb(i))/fact2)**2.0)
END DO
ss2 = ss
ss1 = REAL(s2_min) - .5

! now select the final spin of the residue and the angular momentum
! of the particle with the provision that the final thermal excitation
! energy of the residue is positive
DO
   DO  
       CALL RANDOM_NUMBER(ran)
       ss = ss1 + (ss2-ss1)*RAN
       prob = (2.0*ss+1.0)/eg*EXP(ef-((ss+bb(i))/fact2)**2.0)
       CALL RANDOM_NUMBER(prob2)
       IF (prob2 < prob) EXIT
   END DO
   s2 = NINT(ss)
   l = s2 + lms
   CALL yrast(Z2,A2,REAL(s2),0.0,0.0,E_rot)   ! rotational energy
   U2 = Ex0 - Ba - Ec - E_rot - REAL(l*(l+1))*const
   IF (U2 >= 0.0) EXIT
END DO

!now choose the kinetic energy of the emitted particle

CALL get_little_a(A2,Z2,REAL(s2),U2,bs,shell_corr,e_rot,aden,entropy,&
pair_energy,1.0)

temp = (U2/aden)**0.5
ef = entropy
Ranges = MIN(4.0*temp,U2)

DO 
  CALL RANDOM_NUMBER(ran)
  ek = ranges*RAN
  prob = EXP(2.0*(aden*(U2-ek))**0.5 - ef) &
 *(1.0+U2**2.0)/(1.0+(U2-ek)**2.0)
  CALL RANDOM_NUMBER(prob2)
  IF (prob2 < prob) EXIT
END DO


RETURN
END SUBROUTINE light_particle_quick


SUBROUTINE prepare_evap ()
INTEGER :: i

i = 1
DO
   if (i > N_mode) EXIT
   if (z1(i) >= para%z_imf_min) EXIT 
   list_hf0(i) = i
   i = i + 1
END DO
i = i + 1
list_hf0(i) = -1

END SUBROUTINE prepare_evap
END MODULE evap