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='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(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