revdep/library/silicate/deldir/SavedRatfor/qtest.r

In mdsumner/gibble: Geometry Decomposition

subroutine qtest(h,i,j,k,shdswp,x,y,ntot,eps,nerror)

# Test whether the LOP is satisified; i.e. whether vertex j
# is outside the circumcircle of vertices h, i, and k of the
# quadrilateral.  Vertex h is the vertex being added; i and k are
# the vertices of the quadrilateral which are currently joined;
# j is the vertex which is ``opposite'' the vertex being added.
# If the LOP is not satisfied, then shdswp ("should-swap") is true,
# i.e. h and j should be joined, rather than i and k.  I.e. if j
# is outside the circumcircle of h, i, and k then all is well as-is;
# *don't* swap ik for hj.  If j is inside the circumcircle of h,
# i, and k then change is needed so swap ik for hj.
# Called by swap.

implicit double precision(a-h,o-z)
dimension x(-3:ntot), y(-3:ntot)
integer h
logical shdswp

nerror = -1

# Look for ideal points.
if(i<=0) ii = 1
else ii = 0
if(j<=0) jj = 1
else jj = 0
if(k<=0) kk = 1
else kk = 0
ijk = ii*4+jj*2+kk

# All three corners other than h (the point currently being
# added) are ideal --- so h, i, and k are co-linear; so 
# i and k shouldn't be joined, and h should be joined to j.
# So swap.  (But this can't happen, anyway!!!)
# case 7:
if(ijk==7) {
	shdswp = .true.
	return
}

# If i and j are ideal, find out which of h and k is closer to the
# intersection point of the two diagonals, and choose the diagonal
# emanating from that vertex.  (I.e. if h is closer, swap.)
# Unless swapping yields a non-convex quadrilateral!!!
# case 6:
if(ijk==6) {
	xh = x(h)
	yh = y(h)
	xk = x(k)
	yk = y(k)
	ss = 1 - 2*mod(-j,2)
	test = (xh*yk+xk*yh-xh*yh-xk*yk)*ss
	if(test>0.d0) shdswp = .true.
	else shdswp = .false.
# Check for convexity:
	if(shdswp) call acchk(j,k,h,shdswp,x,y,ntot,eps)
	return
}

# Vertices i and k are ideal --- can't happen, but if it did, we'd
# increase the minimum angle ``from 0 to more than 2*0'' by swapping ...
#
# 24/7/2011 --- I now think that the forgoing comment is misleading,
# although it doesn't matter since it can't happen anyway.  The
# ``2*0'' is wrong.  The ``new minimum angle would be min{alpha,beta}
# where alpha and beta are the angles made by the line joining h
# to j with (any) line with slope = -1.  This will be greater than
# 0 unless the line from h to j has slope = - 1.  In this case h,
# i, j, and k are all co-linear, so i and k should not be joined
# (and h and j should be) so swapping is called for.  If h, i,
# j and j are not co-linear then the quadrilateral is definitely
# convex whence swapping is OK.  So let's say swap.
# case 5:
if(ijk==5) {
	shdswp = .true.
	return
}

# If i is ideal we'd increase the minimum angle ``from 0 to more than
# 2*0'' by swapping, so just check for convexity:
# case 4:
if(ijk==4) {
	call acchk(j,k,h,shdswp,x,y,ntot,eps)
	return
}

# If j and k are ideal, this is like unto case 6.
# case 3:
if(ijk==3) {
	xi = x(i)
	yi = y(i)
	xh = x(h)
	yh = y(h)
	ss = 1 - 2*mod(-j,2)
	test = (xh*yi+xi*yh-xh*yh-xi*yi)*ss
	if(test>0.d0) shdswp = .true.
	else shdswp = .false.
# Check for convexity:
	if(shdswp) call acchk(h,i,j,shdswp,x,y,ntot,eps)
	return
}

# If j is ideal we'd decrease the minimum angle ``from more than 2*0
# to 0'', by swapping; so don't swap.
# case 2:
if(ijk==2) {
	shdswp = .false.
	return
}

# If k is ideal, this is like unto case 4.
# case 1:
if(ijk==1) {
	call acchk(h,i,j,shdswp,x,y,ntot,eps) # This checks
                                              # for convexity.
                                              # (Was i,j,h,...)
	return
}

# If none of the `other' three corners are ideal, do the Lee-Schacter
# test for the LOP.
# case 0:
if(ijk==0) {
	call qtest1(h,i,j,k,x,y,ntot,eps,shdswp,nerror)
	return
}

# default:  # This CAN'T happen.
nerror = 7
return

end