MODULE yrast_mod
 
USE gem_kind
USE sierk_stuff
USE mass_stuff

IMPLICIT NONE
PRIVATE

PUBLIC :: sadfit, read_sad, yrastld, yrast, saddle, def_rot
PRIVATE :: cubic

REAL (KIND=r4), ALLOCATABLE, DIMENSION(:,:,:,:,:), PRIVATE :: C
REAL (KIND=r4), DIMENSION(150), PRIVATE :: sad_array
REAL (KIND=r4), PRIVATE :: Bs,Bc
REAL (KIND=r4), PRIVATE :: br_perp,br_para

 
CONTAINS
SUBROUTINE sadfit(Z0,A0,J0)
!
! FUNCTIONAL DESCRIPTION:
!
! this subroutine returns the mass-asymmetric saddle
! point energy for the nucleus z0,a0 with angular
! momentum j0. the saddle point energy is obtained from
! a fits to saddle pint energies calculated for two
! spheroids with a yukawa-plus-exponential nuclear force
! of Krappe, Sierk, AND Nix Phys. Rev C20 (1979)992
! using the parameters of Moller and Nix (Nucl. Phy.
! A361 (1981)117). the saddle point energies were
! calculated with the program sadnew of sierk for nuclei
! which lie along the line a=2.08*z+0.0029*z**2.  (see
! charity et al lbl22447) the program calculates the
! fissility parameter and y parameter in the finite range
! model. and uses these in a bi-cubic interpolation to
! obtain the saddle point energies. the barriers are
! then normalized so the symmetric barriers reproduce the
! the sierk barriers calculated by barfit which were
! obtained from a more realistic shape parameterization.
! Bob Charity     5 SEPT 1987
!
! DUMMY ARGUMENTS:
!
!input: Z0 (REAL (KIND=r4) ::)- proton number of decaying nucleus
!       A0 (REAL (KIND=r4) ::)- nucleon number of decaying nucleus
!       J0 (REAL (KIND=r4) ::)- spin of decaying nucleus
!
! IMPLICIT INPUTS:
!
!    none
!
! IMPLICIT OUTPUTS:
!
! sad_array(1:N_array) - finite range model conditional saddle-point energies
!                      less the masses of the two fragments and their Coulomb
!                      energy
!
! SIDE EFFECTS:
!
!    none
!
!
REAL (KIND=r4), INTENT(IN) :: Z0,A0,J0



INTEGER :: ii,n_array
INTEGER :: Iy,Ix,Ia
INTEGER :: Ia1,Ia2,Iz1,Iz2
LOGICAL :: first = .TRUE.
REAL (KIND=r4) :: r1,r2,A1,A2,Z1,Z2
REAL (KIND=r4) :: E_coul
REAL (KIND=r4) :: ac,rz,cs,delcs,esz,emz,zsqoa,x,tz,elz,elzohb
REAL (KIND=r4) :: y,sad0,ess,bfis,sadf
REAL (KIND=r4) :: alpha,Dkx,Dky,Dka,Rx,Ry,Ra,k1,k2,k3,k4
REAL (KIND=r4) :: c0,c1,c2,c3
REAL (KIND=r4) :: mass1,mass2

REAL (KIND=r4), PARAMETER :: srznw=1.16
REAL (KIND=r4), PARAMETER :: asnw=21.13
REAL (KIND=r4), PARAMETER :: aknw=2.3
REAL (KIND=r4), PARAMETER :: hbarc=197.32858
REAL (KIND=r4), PARAMETER :: alfinv=137.035982
REAL (KIND=r4), PARAMETER :: bb=0.99
REAL (KIND=r4), PARAMETER :: um=931.5016
REAL (KIND=r4), PARAMETER :: elm=0.51104
REAL (KIND=r4), PARAMETER :: spdlt=2.99792458e23

REAL (KIND=r4), DIMENSION(6) :: ky = (/0.0,0.02,0.05,0.10,0.15,0.20/)
REAL (KIND=r4), DIMENSION(11) :: ka = (/0.0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1./)
REAL (KIND=r4), DIMENSION(8) :: kx = (/0.0,0.08299,0.1632,0.2435,0.3217,0.402&
     ,0.51182,0.617422/)
 
!/*----find fissility of nucleus---------------------------*/
N_array = INT(A0/2.0)
Ia = 0
ac = 0.6*hbarc/alfinv/srznw
rz = srznw*A0**.333333333333333
cs = asnw * (1.0-aknw*((A0-2.0*Z0)/A0)**2.0)
delcs = 45.0*hbarc/(8.0*srznw*alfinv)*&
((bb/(1.4142*srznw))**3.0)*(Z0/A0)**2.0
cs = cs + delcs
esz = cs*A0**(2.0/3.0)
emz = um*A0 - elm*Z0
zsqoa = Z0**2.0/A0
x = ac*zsqoa/(2.0*cs)

IF (x > 0.61743) THEN
   if (first) THEN
      WRITE (UNIT=*, FMT=*) "No barriers available for this nucleus"
      WRITE (UNIT=*, FMT=*) "Z=",Z0,"A=",A0
      WRITE (UNIT=*, FMT=*) "using scaled 194Hg barriers"
      first = .FALSE.
   END IF 
   x = 0.61743  
END IF

!/*----find y fissility parameter-------------------------*/
IF (esz == 0.0) THEN
   WRITE (UNIT=*, FMT=*) "esz=0 in sad, Z=",Z0,"A=",A0
   STOP
END IF

IF (emz/esz <= 0.0) THEN
   WRITE (UNIT=*, FMT=*) "square root of negative for a=",A0,"z=",Z0
   WRITE (UNIT=*, FMT=*) "emz=",emz,"esz=",esz
   STOP
END IF

tz = rz*sqrt(emz/esz)/spdlt
elz = esz*tz
elzohb = elz/6.582173e-22
y = 1.25*(J0/elzohb)**2.0

!/*----------find normalization factor-------------------*/
IF (Z0 >= 19.0) THEN
   sad0 = -5.5286e-2 + 190.03*x + 235.189*x**2.0 &
   - 2471.249*x**3.0 + 4266.905*x**4.0 - 2576.048*x**5.0
   CALL jmax_sierk(Z0,A0) !initialize sierk sub
   CALL bfis_sierk(0.0,bfis)
   sadf = bfis/sad0
ELSE
   ess = 2376*x - 6062.4*x**2.0 + 8418.3*x**3.0
   sadf = (esz/ess)**2.12
END IF

!/*-----select subroutines for cubic-spline interpolation--*/


!/*find nearest knots for interpolation */
Iy = 2
DO
   IF (y <= ky(Iy)) EXIT
   Iy = Iy +1
   IF (Iy > 5) EXIT
END DO
Iy = Iy - 1

ix = 2
DO
   IF (x <= kx(Ix)) EXIT
   Ix = Ix + 1
   IF (Ix > 6) EXIT
END DO
Ix = Ix - 1
! /* now y is between knots ky(iy-1) & ky(iy) and x is  */
! /* between knots kx(ix) & kx(ix+1)       */

Dkx = kx(Ix+1) - kx(Ix)
Rx = (x - kx(Ix))/Dkx
Dky = ky(Iy+1) - ky(Iy)
Ry = (y - ky(iy))/Dky

! /*--- fill array with saddle point energies-------*/
DO ii = N_array,1,-1
   alpha = (A0 - 2.0 * REAL(ii))/A0
   IF (alpha >= ka(Ia + 1)) THEN
       DO
           IF (alpha < ka(Ia +1)) EXIT
           Ia = Ia + 1
       END DO
       k1 = cubic( C(Iy,Ix,1,Ia,1), C(Iy,Ix+1,1,Ia,1),&
        C(Iy,Ix,2,Ia,1), C(Iy,Ix+1,2,Ia,1),Dkx,Rx)
       k2 = cubic( C(Iy,Ix,1,Ia+1,1), C(Iy,Ix+1,1,Ia+1,1),&
        C(Iy,Ix,2,Ia+1,1), C(Iy,Ix+1,2,Ia+1,1),Dkx,Rx)
       k3 = cubic( C(Iy,Ix,1,Ia,2), C(Iy,Ix+1,1,Ia,2),&
        C(Iy,Ix,2,Ia,2), C(Iy,Ix+1,2,Ia,2),Dkx,Rx)
       k4 = cubic( C(Iy,Ix,1,Ia+1,2), C(Iy,Ix+1,1,Ia+1,2),&
        C(Iy,Ix,2,Ia+1,2), C(Iy,Ix+1,2,Ia+1,2),Dkx,Rx)

       Dka = ka(Ia+1) -ka(Ia)

       C0 = k1
       C1 = Dka*k3
       C2 = 3.0*(k2-k1) - (k4+2.0*k3)*Dka
       C3 = 2.0*(k1-k2) + (k3+k4)*Dka

       IF (J0 > 0.1) THEN
           k1 = cubic( C(Iy+1,Ix,1,Ia,1), C(Iy+1,Ix+1,1,Ia,1),&
           C(Iy+1,Ix,2,Ia,1), C(Iy+1,Ix+1,2,Ia,1),Dkx,Rx)
           k2 = cubic( C(Iy+1,Ix,1,Ia+1,1),                   &
             C(Iy+1,Ix+1,1,Ia+1,1), C(Iy+1,Ix,2,Ia+1,1),     &
             C(Iy+1,Ix+1,2,Ia+1,1),Dkx,Rx)
           k3 = cubic( C(Iy+1,Ix,1,Ia,2), C(Iy+1,Ix+1,1,Ia,2),&
             C(Iy+1,Ix,2,Ia,2), C(Iy+1,Ix+1,2,Ia,2),Dkx,Rx)  
           k4 = cubic( C(Iy+1,Ix,1,Ia+1,2),                   &
             C(Iy+1,Ix+1,1,Ia+1,2), C(Iy+1,Ix,2,Ia+1,2),     &
             C(Iy+1,Ix+1,2,Ia+1,2),Dkx,Rx)                   
           C0 = (k1-C0)*Ry + C0
           C1 = (Dka*k3-C1)*Ry + C1
           C2 = (3.0*(k2-k1) - (k4+2.*k3)*Dka - C2)*Ry + C2
           C3 = (2.0*(k1-k2) + (k3+k4)*Dka - C3)*Ry + C3
           END IF
       END IF
   Ra = (alpha - ka(Ia))/Dka
   sad_array(ii) = (C0 + C1*Ra + C2*Ra**2.0 + C3*Ra**3.0)*sadf
END DO

DO ii = para%Z_imf_min,N_array
   A1 = REAL(ii)
   A2 = A0 - A1
   Z1 = A1/A0 * Z0
   Z2 = Z0 - Z1
   ia1 = A1
   ia2 = A2
   iz1 = Z1
   iZ2 = Z2
   CALL getmass(iz1,ia1,Mass_LD=mass1)
   CALL getmass(Iz1+1,ia1,mass_LD=x)
   Mass1 = Mass1 + (x-Mass1)*(Z1-REAL(iz1))
   CALL getmass(iz2,ia2,Mass_LD=mass2)
   CALL getmass(iz2+1,ia2,mass_LD=x)
   Mass2 = Mass2 + (x-Mass2)*(Z2-REAL(iz2))

   r1 = r0*A1**0.33333
   r2 = r0*A2**0.33333
   E_coul = Z1*Z2*1.44/(r1 + r2 + sep)
   sad_array(ii) = sad_array(ii) - Mass1 - Mass2 - E_coul
END DO
RETURN
END SUBROUTINE sadfit


SUBROUTINE read_sad ()
INTEGER :: i1,i2,i3,i4,i5
LOGICAL :: be

ALLOCATE (c(6,8,2,11,2))
INQUIRE (FILE = gemini_home//"sad.gdat", &
 EXIST = be)

if (be) THEN
    OPEN (                &
    UNIT=8,               &
    FILE = gemini_home//"sad.gdat",     &
    ACTION = "READ",      &
    STATUS = "OLD",       &
    ACCESS = "SEQUENTIAL",&
    FORM = "UNFORMATTED")
    READ (UNIT=8) c
    CLOSE (UNIT=8)
    WRITE (*,*) "opnened sad.gdat"
ELSE
    OPEN (          &
    UNIT=8,         &
    FILE= gemini_home//"sad.tbl", &  
    STATUS = "OLD", &
    ACTION = "READ")
    WRITE (*,*) "opened sad.tbl"
    READ (UNIT=8, FMT=*)
    DO i5=1,2
     DO i4=1,11   
        DO i3=1,2 
            READ (UNIT=8,FMT=*)
            do i2=1,8
                READ (UNIT=8, FMT=*) (c(i1,i2,i3,i4,i5),i1=1,6)
            END DO
        END DO
     END DO 
    END DO

    CLOSE (UNIT=8)

    OPEN (                &
    UNIT=8,               &
    FILE = gemini_home//"sad.gdat",     &
    ACTION = "WRITE",      &
    STATUS = "NEW",       &
    ACCESS = "SEQUENTIAL",&
    FORM = "UNFORMATTED")
    WRITE (UNIT=8) c
    CLOSE (UNIT=8)
    WRITE (*,*) "created sad.gdat"
END IF
END SUBROUTINE read_sad
 

SUBROUTINE yrastld (Z,A,l,E_rot)
!
! FUNCTIONAL DESCRIPTION:
!
! Finds, by interpolation, the ROTATING LIQUID DROP MODEL yrast line.
! Modified from routine of Plasil o.r.n.l
!
!
! DUMMY ARGUMENTS:
!
! inputs:
! Z REAL (KIND=r4) :: - proton number of nucleus
! A REAL (KIND=r4) :: - nucleon number of nucleus
! l REAL (KIND=r4) :: - spin of nucleus
! outputs:
! E_rot REAL (KIND=r4) :: - rotational energy of nucleus
!
! IMPLICIT INPUTS:
!
!    none
!
! IMPLICIT OUTPUTS:
!
!    none
!
!
! SIDE EFFECTS:
!
!    none
!
!
 

INTEGER :: ix,iy
REAL (KIND=r4), INTENT(IN) :: Z,A,l
REAL (KIND=r4), INTENT(OUT) :: E_rot
REAL (KIND=r4) :: N
REAL (KIND=r4) :: paren,eso,x,ero,y,cx,bx,dx,by,cy,dy,h1,h2
REAL (KIND=r4) :: Hf

REAL (KIND=r4), DIMENSION(6,11), PARAMETER :: x1h = RESHAPE((/&
0.0,0.0,0.0,0.0,0.0,0.0,-0.0057,-0.0058,-0.006,-0.0061,-0.0062,&
-0.0063,-0.0193,-0.0203,-0.0211,-0.022,-0.023,-0.0245,-0.0402, &
-0.0427,-0.0456,-0.0497,-0.054,-0.0616,-0.0755,-0.0812,-0.0899,&
-0.0988,-0.109,-0.12,-0.1273,-0.1356,-0.147,-0.1592,-0.1745,   &
-0.1897,-0.1755,-0.1986,-0.2128,-0.2296,-0.251,-0.26,-0.255,   &
-0.271,-0.291,-0.301,-0.327,-0.335,-0.354,-0.36,-0.365,-0.372, &
-0.403,-0.42,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/),(/6,11/))

REAL (KIND=r4), DIMENSION(6,11), PARAMETER :: x2h = RESHAPE((/ &
0.0,0.0,0.0,0.0,0.0,0.0,-0.0018,-0.0019,-0.00215,-0.0024,-0.0025,&
-0.003,-0.0063,-0.00705,-0.0076,-0.0083,-0.0091,-0.0095,-0.015,  &
-0.0158,-0.0166,-0.0192,-0.0217,-0.025,-0.0245,-0.0254,-0.029,   &
-0.0351,-0.0478,-0.0613,-0.0387,-0.0438,-0.0532,-0.0622,-0.0845, &
-0.0962,-0.0616,-0.0717,-0.0821,-0.0972,-0.1123,-0.1274,-0.0793, &
-0.1014,-0.1138,-0.1262,-0.1394,-0.1526,-0.12,-0.134,-0.1503,    &
-0.1666,-0.1829,-0.1992,-0.1528,-0.171,-0.1907,-0.2104,-0.2301,  &
-0.2498,0.0,0.0,0.0,0.0,0.0,0.0/),(/6,11/))

REAL (KIND=r4), DIMENSION(10,20), PARAMETER :: x3h = RESHAPE((/          &
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-0.00012,-0.00014,-0.00016,      &
-0.00018,-0.0002,-0.00024,-0.00029,-0.00036,-0.00065,-0.00089,-0.00047,  &
-0.0005,-.00058,-.00065,-.00074,-.00085,-.00101,-.00124,-.00138,-.00178, &
-0.001,-0.00105,-0.00124,-0.00138,-0.00156,-0.00179,-0.00275,-0.00292,   &
-0.003,-0.003,-0.00176,-0.0019,-0.00211,-0.00235,-0.00263,-0.00298,      &
-0.00449,-0.0053,-0.0053,-0.0053,-0.003,-0.00308,-0.00318,-0.00352,      &
-0.00392,-0.00417,-0.0062,-0.0062,-0.0062,-0.0062,-0.00374,-0.0041,      &
-0.00444,-0.00488,-0.00521,-0.00545,-0.0066,-0.0066,-0.0066,-0.0066,     &
-0.0053,-0.0055,-0.00585,-0.0064,-0.00695,-0.007,-0.007,-0.007,-0.007,   &
-0.007,-0.00632,-0.007,-0.00742,-0.00792,-0.00856,-0.009,-0.009,-0.009,  &
-0.009,-0.009,-0.0079,-0.0085,-0.01022,-0.0119,-0.012,-0.012,-0.012,     &
-0.012,-0.012,-0.012,-0.00944,-0.0102,-0.0142,-0.0182,-0.019,-0.019,     &
-0.019,-0.019,-0.019,-0.019,-0.0112,-0.0133,-0.0182,-0.0238,-0.024,      &
-0.024,-0.024,-0.024,-0.024,-0.024,-0.01303,-0.0178,-0.0226,-0.0274,     &
-0.028,-0.028,-0.028,-0.028,-0.028,-0.028,-0.0165,-0.0254,-0.0343,       &
-0.0343,-0.034,-0.034,-0.034,-0.034,-0.034,-0.034,-0.0203,-0.033,-0.04,  &
-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.025,-0.0406,-0.046,-0.047,  &
-0.047,-0.047,-0.047,-0.047,-0.047,-0.047,-0.03036,-0.0482,-0.048,       &
-0.048,-0.048,-0.048,-0.048,-0.048,-0.048,-0.048,-0.0363,-0.0558,        &
-0.056,-0.056,-0.056,-0.056,-0.056,-0.056,-0.056,-0.056,-0.04234,        &
-0.0634,-0.064,-0.064,-0.064,-0.064,-0.064,-0.064,-0.064,-0.064,0.0,     &
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/),(/10,20/))

n = a - z
paren = 1.0 - 1.7826*((n-z)/a)**2
if (paren <= 0.) paren = 0.1  ! what to do here
eso = 17.9439*paren*a**.66667
x = 0.019655*z*(z/a)/paren
ero = 34.548*l**2/a**1.6667
y = 1.9254*l**2/(paren*a**2.3333)

ix = 20.*x + 1.0
cx = ix
bx = 20.*x + 1.0
dx = bx - cx

IF (x <= 0.25) THEN
  by = 10.*y + 1.
  IF (by > 9.0) by = 9.0
  IF (by < 1.0) by = 1.0
  iy = by
  cy = iy
  dy = by - cy
  h1 = (x1h(ix+1,iy)-x1h(ix,iy))*dx+x1h(ix,iy)
  h2 = (x1h(ix+1,iy+1)-x1h(ix,iy+1))*dx+x1h(ix,iy+1)
  hf = (h2-h1)*dy+h1
ELSE IF (x < 0.5) THEN
  by = 20.0*y + 1.0
  IF (by > 11.0) by = 10.0
  IF (by < 1.0) by = 1.0
  ix = ix - 5
  iy = by
  cy = iy
  dy = by - cy
  h1 = (x2h(ix+1,iy)-x2h(ix,iy))*dx+x2h(ix,iy)
  h2 = (x2h(ix+1,iy+1)-x2h(ix,iy+1))*dx+x2h(ix,iy+1)
  hf = (h2-h1)*dy + h1
ELSE
  IF (x > 0.94999999) x=0.949999999
  ix = 20.0*x + 1.0
  ix = ix - 10
  by = 100.0*y + 1.0
  IF (by > 19.0) by = 19.0
  IF (by < 1.0) by = 1.0
  iy = by
  cy = iy
  dy = by - cy
  h1 = (x3h(ix+1,iy)-x3h(ix,iy))*dx+x3h(ix,iy)
  h2 = (x3h(ix+1,iy+1)-x3h(ix,iy+1))*dx+x3h(ix,iy+1)
  hf = (h2-h1)*dy+h1
  END IF
E_rot =  ero + hf*eso
RETURN
END SUBROUTINE yrastld

SUBROUTINE yrast(Z,A,J,K,ratio,E_rot)
       
!+
! 
! FUNCTIONAL DESCRIPTION: 
! 
! There are two options which are hard wired into the code for when
! you have an angular momentum greater then the value at which
! the fission barrier vanishes. (jmax)
! if kink=.true. the moment-of-inertia at jmax is used to extrapolate
! up from the E_rot value at jmax. As Sierks momentum of inertia
! changes rapidly in the last few hbar, this option produces a kink
! in the "yrast line". Alternatively, if kink=.false. one can
! start the extrapolation at jmax-dj, in a region where the
! moment-of-inertia is more constant and so the "yrast line"
! has no kink and parabolic (nearly) for all angular momenta.
! These yrast lines are only useful for evaporation, as
! equilibrium fission is not defined as the fission barrier
! is zero. If one has a fission delay turned on, then the
! code will ignore fission, otherwise the fission yield
! predicted by gemini is not well defined.
! The kink=.true. option may be relavent if one assumed
! the object evaporating particles is deformed or di-nuclear
! formed in an initially  fast-fission or deep-inelastic
! reaction, the evaporation of light particles during a fission delay
! may remove enough angular momentum for the fission competition to
! be treated statistically after the delay.
! The kink=.FAlse. option may be more applicable for a asymmetric
! entrance channel where a more compact (spherical) object is formed.
! again if enough angular momentum is removed during the delay,
! the fission can be considered statistically.
!  When the deformed delay option is used, the yrast line is the
!  rotation plus deformation energy of a spheroid, ratio is the
!  ratio of the symmetry axis to the other axis.
!-
REAL (KIND=r4), INTENT(IN) :: Z,A,J,K,ratio
REAL (KIND=r4), INTENT(OUT) :: E_rot

REAL (KIND=r4) :: mom_inertia_1
REAL (KIND=r4) :: E_def,E_yrast,jmax
REAL (KIND=r4), PARAMETER :: dj = 3
LOGICAL, SAVE :: kink = .true.


IF (ratio > 0.0)THEN
!  rotational and deformation energy for rotating spheroid
!  taken from rotating liquid drop model.
!   The Bs and Bc quantities were calculated with formulas from
!   R. Nix"s thesis.
    IF (ratio /= 1.0) THEN
        E_def = (Bs-1.0)*17.9439* &
       (1.0-1.7826*((A-2.0*Z)/A)**2.0)*A**0.66666666 &
        + (bc-1.0)*0.7053*Z**2.0/A**0.33333333

    ELSE
        E_def = 0.0
    END IF
    E_rot = (Br_perp*(J*(J+1.0)-K**2.0)+Br_para*K**2.0)&
      *34.540/A**(5.0/3.0)

    E_yrast = E_rot + e_def
    E_rot = E_yrast
    RETURN
END IF

IF (z < 19) THEN !rotating liquid drop yrast line
    CALL yrastld(Z,A,J,E_rot)
    RETURN
ELSE   !rotating finite range model
    CALL jmax_sierk(Z,A,Jmax) !sierk's yrast line
    IF (J <= Jmax-dj) THEN
        CALL yrast_sierk(j,E_rot)
    ELSE IF (kink .AND. J <= jmax) THEN
        CALL yrast_sierk(j,E_rot)
    ELSE  !there is no well defined yrast line here, but we approx one.
          !by using the moment of inertia at Jmax to extrapolate up.
        CALL m_inertia_sierk(jmax,mom_inertia_1)

        IF (kink) THEN
            CALL yrast_sierk(jmax,E_rot)
            mom_inertia_1 = mom_inertia_1*0.4*r0ld**2.0*A**1.6666666666
            E_rot = E_rot + (J**2.0-jmax**2.0)/(2.0*Mom_inertia_1)*40.848
        ELSE
!           !using the moment of inertia at jmax to extrapolate up,
!           !produces a kink in the yrast line, a less kinked yrast
!           !line can be obtained using the reduced moment of inertia
!           !in the following line - of course it is not clear what the
!           correct yrast line is for such large angular momentum       
            CALL yrast_sierk(jmax-dj,E_rot)
            mom_inertia_1 = (5.0 + mom_inertia_1)/6.0
            mom_inertia_1 = mom_inertia_1*0.4*r0ld**2.0*A**1.6666666666
            E_rot = E_rot + (J**2.0-(jmax-dj)**2.0)/(2.0*Mom_inertia_1)*40.848
        END IF
    END IF
END IF
RETURN
END SUBROUTINE yrast

    PURE FUNCTION cubic(a,b,c,d,e,f) RESULT(funct_out)
    !
    ! FUNCTIONAL DESCRIPTION:
    !
    !    Cubic equation for cubic spline interpolation used in subroutine SADFIT
    !
    ! DUMMY ARGUMENTS:
    !
    !    input: a,b,c,d,e (REAL (KIND=4) ::) coefficients
    !           f (REAL (KIND=4) ::) value at which cubic polynomial is evaluated
    !
    REAL (KIND=r4), INTENT(IN) :: a,b,c,d,e,f
    REAL (KIND=r4) :: funct_out

    funct_out = a + f*(e*c+f*(3.0*(b-a)-(d+2.0*c)*e &
    +f*(2.0*(a-b)+(c+d)*e)))

    RETURN
    END FUNCTION cubic


    SUBROUTINE saddle(Z1,A1,Z2,A2,E_bar,m_inertia_eff)
!
! FUNCTIONAL DESCRIPTION:
!
! calculates the condition saddle-point deformation + rotational energy for a
! particular binary decay configuration.
!
! DUMMY ARGUMENTS:
!
! INPUT:
!      Z1 (REAL (KIND=r4) ::) proton number of lighter emitted complex fragment
!      A1 (REAL (KIND=r4) ::) nucleon number of lighter emitted complex fragment
!
!      Z2 (REAL (KIND=r4) ::) proton number of heavier emitted complex fragment
!      A2 (REAL (KIND=r4) ::) nucleon number of heavier emitted complex fragment
!
!  output:
!      E_bar (REAL (KIND=r4) ::) - conditional saddle-point deformation plus rotational
!                     energy
!
! IMPLICIT INPUTS:
!
! sad_array - array containing saddle-point energies from finite range model
!             less the mass of the emitted fragments and their coulomb energy
!
!
    REAL (KIND=r4), INTENT(IN) :: A1,A2,Z1,Z2
    REAL (KIND=r4), INTENT(OUT) :: E_bar,m_inertia_eff

    INTEGER :: Ia1,Ia2,Iz1,Iz2,i
    REAL (KIND=r4) :: Mass1,Mass2
    REAL (KIND=r4) :: r1,r2,A_red,E_coul
    REAL (KIND=r4) :: M_inertia_1,M_inertia_2,M_inertia_orbit
    REAL (KIND=r4) :: M_inertia_saddle,m_inertia_k

!    find  the saddle point energy for the configuration of same mass
!    asymmetry but where the z/a of the fragments is the same
!    as the parent nucleus

    ia1 = A1
    ia2 = A2
    iz1 = Z1
    Iz2 = Z2
    i = MIN(ia1,ia2)

!   calculate the mass of the fragments
    CALL getmass(iz1,ia1,MASS_LD=mass1)

    CALL getmass(iz2,ia2,MASS_LD=mass2)

    ! calculate coulomb energy
    R1 = r0*A1**0.3333333
    R2 = r0*A2**0.3333333
    M_inertia_1 = 0.4*A1*R1**2
    M_inertia_2 = 0.4*A2*R2**2
    A_red = A1*A2/(A1+A2)
    M_inertia_orbit = A_red*(R1+R2+Sep)**2
    M_inertia_saddle = M_inertia_1 + M_inertia_2 + M_inertia_orbit
    E_Coul = 1.44*Z1*Z2/(R1+R2+sep)


    E_bar =  sad_array(i) + Mass1 + Mass2 + E_coul
 
    m_inertia_k = m_inertia_1 + m_inertia_2
    M_Inertia_eff = 1.0/(1.0/M_inertia_k - 1.0/M_inertia_saddle )

    RETURN
    END SUBROUTINE saddle

    SUBROUTINE def_rot (ratio,bss)
    REAL (KIND=r4), INTENT(IN) :: ratio
    REAL (KIND=r4), INTENT(OUT) :: bss
    REAL (KIND=r4) :: ex

    IF (ratio == 1.0) THEN
        bs = 1.0
        bc = 1.0
        br_perp = 1.0
        br_para = 1.0
    ELSE
        IF (ratio > 1.0) THEN !Prolate 
            ex = (1.0-1.0/ratio**2.0)**0.5
            Bs = 0.5*(1.0-ex**2.0)**0.3333333*&
           (1.0+ASIN(ex)/ex/(1.0-ex**2.0)**0.5)
            Bc = 0.5*(1.0-ex**2.0)**0.3333333*LOG((1.0+ex)/(1.0-ex))/ex
        ELSE !oblate
            ex = (1.0/ratio**2.0-1.0)**0.5
            Bs = 0.5*(1.0+ex**2.0)**0.333333333*&
           (1.0+LOG(ex + (1.0+ex**2.0)**0.5)/ex/(1.0+ex**2.0)**0.5)
            Bc = (1.0+EX**2.0)**0.33333*ATAN(ex)/ex
        END IF 
        Br_perp = 2.0/(1.0/ratio**0.66666666+ratio**1.33333)
        Br_para = ratio**0.6666666
    END IF    
    bss = bs
    RETURN
    END SUBROUTINE def_rot 


END MODULE yrast_mod