MODULE tl_stuff
USE gem_kind
USE trig
IMPLICIT NONE
PRIVATE

PUBLIC :: get_tl,get_tl_half, prepare_deformed_emission, read_tl
PUBLIC :: clear_look_up, prepare_plm2, prepare_interpolation, get_tll
PUBLIC :: sle,inverse,prepare_tl


INTEGER, ALLOCATABLE, DIMENSION(:), PRIVATE :: tl1
LOGICAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE :: look_up
REAL (KIND=r4), ALLOCATABLE, DIMENSION(:,:), PRIVATE    :: tl_table

INTEGER, PARAMETER, PRIVATE :: n = 3
INTEGER, PRIVATE :: i,I_def,ix
INTEGER, PRIVATE :: kk_max,kk
INTEGER, PARAMETER, PRIVATE :: N_tl = 22

REAL (KIND=r4),PRIVATE :: eek
REAL (KIND=r4),PRIVATE :: ee_tl_up,ee_tl_down
REAL (KIND=r4),PRIVATE :: r_max,eee_tl

REAL (KIND=r4), DIMENSION(4), PRIVATE :: eetl
REAL (KIND=r4), DIMENSION(3,3), PRIVATE :: mat2,m_inv_low,m_inv_high
REAL (KIND=r4), DIMENSION(3), PRIVATE :: coef_low, coef_high

!REAL (KIND=r4),PARAMETER,PRIVATE :: x1 = 0.693 
!REAL (KIND=r4),PARAMETER,PRIVATE :: x3 = 1.58 
!REAL (KIND=r4),PARAMETER,PRIVATE :: x4 = 2.08 

REAL (KIND=r4),PARAMETER,PRIVATE :: x1 = 0.8
REAL (KIND=r4),PARAMETER,PRIVATE :: x2 = 1.00
REAL (KIND=r4),PARAMETER,PRIVATE :: x3 = 1.3 
REAL (KIND=r4),PARAMETER,PRIVATE :: x4 = 1.7 
REAL (KIND=r4),PARAMETER,PRIVATE :: d_beta_k = 1.0
REAL (KIND=r4),PARAMETER,PRIVATE :: d_beta_d = 10.0

REAL (KIND=r4), ALLOCATABLE,DIMENSION(:,:),PRIVATE,TARGET :: lm_weight,ds
REAL (KIND=r4), ALLOCATABLE,DIMENSION(:,:),PRIVATE :: plm2
REAL (KIND=r4), ALLOCATABLE,DIMENSION(:), PRIVATE :: factor
REAL (KIND=r4), ALLOCATABLE,DIMENSION(:,:,:), PRIVATE,TARGET :: ratio_r
REAL (KIND=r4), POINTER, DIMENSION(:),PRIVATE :: ds0,beta0
REAL (KIND=r4), POINTER, DIMENSION(:,:),PRIVATE :: ratio_r0,lm_weight0
REAL (KIND=r4), ALLOCATABLE,DIMENSION(:),PRIVATE :: ee_tl,xi,tl_s
REAL (KIND=r4), ALLOCATABLE,DIMENSION(:,:),PRIVATE,TARGET :: beta


REAL (KIND=r4), ALLOCATABLE,DIMENSION(:,:,:,:),PRIVATE :: tl_def,tl_pol
REAL (KIND=r4), ALLOCATABLE,DIMENSION(:,:,:), PRIVATE :: tl_coef

REAL (KIND=r4), DIMENSION(0:l_max_quantum,3), PRIVATE :: c1,c2,c3,c4

INTEGER, PRIVATE :: n_points
INTEGER, PRIVATE :: i_beta
INTEGER, PRIVATE :: i_beta_max,mode0
INTEGER, DIMENSION(10), PRIVATE :: n_points_b
INTEGER, DIMENSION(0:l_max_quantum),PRIVATE :: n_points_l,Lbeta 

REAL (KIND=r4), DIMENSION(0:L_max_quantum), PRIVATE :: d_beta
REAL (KIND=r4), PRIVATE :: al1,al2


REAL (KIND=r4), PRIVATE :: scale,expand0

LOGICAL, PRIVATE :: deformed, expanded, modified

CONTAINS

SUBROUTINE get_tl(l,ek,tl,tl_up,tl_down)
  
INTEGER, INTENT(IN) :: l
REAL (KIND=r4), INTENT(IN) :: ek
REAL (KIND=r4), DIMENSION(0:l_max_quantum), INTENT(OUT) :: tl
REAL (KIND=r4), INTENT(OUT) :: tl_up,tl_down
REAL (KIND=r4) :: tl_min
!check to see to transmission coefficients are in look_up table

if (l > l_max_quantum) THEN
     WRITE (*,*) "l>l_amx_quantum",l    !rjc
END IF 
 


if (modified) THEN
    ix = INT(ek+1)
    if (ix > tl1(l)) THEN
        tl(0:l) = 1.0
        RETURN
    ELSE if (ix <= 50 .AND. l <= l_max_quantum) THEN
        IF (look_up(ix,l)) THEN
            IF (para%k_sum) THEN
                kk = l*(l+1)/2 + 1                           
                tl(0:l) = tl_table(kk:kk+l,ix)
            ELSE
                tl(0) = tl_table(l+1,ix)
            END IF 
        RETURN
        END IF 
    END IF 
END IF       
eek = ek
IF (Mode0 <= 1) eek = ek + 1.0
eek = eek*scale

IF (modified) THEN

! ******** ratio=x1 *********************
    call get_tll(c1(l,:),eek,eetl(1))

!***** ratio=x2 ***************************
    call get_tll(c2(l,:),eek,eetl(2))


!************ ratio=x3 **********************
    call get_tll(c3(l,:),eek,eetl(3))

!    prepare for quadratic interpolation
!    for ratio <1.3 low
!    for ratio >1.3 high


     coef_low = MATMUL(m_inv_low,eetl(1:3))

!************ ratio=x4 ***********************
    if (r_max > x3) THEN

        call get_tll(c4(l,:),eek,eetl(4))
        coef_high= MATMUL(m_inv_high,eetl(2:4))
     ELSE
        eetl(4) = -7.0
     END IF 

    if (deformed) THEN

 
         if (para%k_sum) THEN
            i_beta = Lbeta(l)
         ELSE
            i_beta = 1
         END IF
    
         ALLOCATE (ee_tl(n_points_b(i_beta)))
         ALLOCATE (tl_s(n_points_b(i_beta)))

         ratio_r0 => ratio_r(1:n_points_b(i_beta),:,i_beta)
!         ds0 => ds(1:n_points_b(i_beta),i_beta)

!        quadratic interpolation

         WHERE (ratio_r0(:,2) > x3)
             ee_tl = MATMUL(ratio_r0,coef_high)
         ELSEWHERE
             ee_tl = MATMUL(ratio_r0,coef_low)
         END WHERE          


         WHERE (ee_tl > 65.0)
             tl_s = 0.0
         ELSEWHERE
             tl_s = 1.0/(1.0+exp(ee_tl))
         END WHERE


!         WRITE (UNIT=17, FMT=*) eetl
!
!         DO i=1,9
!           WRITE (UNIT=17, FMT=*) ratio_r0(i,2),ee_tl(i),tl_s(i)
!         END DO 

         if (para%k_sum) THEN
            kk = l*(l+1)/2 + 1                           
            lm_weight0 => lm_weight(1:n_points_b(i_beta),kk:kk+l)
            tl(0:l) = MATMUL(tl_s,lm_weight0(:,1:l+1))
            tl_min = MINVAL(tl(0:l))
         ELSE
            lm_weight0 => lm_weight(1:n_points_b(i_beta),1:1)
             tl(0:0) = MATMUL(tl_s,lm_weight0(:,1:1))
            tl_min = tl(0)
         END IF 

        DEALLOCATE (ee_tl)
        DEALLOCATE (tl_s)

    ELSE  if (expanded) THEN ! no deformation, just expanded
         IF (expand0 <= x3) THEN
             eee_tl = coef_low(1)*expand0**2 &
             + coef_low(2)*expand0 &
             + coef_low(3)
         ELSE
             eee_tl = coef_high(1)*expand0**2 &
             + coef_high(2)*expand0 &
             + coef_high(3)
         END IF              

        IF (eee_tl  > 65.0) THEN
           tl(0) = 0.0
        ELSE 
            tl(0) = 1.0/(1.0+exp(eee_tl))
        END IF
        tl_min = tl(0)
        if (para%k_sum .AND. l>0 ) tl(1:l) = tl(0)
    END IF                             

!   load transmission coefficients in lookup table
    if (tl_min >= 0.99) THEN
         tl1(l) = ix
    ELSE IF (ix <= 50) THEN
         if (para%k_sum) THEN
           kk = l*(l+1)/2 + 1                           
           tl_table(kk:kk+l,ix) = tl(0:l)
         ELSE
           tl_table(l+1,ix) = tl(0)
         END IF 
         look_up(ix,l) = .TRUE.
     END IF 


ELSE  ! not deformed or expanded
    call get_tll(c2(l,:),eek,eee_tl)

    IF (eee_tl  > 65.0) THEN
        tl(0) = 0.0
    ELSE 
        tl(0) = 1.0/(1.0+exp(eee_tl))
    END IF

    if (para%k_sum .AND. l>0) THEN
        tl(1:l) = tl(0)
        RETURN
    END IF 


    IF ((mode0 == 2 .or. mode0 == 1) .AND. para%polarization) THEN

        call get_tll(c1(l,:),eek,ee_tl_up)
        call get_tll(c3(l,:),eek,ee_tl_down)


        IF (ee_tl_up  > 65) THEN
            tl_up = 0.0
        ELSE 
            tl_up = 1.0/(1.0+exp(ee_tl_up))
        END IF

        IF (ee_tl_down  > 65) THEN
            tl_down = 0.0
        ELSE 
            tl_down = 1.0/(1.0+exp(ee_tl_down))
        END IF
    END IF
END IF
RETURN
END SUBROUTINE get_tl

SUBROUTINE get_tl_half(mode,iz2,l_min,ek_min)
INTEGER, INTENT(IN) :: mode,iz2, l_min
REAL   (KIND=r4), INTENT(OUT) :: ek_min
REAL (KIND=r4) :: c1,c2,c3
!   find barrier i.e. where tl = .5

    IF (TL_COEF(mode,IZ2,7) == 0) THEN
        WRITE (UNIT=*, FMT=*) "mode=",mode,"IZ2=",IZ2
    END IF

    c1 = tl_coef(mode,iz2,7)

    c2 = tl_coef(mode,iz2,1)              &
    +REAL(l_min+1)*tl_coef(mode,iz2,2)   &
    +REAL(l_min+1)**2*tl_coef(mode,iz2,3)

    c3 = tl_coef(mode,iz2,4)               &
    +REAL(l_min+1)*tl_coef(mode,iz2,5)    &
    +REAL(l_min+1)**2*tl_coef(mode,iz2,6)
    ek_min = c2**2-4.0*c1*c3 

    IF (c1 == 0) THEN
        ek_min = -c3/c2
    ELSE IF (ek_min < 0) THEN  ! Tl never gets to .5, This should never
                             ! happen, but it occasionally does for
                             ! light daughter nuclei and very large
                             ! orbital angular momentum. Such large
                             ! angular momentum contribute little
                             ! to the decay width, so I don"t think this
                             ! is a problem, and this IF statement is
                             ! really here to stop the code from
                             ! crashing
        ek_min = (c3/c1)**0.5
    ELSE
        ek_min = (-c2-ek_min**0.5)/(2.0*c1)   ! Tl=.5 at this kinetic
    END IF
    if (mode <= 1) ek_min = ek_min - 1
    if (ek_min <= 1) ek_min = 1.0

RETURN
END SUBROUTINE get_tl_half
    

SUBROUTINE prepare_deformed_emission (ratio,expand)
!get radius at angle beta to symmetry axis

REAL (KIND=r4), INTENT(IN) :: ratio,expand

INTEGER :: j_beta,ell,m

deformed = .FALSE.
expand0 = 0.0
expanded = .FALSE.

IF (expand /= 0.0 .AND. expand /= 1.) THEN
   expand0 = expand
   expanded = .TRUE. 
END IF 

if (ratio == 1.0 .OR. ratio == 0.0) THEN
  r_max = expand
  RETURN
END IF 

deformed = .TRUE.

IF (ratio < 0.0) STOP !not possible

if (para%k_sum) THEN
  i_beta_max = 10
ELSE
  i_beta_max = 1
END IF 
 
DO j_beta=1,i_beta_max
    if (n_points_b(j_beta) == 0) CYCLE    
    ALLOCATE (xi(n_points_b(j_beta)))
    ALLOCATE (factor(n_points_b(j_beta)))
    beta0 => beta(1:n_points_b(j_beta),j_beta)
    ds0 => ds(1:n_points_b(j_beta),j_beta)
    ratio_r0 => ratio_r(1:n_points_b(j_beta),:,j_beta)

    !the angle xi is defined by
    !x = ratio**(2.0/3.0)*COS(xi)    = a*COS(xi)
    !y = 1.0/ratio**(1.0/3.0)*SIN(xi) = c*SIN(xi)

    
    xi = ATAN(ratio*TAN(beta0))

    !radius
    factor = (SIN(xi)**2+&
    (ratio*COS(xi))**2)**0.5/ratio**(1.0/3.0)

    !surface area
    ds0 = 2.0*pi*ratio**(1.0/3.0)*SIN(xi)*COS(xi)**2&
    /COS(beta0)**2*(COS(xi)**2+(ratio*SIN(xi))**2)**0.5 

    if (expanded) factor = factor*expand

    ratio_r0(:,1) = factor**2
    ratio_r0(:,2) = factor
!   ratio_r0(:,3) = 1.0 this is now done in prepare_interpolation

    DEALLOCATE (xi)
    DEALLOCATE (factor)
END DO

r_max = MAXVAL(ratio_r(:,2,1))

if ( .NOT. para%k_sum) THEN
    lm_weight(:,1) = ds(:,1)/SUM(ds)
    RETURN
END IF 

!neeed to calculate tl are a function of projection
!on symmetry axis
!now need to get |d(l,m,m")|**2 coefficients
!to get m projection 

kk = 0
DO  ell=0,l_max_quantum
    n_points = n_points_l(ell)
    DO m=0,ell
        kk = kk + 1
        lm_weight(1:n_points,kk) = plm2(1:n_points,kk)&
        *ds(1:n_points,Lbeta(ell))
        lm_weight(1:n_points,kk) = lm_weight(1:n_points,kk)&
        /SUM(lm_weight(1:n_points,kk))
    END DO
END DO

RETURN
END SUBROUTINE prepare_deformed_emission


SUBROUTINE read_tl ()
LOGICAL :: be
INTEGER :: mode,iii,jj
!----- normal transmission coeff ----------

ALLOCATE (tl_coef(0:n_tl,100,7))

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

if (be) THEN
    OPEN (                 &
    UNIT=8,                &
    ACCESS = "SEQUENTIAL", &
    FORM = "UNFORMATTED",  &
    ACTION = "READ",       &
    FILE = gemini_home//"tl.gdat",       &
    STATUS = "OLD")
    WRITE (*,*) "opened tl.gdat"
    READ (UNIT=8) tl_coef

    CLOSE (UNIT=8)
ELSE
    OPEN (          &
    UNIT=8,         &
    FILE= gemini_home//"tl.tbl", &  
    STATUS = "OLD", &
    ACTION = "READ")
    WRITE (*,*) "opened tl.tbl"

    READ (UNIT=8, FMT=*)
    DO mode=0,n_tl
       READ (UNIT=8, FMT=*) 
       DO jj=1,100
          READ (UNIT=8, FMT="(7(x,F10.3))") tl_coef(mode,jj,:)
       END DO
       READ (UNIT=8, FMT=*)
    END DO 
    close (UNIT=8)

    OPEN (                 &
    UNIT=8,                &
    ACCESS = "SEQUENTIAL", &
    FORM = "UNFORMATTED",  &
    ACTION = "WRITE",       &
    FILE = gemini_home//"tl.gdat",       &
    STATUS = "NEW")
    
    WRITE (UNIT=8) tl_coef

    CLOSE (UNIT=8)
    WRITE (*,*) "created tl.gdat"
END IF 

!-----  transmission coeff for deformed nuclei------

IF (para%ratio /= 1 .OR. para%K_sum .or. para%expand) THEN
    ALLOCATE (tl_def(n_tl,100,7,3))

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

    if (be) THEN
        OPEN (                 &
        UNIT=8,                &
        ACCESS = "SEQUENTIAL", &
        FORM = "UNFORMATTED",  &
        ACTION = "READ",      &
        FILE = gemini_home//"tl_def.gdat",   &
        STATUS = "OLD")
        WRITE (*,*) "opened tl_def.gdat"
        READ (UNIT=8) tl_def

        CLOSE (UNIT=8)
    ELSE
        OPEN (          &
        UNIT=8,         &
        FILE= gemini_home//"tl_def.tbl", &  
        STATUS = "OLD", &
        ACTION = "READ")
        WRITE (*,*) "opened tl_def.tbl"
         READ (UNIT=8, FMT=*)
         DO mode=1,n_tl
             READ (UNIT=8, FMT=*) 
             DO iii = 1,3
             READ (UNIT=8,FMT=*) 
             DO jj=1,100
                READ (UNIT=8, FMT="(7(x,F10.3))") tl_def(mode,jj,:,iii)
             END DO
             READ (UNIT=8, FMT=*)
             END DO
         END DO 
        close (UNIT=8)

        OPEN (                 &
        UNIT=8,                &
        ACCESS = "SEQUENTIAL", &
        FORM = "UNFORMATTED",  &
        ACTION = "WRITE",       &
        FILE = gemini_home//"tl_def.gdat",       &
        STATUS = "NEW")
 
        WRITE (UNIT=8) tl_def

        CLOSE (UNIT=8)
        WRITE (*,*) "created tl_def.gdat"
    END IF 
END IF 
!----- transmission coeff with l.s term for neutons and protons

IF (para%polarization) THEN
    ALLOCATE (tl_pol(2,100,7,2))

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

    if (be) THEN
        OPEN (                 &
        UNIT=8,                &
        ACCESS = "SEQUENTIAL", &
        FORM = "UNFORMATTED",  &
        ACTION = "READ",      &
        FILE = gemini_home//"tl_pol.gdat",   &
        STATUS = "OLD")
        WRITE (*,*) "opened tl_pol.gdat"
        READ (UNIT=8) tl_pol

        CLOSE (UNIT=8)
    ELSE
        OPEN (          &
        UNIT=8,         &
        FILE= gemini_home//"tl_pol.tbl", &  
        STATUS = "OLD", &
        ACTION = "READ")
        WRITE (*,*) "opened tl_pol.tbl"
         READ (UNIT=8, FMT=*)
         DO mode=1,2
             READ (UNIT=8, FMT=*)
             DO iii = 1,2
             READ (UNIT=8,FMT=*) 
             DO jj=1,100
                READ (UNIT=8, FMT=*) tl_pol(mode,jj,:,iii)
             END DO
             READ (UNIT=8, FMT=*)
             END DO
         END DO 
        close (UNIT=8)

        OPEN (                 &
        UNIT=8,                &
        ACCESS = "SEQUENTIAL", &
        FORM = "UNFORMATTED",  &
        ACTION = "WRITE",       &
        FILE = gemini_home//"tl_pol.gdat",       &
        STATUS = "NEW")

        WRITE (UNIT=8) tl_pol

        CLOSE (UNIT=8)
        WRITE (*,*) "created tl_pol.gdat"
    END IF 
END IF

if (para%ratio < 0) r_max = -para%ratio

if (para%k_sum .OR. para%ratio /= 1.0 .or. para%expand) THEN
    call prepare_interpolation ()
    ALLOCATE (look_up(50,0:l_max_quantum))
    ALLOCATE (tl1(0:l_max_quantum))
    ALLOCATE (tl_table(kk_max,50))



    if (para%k_sum) THEN
        ALLOCATE (lm_weight(i_def,kk_max))
        ALLOCATE (plm2(i_def,kk_max))
        CALL prepare_plm2 ()
    ELSE
        ALLOCATE (lm_weight(i_def,1))
    END IF
END IF 
RETURN
END SUBROUTINE read_tl

SUBROUTINE clear_look_up ()
if (.not. modified) RETURN
look_up = .FALSE.
tl1 = 10000
RETURN
END SUBROUTINE clear_look_up


SUBROUTINE prepare_plm2 ()
LOGICAL :: be
INTEGER :: lll,ii_def,ell,m,kk,k

REAL (KIND=r8), DIMENSION(0:l_max_quantum) :: coef, coef_old, coef_old_old
REAL (KIND=r8), DIMENSION(:), ALLOCATABLE :: add


REAL (KIND=r8) :: yy,xx



!see if file containing plm2 values exists
INQUIRE (FILE = gemini_home//"plm2.gdat", &
 EXIST = be)

!open file and read plm2 numbers
if (be) THEN
   OPEN (             &
   UNIT=2,            &
   FILE = gemini_home//"plm2.gdat", &
   STATUS = "OLD",    &
   ACTION = "READ",    &
   FORM = "UNFORMATTED")   
   WRITE (*,*) "opened plm2.gdat"
   READ (UNIT=2) lll,ii_def
        if (lll == l_max_quantum&
        .AND. ii_def == i_def)THEN
        READ (UNIT=2) plm2
        CLOSE (UNIT=2)
        RETURN
   END IF  
   CLOSE (UNIT=2)
END IF 


ALLOCATE (add(i_def))



plm2 = 0.0
kk = 0
DO ell=0,l_max_quantum
   i_beta = Lbeta(ell)
    DO m=0,ell
        kk = kk + 1
        IF (m == 0) THEN
            SELECT CASE (ell)
            CASE (0) 
                coef(0) = 1.0
            CASE (1) 
                coef(0) = 0.0
                coef(1) = 1.0
            CASE DEFAULT
                DO k=0,ell
                    coef(k) = -REAL(ell-1,KIND=r8)*coef_old_old(k)&
                    /REAL(ell,KIND=r8)
                    IF (k > 0)&
                    coef(k) = coef(k) +&
                    REAL(2*ell-1,KIND=r8)*coef_old(k-1)/REAL(ell,KIND=r8)
                END DO
            END SELECT
            coef_old_old = coef_old
            coef_old = coef
        ELSE !derivative of Pl
            DO k=0,ell-m
                coef(k) = coef(k+1)*REAL(k+1,KIND=r8)
            END DO
            coef(ell-m+1) = 0.0
        END IF

    
        loop_beta4: DO i=1,n_points_l(ell)
            xx = COS(beta(i,i_beta))
            yy = 0.0
            DO k=0,ell-m
                yy = yy + coef(k)*xx**REAL(k,KIND=r8)
            END DO
            add(i) = REAL(SIN(beta(i,i_beta)),KIND=r8)&
            **REAL(2*m,KIND=r8)*yy**2
        END DO loop_beta4
        plm2(:,kk) =  REAL(add(:)/SUM(add),KIND=r4)
    END DO
END DO

DEALLOCATE (add)


!having calculated the plm2"s, put them into a file
OPEN (             &
UNIT=2,            &
FILE = gemini_home//"plm2.gdat", &
STATUS = "NEW",    &
ACTION = "WRITE",    &
FORM = "UNFORMATTED")   
WRITE (UNIT=2) l_max_quantum,i_def
WRITE (UNIT=2) plm2
CLOSE (UNIT=2)
WRITE (*,*) "created plm2.gdat"
RETURN
END SUBROUTINE prepare_plm2

SUBROUTINE prepare_interpolation ()
INTEGER :: ell

!quadratic interpolation
mat2(1,1) =  x1**2
mat2(1,2) =  x1
mat2(1,3) =  1.0
mat2(2,1) = x2**2
mat2(2,2) = x2
mat2(2,3) = 1.0
mat2(3,1) = x3**2
mat2(3,2) = x3
mat2(3,3) = 1.0


CALL inverse(mat2,m_inv_low)

mat2(1,1) =  x2**2
mat2(1,2) =  x2
mat2(1,3) =  1.0
mat2(2,1) = x3**2
mat2(2,2) = x3
mat2(2,3) = 1.0
mat2(3,1) = x4**2
mat2(3,2) = x4
mat2(3,3) = 1.0

CALL inverse(mat2,m_inv_high)

!cubic interpolation
!mat3(1,1) =  x1**3.0
!mat3(1,2) =  x1**2
!mat3(1,3) =  x1
!mat3(1,4) = 1.0
!mat3(2,1) = x2**3.0
!mat3(2,2) = x2**2
!mat3(2,3) = x2
!mat3(2,4) = 1.0
!mat3(3,1) = x3**3.0
!mat3(3,2) = x3**2
!mat3(3,3) = x3
!mat3(3,4) = 1.0
!mat3(4,1) = x4**3.0
!mat3(4,2) = x4**2
!mat3(4,3) = x4
!mat3(4,4) = 1.0

!CALL inverse(mat3,m3_inv)

IF (para%k_sum) THEN
    n_points_b = 0
    DO ell=0,l_max_quantum
       IF (ell == 0) THEN
           i_beta = 10
       ELSE
           i_beta = INT(50.0/REAL(ell))
       END IF 
       if (i_beta > 10) i_beta = 10
       if (i_beta <= 1) i_beta = 1
       d_beta(ell) = REAL(i_beta)
       Lbeta(ell) = i_beta
       n_points_l(ell) = INT((90.0-d_beta(ell)/2.0)/d_beta(ell)) + 1
       n_points_b(i_beta) =  n_points_l(ell)
    END DO
    I_def = n_points_l(l_max_quantum)
    ALLOCATE (beta(i_def,10))
    ALLOCATE (ds(i_def,10))
    ALLOCATE (ratio_r(i_def,3,10))

     DO i_beta=1,10
        DO i=1,n_points_b(i_beta)
           beta(i,i_beta) = ((REAL(i)-0.5)*REAL(i_beta)-0.01)*pi/180.0
        END DO 
     END DO

    kk_max = l_max_quantum*(l_max_quantum+1)/2+l_max_quantum+1
ELSE
    IF (para%ratio > 0.0) THEN
        i_def = INT((90.0-d_beta_d/2.0)/d_beta_d) + 1    
        n_points_b(1) = i_def
        ALLOCATE (beta(i_def,1))
        ALLOCATE (ds(i_def,1))
        ALLOCATE (ratio_r(i_def,3,1))

        DO i=1,i_def
           beta(i,1) = ((REAL(i)-0.5)*d_beta_d-0.01)*pi/180.0
        END DO 
    ELSE  ! gaussian distribution of Coulomb barriers
        i_def = 11
        n_points_b(1) = i_def
        ALLOCATE (ds(i_def,1))
        ALLOCATE (ratio_r(i_def,3,1))
    END IF 
    kk_max = l_max_quantum + 1
    
END IF     

ratio_r(:,3,:) = 1.0

RETURN
END SUBROUTINE prepare_interpolation

PURE SUBROUTINE get_tll(c,ek,eetl)
   
![subroutine_header_comments]
REAL (KIND=r4), DIMENSION(3), INTENT(IN) :: c
REAL (KIND=r4), INTENT(IN) :: ek
REAL (KIND=r4), INTENT(OUT) :: eetl
REAL (KIND=r4) :: emax,d2e
!the fits to the transmission coefficients can occosionally
! have problems. check to see if the fit has a maximum
! for positive ek, if so set tl=0 for ek below this value

if (c(1) == 0.0) THEN
    emax = 0.0
    d2e = 0.0
ELSE
    emax = c(3)/c(1)
    IF (emax > 0.0) THEN     !real stationary points
        emax = emax**0.5
        d2e = 2.0*c(3)/emax**3.0
    ELSE
        emax = 0.0              !no real maximum
        d2e = 0.0
    END IF
END IF                          

if (ek < emax .AND. d2e < 0.0) THEN !fit has a maximum at
                                         !positive ek, set tl=0 below
                                         !this energy
  eetl = 100.0
ELSE IF (ek > emax .and. d2e > 0.0) THEN !fit has minimum at
                                          !positive ek, set eetl at
                                          !it minimum value above this
                                          !energy
  eetl = c(1)*emax + c(2) + c(3)/emax
ELSE
  eetl = c(1)*ek + c(2) + c(3)/ek
END IF 

RETURN
END SUBROUTINE get_tll


PURE SUBROUTINE inverse(m,m_inv)
!finds the inverve of a square matrix

REAL (KIND=r4), INTENT(IN), DIMENSION(:,:)  :: m
REAL (KIND=r4), INTENT(OUT), DIMENSION(:,:) :: m_inv

INTEGER :: i,n
REAL (KIND=r4), DIMENSION(SIZE(m,1)) :: y
REAL (KIND=r4), DIMENSION(SIZE(m,1),SIZE(m,2)) :: mm

n = SIZE(m,1)

DO i=1,n
    y = 0.0
    y(i) = 1.0
    mm = m
    CALL sle(mm,y)
    m_inv(:,i) = y
END DO

RETURN
END SUBROUTINE inverse


PURE SUBROUTINE sle(m,y)
         
!solves a system of linear equations by row reduction method
REAL (KIND=r4), INTENT(INOUT), DIMENSION(:,:) :: m
REAL (KIND=r4), INTENT(INOUT), DIMENSION(:)   :: y

INTEGER :: i,j,n
REAL (KIND=r4) :: const
n = SIZE(y)

DO j=1,n
    const = m(j,j)
    m(j,j:n) = m(j,j:n)/const
    y(j) = y(j)/const
    DO i=1,n
        IF (i /= j) THEN
            const = m(i,j)
            m(i,j:n) = m(i,j:n) - const*m(j,j:n)
            y(i) = y(i) - const*y(j)
        END IF
    END DO
END DO
RETURN
END SUBROUTINE sle

SUBROUTINE prepare_tl(mode,iz2,ia2)
!transmission were fit to calculations for residues
!which lie on the evaporation attracter line (a_eal) for each Z value.
!However there is a small A dependence on the Coulomb barriers and
!transmission coefficients
!and this subroutine is a simple scaling of kinetic energy which
!I a find accounts for this reasonably well

INTEGER, INTENT(IN) :: mode,iz2,ia2
REAL (KIND=r4) :: a_eal,delta_a
INTEGER :: l


mode0 = mode
scale = 1.0
!correction for A (daughter) dependence of coulomb barriers
IF (iz2 > 20) THEN
SELECT CASE (mode)
 CASE(2)
     a_eal = 2.045*REAL(iz2) + 3.57e-3*REAL(iz2)**2
     delta_a = REAL(ia2) - a_eal
     scale = 1.+delta_a*(-5.412976E-04-9.176213E-02/REAL(iz2)&
     + 3.916151E-06*REAL(iz2))
     scale = 1./scale
  CASE(7)
     a_eal = 2.045*REAL(iz2) + 3.57e-3*REAL(iz2)**2
     delta_a = REAL(ia2) - a_eal
     scale = 1.+delta_a*(-8.185244E-04-8.234510E-02/REAL(iz2)&
     + 1.018610E-05*REAL(iz2))
     scale = 1./scale
  END SELECT
END IF 

modified = deformed .OR. expanded
if (mode < 1) modified = .FALSE.
if (iz2 < 20) modified = .FALSE.


   DO l=0,l_max_quantum
    al1 = REAL(l+1)
    al2 = al1**2
   
    c2(l,1) = tl_coef(mode,iz2,7)

    c2(l,2) = tl_coef(mode,iz2,1) &
    +tl_coef(mode,iz2,2)*al1      &
    +tl_coef(mode,iz2,3)*al2

    c2(l,3) = tl_coef(mode,iz2,4) &
    +tl_coef(mode,iz2,5)*al1      &
    +tl_coef(mode,iz2,6)*al2
   END DO 

    IF (modified) THEN
       DO l=0,l_max_quantum
        al1 = REAL(l+1)
        al2 = al1**2
   
        c1(l,1) = tl_def(mode,iz2,7,1)

        c1(l,2) = tl_def(mode,iz2,1,1) &
        +tl_def(mode,iz2,2,1)*al1     &
        +tl_def(mode,iz2,3,1)*al2

        c1(l,3) = tl_def(mode,iz2,4,1)&
        +tl_def(mode,iz2,5,1)*al1    &
        +tl_def(mode,iz2,6,1)*al2

        c3(l,1) = tl_def(mode,iz2,7,2)

        c3(l,2) = tl_def(mode,iz2,1,2) &
        +tl_def(mode,iz2,2,2)*al1     &
        +tl_def(mode,iz2,3,2)*al2

        c3(l,3) = tl_def(mode,iz2,4,2)&
        +tl_def(mode,iz2,5,2)*al1     &
        +tl_def(mode,iz2,6,2)*al2

        c4(l,1) = tl_def(mode,iz2,7,3)

        c4(l,2) = tl_def(mode,iz2,1,3) &
        +tl_def(mode,iz2,2,3)*al1     &
        +tl_def(mode,iz2,3,3)*al2

        c4(l,3) = tl_def(mode,iz2,4,3)&
        +tl_def(mode,iz2,5,3)*al1     &
        +tl_def(mode,iz2,6,3)*al2
     END DO 
   END IF 

   IF ((mode == 1 .OR. mode==2) .AND. para%polarization) THEN
       DO l=0,l_max_quantum
        al1 = REAL(l+1)
        al2 = al1**2
   
        c1(l,1) = tl_pol(mode,iz2,7,1)

        c1(l,2) = tl_pol(mode,iz2,1,1) &
        +tl_pol(mode,iz2,2,1)*al1     &
        +tl_pol(mode,iz2,3,1)*al2

        c1(l,3) = tl_pol(mode,iz2,4,1)&
        +tl_pol(mode,iz2,5,1)*al1    &
        +tl_pol(mode,iz2,6,1)*al2

        c3(l,1) = tl_pol(mode,iz2,7,2)

        c3(l,2) = tl_pol(mode,iz2,1,2) &
        +tl_pol(mode,iz2,2,2)*al1     &
        +tl_pol(mode,iz2,3,2)*al2

        c3(l,3) = tl_pol(mode,iz2,4,2)&
        +tl_pol(mode,iz2,5,2)*al1     &
        +tl_pol(mode,iz2,6,2)*al2
   END DO
END IF
RETURN
END SUBROUTINE prepare_tl
END MODULE tl_stuff