subroutine dirseg(dirsgs,ndir,nadj,madj,x,y,ntot,rw,eps,ind,nerror)

# Output the endpoints of the segments of boundary of Dirichlet
# tiles.  (Do it economically; each such segment once and only once.)
# Called by master.

implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
dimension dirsgs(8,ndir), rw(4), ind(1)
logical collin, adjace, intfnd, bptab, bptcd, goferit

nerror = -1

# Add in some dummy corner points, outside the actual window.
# Far enough out so that no resulting tile boundaries intersect the
# window.

# Note that these dummy corners are needed by the routine `dirout'
# but will screw things up for `delseg' and `delout'.  Therefore
# this routine (`dirseg') must be called ***before*** dirout, and
# ***after*** delseg and delout.

# Dig out the corners of the rectangular window.
xmin = rw(1)
xmax = rw(2)
ymin = rw(3)
ymax = rw(4)

a = xmax-xmin
b = ymax-ymin
c = sqrt(a*a+b*b)

npd  = ntot-4
nstt = npd+1
i = nstt
x(i) = xmin-c
y(i) = ymin-c
i = i+1
x(i) = xmax+c
y(i) = ymin-c
i = i+1
x(i) = xmax+c
y(i) = ymax+c
i = i+1
x(i) = xmin-c
y(i) = ymax+c

do j = nstt,ntot {
	call addpt(j,nadj,madj,x,y,ntot,eps,nerror)
	if(nerror > 0) return
}

# Put the segments into the array dirsgs.

# For each distinct pair of (genuine) data points, find out if they are
# adjacent.  If so, find the circumcentres of the triangles lying on each
# side of the segment joining them.
kseg = 0
do i1 = 2,npd {
	i = ind(i1)
        do j1 = 1,i1-1 {
		j = ind(j1)
                call adjchk(i,j,adjace,nadj,madj,ntot,nerror)
		if(nerror > 0) return
                if(adjace) {
                        xi = x(i)
                        yi = y(i)
                        xj = x(j)
                        yj = y(j)
                        # Let (xij,yij) be the midpoint of the segment joining
                        # (xi,yi) to (xj,yj).
                        xij = 0.5*(xi+xj)
                        yij = 0.5*(yi+yj)
                        call pred(k,i,j,nadj,madj,ntot,nerror)
			if(nerror > 0) return
                        call circen(i,k,j,a,b,x,y,ntot,eps,collin,nerror)
			if(nerror > 0) return
                        if(collin) {
				nerror = 12
				return
			}

                        # If a circumcentre is outside the rectangular window
                        # of interest, draw a line joining it to the mid-point
                        # of (xi,yi)->(xj,yj).  Find the intersection of this
                        # line with the boundary of the window; for (a,b),
                        # call the point of intersection (ai,bi).  For (c,d),
                        # call it (ci,di).
                        call dldins(a,b,xij,yij,ai,bi,rw,intfnd,bptab)
			if(!intfnd) {
				nerror = 16
				return
			}
                        call succ(l,i,j,nadj,madj,ntot,nerror)
			if(nerror > 0) return
                        call circen(i,j,l,c,d,x,y,ntot,eps,collin,nerror)
			if(nerror > 0) return
                        if(collin) {
				nerror = 12
				return
			}
                        call dldins(c,d,xij,yij,ci,di,rw,intfnd,bptcd)
			if(!intfnd) {
				nerror = 16
				return
			}
			goferit = .false.
			if(bptab & bptcd) {
				xm = 0.5*(ai+ci)
				ym = 0.5*(bi+di)
				if(xmin<xm&xm<xmax&ymin<ym&ym<ymax) {
					goferit = .true.
				}
			}
			if((!bptab)|(!bptcd)) goferit = .true.
			if(goferit) {
				kseg = kseg + 1
				if(kseg > ndir) {
					nerror = 15
					return
				}
				dirsgs(1,kseg) = ai
				dirsgs(2,kseg) = bi
				dirsgs(3,kseg) = ci
				dirsgs(4,kseg) = di
				dirsgs(5,kseg) = i1
				dirsgs(6,kseg) = j1
				if(bptab) dirsgs(7,kseg) = 1.d0
				else dirsgs(7,kseg) = 0.d0
				if(bptcd) dirsgs(8,kseg) = 1.d0
				else dirsgs(8,kseg) = 0.d0
			}
                }
        }
}
ndir = kseg

return
end