!-------------------------------------------------------------------
      ! function  : search
      !
      ! package   : Li3 Potential
      !
      ! Language  : Fortran 77
      !
      ! author    : M. Salazar     (msalazar@uu.edu)
      !
      ! date      : 05/10/02       version: 0.1
      ! revision  :                version:
      !
      ! purpose   :  Check if the point {s1,s2,s3} is inside one of
      !              the tetrahydra formed by ab initio points.
      !
      ! input     :  s1 -> s1 triangular coordinate
      !              s2 -> s2 triangular coordinate
      !              s3 -> s3 triangular coordinate
      ! output    :  search -> true if the point lies in a tetrahedra.
      !
      !
      ! notes     : Needs: li3_tetra.dat
      !
      !-------------------------------------------------------------------
      logical function search(s1,s2,s3)
      implicit none
      real*8 s1,s2,s3
      integer ntmax
      parameter (ntmax = 20000)
      integer  i,numtetra,j,q(4),ql(4),co,l,k
      double precision x(4),y(4),z(4),e(4),s(3),b(4)
      integer qv(4,ntmax),qlv(4,ntmax)
      double precision xv(4,ntmax),yv(4,ntmax),zv(4,ntmax),ev(4,ntmax)
      logical firstcall,debug
      data firstcall /.true./
      data debug /.true./

      s(1) = s1
      s(2) = s2
      s(3) = s3

      numtetra=19819
      if(firstcall)  then
          firstcall = .false.
          open(unit=50,access='direct',file='./Potentials/li3web/
     +li3_tets.dat',
     +         status='old',form='formatted',recl=400)
          do i=1,numtetra
              read(50,350,rec=i) j,qv(1,i),qv(2,i),qv(3,i),qv(4,i),
     +            xv(1,i),yv(1,i),zv(1,i),ev(1,i),
     +            xv(2,i),yv(2,i),zv(2,i),ev(2,i),
     +            xv(3,i),yv(3,i),zv(3,i),ev(3,i),
     +            xv(4,i),yv(4,i),zv(4,i),ev(4,i),
     +            qlv(1,i),qlv(2,i),qlv(3,i),qlv(4,i)
          end do
          close(50)
      end if



      do i=1,numtetra
          do k=1,4
             x(k) = xv(k,i)
             y(k) = yv(k,i)
             z(k) = zv(k,i)
             e(k) = ev(k,i)
             q(k) = qv(k,i)
             ql(k) = qlv(k,i)
          end do
          !   compute the barycentric coordinates for s1,s2,s3
          call barycentric(x,y,z,s,b)

          co=0
          do l=1,4
              if (b(l).gt.-1.0d-5) then
                  co=co+1
              end if
          end do

          if (co.eq.4) then
              if(debug) write(6,200) s(1),s(2),s(3),i
              search = .true. 
              return
          end if
      end do
      if (debug) write(6,300) s(1),s(2),s(3)
      search =.false. 
  200	format('the point ',3f10.6,' is in tetrahedra ',i8)
  300	format('the point ',3f10.6,' is outside the range')
  350 format(i6,i6,i6,i6,i6,d12.6,d12.6,d12.6,d12.6,
     +     d12.6,d12.6,d12.6,d12.6,d12.6,d12.6,d12.6,d12.6,
     +     d12.6,d12.6,d12.6,d12.6,i6,i6,i6,i6)
      
      return
      end


      !-------------------------------------------------------------------
      ! subroutine: barycentric
      !
      ! package   : Li3 Potential
      !
      ! Language  : Fortran 77
      !
      ! author    : M. Salazar     (msalazar@uu.edu)
      !
      ! date      : 05/10/02       version: 0.1
      ! revision  :                version:
      !
      ! purpose   :  Computes the barycentric coordinates for the point 
      !              {s(3)} 
      !
      ! input     :  x -> \
      !              y -> - tetrahedra data
      !              z -> /
      !              s -> point in triangular coordinates
      !             
      ! output    :  b -> barycentric coordinates
      !
      !
      !-------------------------------------------------------------------
      subroutine barycentric(x,y,z,s,b)

      double precision x(4),y(4),z(4),s(3),b(4)
      double precision det,a11,a12,a13,a14
      double precision a21,a22,a23,a24,a31,a32,a33,a34
      double precision a41,a42,a43,a44,al1,al2,al3,al4
      double precision bt1,bt2,bt3,bt4,gm1,gm2,gm3,gm4
      double precision dl1,dl2,dl3,dl4
*
* this portion of the program computes the matrix elements.
*
      a11=x(2)*(y(3)*z(4)-y(4)*z(3))-x(3)*(y(2)*z(4)-
     +	y(4)*z(2))+x(4)*(y(2)*z(3)-y(3)*z(2))
      a12=-x(1)*(y(3)*z(4)-y(4)*z(3))+x(3)*(y(1)*z(4)-
     +	y(4)*z(1))-x(4)*(y(1)*z(3)-y(3)*z(1))
      a13=x(1)*(y(2)*z(4)-y(4)*z(2))-x(2)*(y(1)*z(4)-
     +	y(4)*z(1))+x(4)*(y(1)*z(2)-y(2)*z(1))
      a14=-x(1)*(y(2)*z(3)-y(3)*z(2))+x(2)*(y(1)*z(3)-
     +	y(3)*z(1))-x(3)*(y(1)*z(2)-y(2)*z(1))

      a21=y(3)*(z(2)-z(4))+y(2)*(z(4)-z(3))+y(4)*(z(3)-z(2))
      a22=y(3)*(z(4)-z(1))+y(4)*(z(1)-z(3))+y(1)*(z(3)-z(4))
      a23=y(2)*(z(1)-z(4))+y(4)*(z(2)-z(1))+y(1)*(z(4)-z(2))
      a24=y(2)*(z(3)-z(1))+y(3)*(z(1)-z(2))+y(1)*(z(2)-z(3))

      a31=x(3)*(z(4)-z(2))+x(4)*(z(2)-z(3))+x(2)*(z(3)-z(4))
      a32=x(3)*(z(1)-z(4))+x(4)*(z(3)-z(1))+x(1)*(z(4)-z(3))
      a33=x(2)*(z(4)-z(1))+x(4)*(z(1)-z(2))+x(1)*(z(2)-z(4))
      a34=x(3)*(z(2)-z(1))+x(2)*(z(1)-z(3))+x(1)*(z(3)-z(2))

      a41=x(3)*(y(2)-y(4))+x(4)*(y(3)-y(2))+x(2)*(y(4)-y(3))
      a42=x(3)*(y(4)-y(1))+x(4)*(y(1)-y(3))+x(1)*(y(3)-y(4))
      a43=x(2)*(y(1)-y(4))+x(4)*(y(2)-y(1))+x(1)*(y(4)-y(2))
      a44=x(2)*(y(3)-y(1))+x(3)*(y(1)-y(2))+x(1)*(y(2)-y(3))
*
* det is the determinant of the matrix and alpha, beta,
* delta and gamma are the coefficiants of the barycentric
* coordinates.
*
      det=a11+a12+a13+a14

      al1=a11/det
      bt1=a21/det
      gm1=a31/det
      dl1=a41/det

      al2=a12/det
      bt2=a22/det
      gm2=a32/det
      dl2=a42/det

      al3=a13/det
      bt3=a23/det
      gm3=a33/det
      dl3=a43/det

      al4=a14/det
      bt4=a24/det
      gm4=a34/det
      dl4=a44/det
*
* b1,b2,b3, and b4 are the barycentric coordinates in three
* dimensions.
*
      b(1)=al1+bt1*s(1)+gm1*s(2)+dl1*s(3)
      b(2)=al2+bt2*s(1)+gm2*s(2)+dl2*s(3)
      b(3)=al3+bt3*s(1)+gm3*s(2)+dl3*s(3)
      b(4)=al4+bt4*s(1)+gm4*s(2)+dl4*s(3)

      return

      end