!-------------------------------------------------------------------
    !subroutine : extrapolate
    !
    ! package   : Li3 Potential
    !
    ! Language  : Fortran 77
    !
    ! author    : F. Colavecchia (flavioc@lanl.gov)
    !
    ! date      : 05/10/02       version:
    ! revision  :                version:
    !
    ! purpose   :  Compute the three body term of the spin aligned
    !              Li3 potential by  extrapolating ab-initio data
    !
    ! input     :  s1  -> s1 symmetric triangular coordinate in atomic units  
    !              s2  -> s2 symmetric triangular coordinate in atomic units  
    !              s3  -> s3 symmetric triangular coordinate in atomic units  
    !
    ! output    :  pot -> three body term  
    !
    ! modules   :
    !
    ! common    :
    !
    ! notes     :  Needs: parameters.dat
    !
    !-------------------------------------------------------------------
      subroutine extrapolate(s1,s2,s3,v)
      implicit none
      real*8 s1,s2,s3,v
      double precision x(20)
!
!     x = vector of fitting parameters
!      
      logical debug,firstcall
      integer n,i
      real*8 tmp1,tmp2
      real*8 ts(3)
      real*8 Tsymm,Psymm
      data debug /.false./
      data firstcall /.true./
      save

      !write(*,*) 'x =',x

      if(firstcall) then
          firstcall=.false.
          !
          !   Read the parameters of the extrapolation
          !
          open(unit=1200,file='parameters.dat',status='old')
          read(1200,*) n
          do i=1,n
              read(1200,*) x(i)
          end do
          close(1200)
      end if
!
      ts(1) =s1
      ts(2) =s2**2+s3**2
      ts(3) =s3**3-3d0*s3*s2**2
!
!     Compute the function
!
      tmp1  = Tsymm(ts,x,n)
      tmp2  = Psymm(ts,x,n)
      v = x(1)*tmp1*tmp2
      if(debug) write(*,'(3f12.5,e14.5)') s1,s2,s3,v
      return
      end


      real*8 function Tsymm(t,x,n)
      implicit none
      integer n
      real*8 t(3),x(n)
      
      real*8 tmp1
      integer i
      
      tmp1 = 1d0-tanh(t(1)*x(2)/2d0)
      Tsymm = tmp1
      
      return
      end
!
!     Linear function of the parameters, quadric in the coordinates.
!
!
      real*8 function Psymm(t,x,n)
      implicit none
      integer n
      real*8 t(3),x(n)

      logical debug
      integer i,j,k,l
      integer i1,i2,i3,i4      
      real*8 tmp1
      data debug /.false./
!     0th order
      tmp1 = 1d0 
!     1st order
      tmp1 = tmp1 + t(1)*x(3)
!     2nd order
      tmp1 = tmp1 + x(4)*t(1)**2+x(5)*t(2)
!     3rd order
      tmp1 = tmp1 + x(6)*t(1)**3+x(7)*t(2)*t(1)+x(8)*t(3)
!     4th order
      tmp1 = tmp1 + x(9)*t(1)**4+x(10)*t(2)*t(1)**2+x(11)*t(1)*t(3)
     >            + x(12)*t(2)**2
!     5th order
      tmp1 = tmp1 + x(13)*t(1)**5+x(14)*t(2)*t(1)**3+x(15)*t(1)**2*t(3)
     >            + x(16)*t(1)*t(2)**2+x(17)*t(2)*t(3)
      
      Psymm = tmp1
      
      return
      end
!
!     Gradients of the functions
!
      subroutine dTsymm(t,x,grad,n)
      implicit none
      integer n
      real*8 t(3),x(n),grad(n)
      
      real*8 dtmp
      integer i
      
      dtmp = -t(1)/2d0*(1d0/cosh(t(1)*x(2)/2d0))**2d0
      grad(2) = dtmp    
      
      return
      end


      subroutine  dPsymm(t,x,grad,n)
      implicit none
      integer n
      real*8 t(3),x(n),grad(n)
      
      integer i,j,k,l
      integer i1,i2,i3,i4
!
!     x(2-4) = C(1-3)
!     x(5-10) = C11,C21,C22,C31,C32,C33
!     x(11-20) = C111,....
!
!     Unfold x
!
      grad(3) = t(1)
      grad(4) = t(1)**2
      grad(5) = t(2)
      grad(6) = t(1)**3
      grad(7) = t(2)*t(1)
      grad(8) = t(3)
      grad(9) = t(1)**4
      grad(10)= t(2)*t(1)**2
      grad(11)= t(1)*t(3)
      grad(12)= t(2)**2         
      grad(13)= t(1)**5         
      grad(14)= t(2)*t(1)**3         
      grad(15)= t(1)**2*t(3)       
      grad(16)= t(1)*t(2)**2         
      grad(17)= t(2)*t(3)
      
      return
      end