!-------------------------------------------------------------------
    !subroutine : extrapolate
    !
    ! package   : Li3 Potential
    !
    ! Language  : Fortran 77 - freeform
    !
    ! 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,
    !              see section I.2 of the paper JCP 118 5484 (2003).
    !
    ! 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 through extrapolation
    !
    ! modules   :
    !
    ! common    :
    !
    ! notes     :  Needs: parameters.dat (Table VIII of the paper)
    !
    !-------------------------------------------------------------------
      subroutine extrapolate(s1,s2,s3,v)
      USE FileUnits_Module
      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 Pfunc,Qfunc
      data debug /.false./
      data firstcall /.true./
      save

      if(firstcall) then
          firstcall=.false.
          !
          !   Read the parameters of the extrapolation on the first call
          !
          open(unit=1200,file='./Potentials/li3web/parameters.dat'
     + ,status='old')
          read(1200,*) n
          do i=1,n
              read(1200,*) x(i)
          end do
          close(1200)
      end if
      !
      !     These is the integrity basis, eq. !16) of the paper.
      !
      ts(1) =s1
      ts(2) =s2**2+s3**2
      ts(3) =s3**3-3d0*s3*s2**2
      !
      !     Compute the function eq. (18)
      !
      tmp1  = Pfunc(ts,x,n)
      !
      !     Compute the function eq. (19)/Vo
      !
      tmp2  = Qfunc(ts,x,n)
      !
      !     This is eq. (17)
      !      
      v = x(1)*tmp1*tmp2
      if(debug) WRITE(Msg_Unit,'(3f12.5,e14.5)') s1,s2,s3,v
      return
      end