subroutine trifnd(j,tau,nedge,nadj,madj,x,y,ntot,eps,nerror)

# Find the triangle of the extant triangulation in which
# lies the point currently being added.
# Called by initad.

implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot), xt(3), yt(3)
integer tau(3)
logical adjace

nerror = -1

# The first point must be added to the triangulation before
# calling trifnd.
if(j==1) {
	nerror = 11
	return
}

# Get the previous triangle:
j1     = j-1
tau(1) = j1
tau(3) = nadj(j1,1)
call pred(tau(2),j1,tau(3),nadj,madj,ntot,nerror)
if(nerror > 0) return
call adjchk(tau(2),tau(3),adjace,nadj,madj,ntot,nerror)
if(nerror>0) return
if(!adjace) {
        tau(3) = tau(2)
	call pred(tau(2),j1,tau(3),nadj,madj,ntot,nerror)
	if(nerror > 0) return
}

# Move to the adjacent triangle in the direction of the new
# point, until the new point lies in this triangle.
1       continue
ntau  = 0 # This number will identify the triangle to be moved to.
nedge = 0 # If the point lies on an edge, this number will identify that edge.
do i = 1,3 {
        ip = i+1
        if(ip==4) ip = 1 # Take addition modulo 3.

# Get the coordinates of the vertices of the current side,
# and of the point j which is being added:
	xt(1) = x(tau(i))
	yt(1) = y(tau(i))
	xt(2) = x(tau(ip))
	yt(2) = y(tau(ip))
	xt(3) = x(j)
	yt(3) = y(j)

# Create indicator telling which of tau(i), tau(ip), and j
# are ideal points.  (The point being added, j, is ***never*** ideal.)
	if(tau(i)<=0) i1 = 1
	else i1 = 0
	if(tau(ip)<=0) j1 = 1
	else j1 = 0
	k1 = 0
	ijk = i1*4+j1*2+k1

# Calculate the ``normalized'' cross product; if this is positive
# then the point being added is to the left (as we move along the
# edge in an anti-clockwise direction).  If the test value is positive
# for all three edges, then the point is inside the triangle.  Note
# that if the test value is very close to zero, we might get negative
# values for it on both sides of an edge, and hence go into an
# infinite loop.
	call cross(xt,yt,ijk,cprd)
	if(cprd >= eps) continue
	else if(cprd > -eps) nedge = ip
	else {
		ntau = ip
		break
	}
}

# We've played ring-around-the-triangle; now figure out the
# next move:
switch(ntau) {
        case 0: # All tests >= 0.; the point is inside; return.
                return

# The point is not inside; work out the vertices of the triangle to which
# to move.  Notation: Number the vertices of the current triangle from 1 to 3,
# anti-clockwise. Then "triangle i+1" is adjacent to the side from vertex i to
# vertex i+1, where i+1 is taken modulo 3 (i.e. "3+1 = 1").
        case 1:
                # move to "triangle 1"
                #tau(1) = tau(1)
                tau(2)  = tau(3)
		call succ(tau(3),tau(1),tau(2),nadj,madj,ntot,nerror)
		if(nerror > 0) return
        case 2:
                # move to "triangle 2"
                #tau(1) = tau(1)
                tau(3)  = tau(2)
		call pred(tau(2),tau(1),tau(3),nadj,madj,ntot,nerror)
		if(nerror > 0) return
        case 3:
                # move to "triangle 3"
                tau(1)  = tau(3)
                #tau(2) = tau(2)
		call succ(tau(3),tau(1),tau(2),nadj,madj,ntot,nerror)
		if(nerror > 0) return
}

# We've moved to a new triangle; check if the point being added lies
# inside this one.
go to 1

end