c parep2:$MDOCS/delaunay/triangulate.f 5/26/95

	program triangulate_test
c	triangulate a map
c	assume polygon 1 is external boundary
c	do not triangulate clockwise negative-area polygon (external boundary)
c	assume all other polygons described counterclockwise

	integer maxnp			! max polygons in map
	data maxnp /1000/

	integer maxstot			! max segments in all polygons
	data maxstot /10000/

	integer maxntri			! max triangles in map
	data maxntri /5000/

	integer np			! number of polygons in map
	data np /3/

	integer nstot			! total segments in all polygons
	data nstot /30/

	integer ns(1:1000)		! number of segments in each polygon
	real xys(1:2,1:10000)		! coords of all poly boundary segments
					! no repeated endpoints
					! poly 1: 1 to ns(1)
					! poly 2: ns(1)+1 to ns(1)+ns(2)
	integer ntri			! output: number of triangles
	real xytri(1:2,1:3,1:5000)	! output: coords of all triangles
	integer iotmsg			! file for error messages

	ns(1)=13
	ns(2)=7
	ns(3)=10

c	polygon 1 (external boundary, clockwise)
	xys(1,1) = 0.
	xys(2,1) = 0.5
	xys(1,2) = 0.
	xys(2,2) = 1.
	xys(1,3) = -1.
	xys(2,3) = 1.
	xys(1,4) = 0.
	xys(2,4) = 2.
	xys(1,5) = 2.
	xys(2,5) = 2.
	xys(1,6) = 2.
	xys(2,6) = 1.
	xys(1,7) = 3.
	xys(2,7) = 1.
	xys(1,8) = 3.
	xys(2,8) = 0.
	xys(1,9) = 2.
	xys(2,9) = 0.
	xys(1,10) = 2.
	xys(2,10) = -1.
	xys(1,11) = 1.
	xys(2,11) = -1.
	xys(1,12) = 1.
	xys(2,12) = 0.
	xys(1,13) = 0.
	xys(2,13) = 0.
c	polygon 2 (7 points)
	xys(1,14) = 0.
	xys(2,14) = 0.
	xys(1,15) = 1.
	xys(2,15) = 0.
	xys(1,16) = 1.
	xys(2,16) = 1.
	xys(1,17) = 0.
	xys(2,17) = 2.
	xys(1,18) = -1.
	xys(2,18) = 1.
	xys(1,19) = 0.
	xys(2,19) = 1.
	xys(1,20) = 0.
	xys(2,20) = 0.5
c	polygon 3 (10 points)
	xys(1,21) = 0.
	xys(2,21) = 2.
	xys(1,22) = 1.
	xys(2,22) = 1.
	xys(1,23) = 1.
	xys(2,23) = 0.
	xys(1,24) = 1.
	xys(2,24) = -1.
	xys(1,25) = 2.
	xys(2,25) = -1.
	xys(1,26) = 2.
	xys(2,26) = 0.
	xys(1,27) = 3.
	xys(2,27) = 0.
	xys(1,28) = 3.
	xys(2,28) = 1.
	xys(1,29) = 2.
	xys(2,29) = 1.
	xys(1,30) = 2.
	xys(2,30) = 2. 

	call triangulate(np,nstot,ns,xys,ntri,xytri,iotmsg)
	write(6,601) np,nstot,ntri
601	format(3I6)
	write(6,602) (k,((xytri(i,j,k),i=1,2),j=1,3),k=1,ntri)
602	format(i3,2f8.2,4x,2f8.2,4x,2f8.2)
	call exit(0)
	end

	subroutine triangulate(np,nstot,ns,xys,ntri,xytri,iotmsg)
	integer np			! number of polygons in map
	integer nstot			! total segments in all polygons
	integer ns(1:np)		! number of segments in each polygon
	real xys(1:2,1:nstot)		! coords of all poly boundary segments
					! no repeated endpoints
					! poly 1: 1 to ns(1)
					! poly 2: ns(1)+1 to ns(1)+ns(2)
	integer ntri			! output: number of triangles
	real xytri(1:2,1:3,1:ntri)	! output: coords of all triangles
	integer iotmsg			! file for error messages

c	internal variables

	integer maxinpoly		! max size of any one polygon
	data maxinpoly /500/
c	ixy = local pointer to global xys array
	integer ixy(500)

	integer maxntri			! max triangles in map
	data maxntri /5000/
	integer ixy3(3,5000)		! xys pointers of triangle vertices

	integer itri1,itri2,itri3
	integer i,k
	integer idebug
	data idebug / 1 /		! 1 for debug output

c	global pointer to xys array
	k=0
c	global number of triangles
	ntri=0
c	polygon counter
	do 10 i=1,np

c	nsi = points in this polygon
	nsi=ns(i)
	do 15 j=1,nsi
	k=k+1
c	ixy = pointer to global xys array, local within polygon
	ixy(j)=k
15	continue

c	skip external boundary
	if(i.eq.1) goto 10

	call deldivide (nsi,ixy,xys,m,ixy3,idebug)

c	on return, m=no. of triangles
	do 20 j=1,m
	ntri=ntri+1
c	triangle pointers to xys array
	itri1=ixy3(1,j)
	itri2=ixy3(2,j)
	itri3=ixy3(3,j)
	xytri(1,1,ntri)=xys(1,itri1)
	xytri(1,2,ntri)=xys(1,itri2)
	xytri(1,3,ntri)=xys(1,itri3)
	xytri(2,1,ntri)=xys(2,itri1)
	xytri(2,2,ntri)=xys(2,itri2)
	xytri(2,3,ntri)=xys(2,itri3)
20	continue
10	continue

	return
	end
c parep2:$MDOCS/delaunay/delaunay.f 5/26/95
c from seedis:[seedis.map_routines]delaunay.sub

ccc 	program deldivide_test
ccc 
ccc c	split a n-sided polygon into a n-2 Delaunay triangles
ccc c	on input, n = no of polygon vertices
ccc c	on input, ixy(1:n) = positions of polygon vertices, counterclockwise
ccc c	xy(1:2,1:n) = coords of polygon vertices.  xy is not changed.
ccc c	on output, m is number of triangles
ccc c	ixy3(1:3,1:m) = positions of triangle vertices, counterclockwise
ccc c	idebug = 1 to get debug output
ccc 
ccc c	n = number of points in polygon
ccc 	integer n
ccc 	data n /7/
ccc 
ccc c 	ixy = positions of polygon coordinates in the array xy
ccc 	integer ixy(7)
ccc 	data ixy /1,2,3,4,5,6,7/
ccc 
ccc c	xy = array of coordinates
ccc 	real xy(2,7)
ccc 	
ccc c	number of triangles
ccc 	integer m
ccc 	data m /0/
ccc 
ccc c	ixy3 = positions of triangle coordinates in the array xy
ccc 	integer ixy3(3,5)
ccc 
ccc c	idebug = 1 to turn debug on
ccc c	integer idebug
ccc 	data idebug /0/
ccc 
ccc c	coordinates counterclockwise
ccc c	#1 is concave and yes 1
ccc c	#2 is concave and ambiguous 2
ccc c	#3 is concave and no 3
ccc c	#4 is concave and no 2
ccc c	#5 is concave and yes 1
ccc c	#6 is convex (no) 5
ccc c	#7 is colinear 4
ccc 	xy(1,1) = 0.
ccc 	xy(2,1) = 0.
ccc 	xy(1,2) = 1.
ccc 	xy(2,2) = 0.
ccc 	xy(1,3) = 1.
ccc 	xy(2,3) = 1.
ccc 	xy(1,4) = 0.
ccc 	xy(2,4) = 2.
ccc 	xy(1,5) = -1.
ccc 	xy(2,5) = 1.
ccc 	xy(1,6) = 0.
ccc 	xy(2,6) = 1.
ccc 	xy(1,7) = 0.
ccc 	xy(2,7) = 0.5
ccc 
ccc 	write(6,601) n,ixy,m
ccc 601	format(i3,5x,7i3,5x,i3)
ccc 	write(6,602) (j, (ixy3(i,j),i=1,3),j=1,5)
ccc 602	format(i3,5x,3i3)
ccc 
ccc 	call deldivide (n,ixy,xy,m,ixy3,idebug)
ccc 
ccc 	write(6,601) n,ixy,m
ccc 	write(6,602) (j, (ixy3(i,j),i=1,3),j=1,m)
ccc 
ccc 	call exit(0)
ccc 	end

	subroutine deldivide (n,ixy,xy,m,ixy3,idebug)
c	split an n-sided polygon into n-2 Delaunay triangles
c	on input, n = no of polygon vertices
c	on input, ixy(1:n) = positions of polygon vertices, counterclockwise
c	xy(1:2,1:n) = coords of polygon vertices.  xy is not changed.
c	on output, m is number of triangles
c	ixy3(1:3,1:m) = positions of triangle vertices, counterclockwise
c	idebug = 1 to get debug output
	
c	n = no of points in polygon
	integer n
c 	ixy = positions of polygon coordinates in the array xy
	integer ixy(n)
c	xy = coordinates
	real     xy(2,n)
c	m = no of delaunay triangles found (should be n-2)
	integer m
c	ixy3 = positions of triangle coordinates in the array xy
	integer ixy3(3,m)
c	idebug = 1 for debug output
	integer idebug

	m=0
c	call delsplit with n=n,n-1,...3
	do 10 i=n,3,-1
	ii = i
	m=m+1
c	write(6,601) ii,(ixy(j),j=1,7),(ixy3(j,m),j=1,3)
601	format(i5,5x,7i3,5x,3i3)
	call delsplit (ii,ixy,xy,ixy3(1,m),idebug)
c	write(6,601) ii,(ixy(j),j=1,7),(ixy3(j,m),j=1,3)
c	write(6,602)
602	format(/1x)
10	continue

	return
	end

ccc 	program delsplit_test
ccc 
ccc c	split a polygon into a Delaunay triangle
ccc c	and a smaller polygon
ccc c	on return (if n>3), n is reduced by one and ixy is shortened
ccc c	ixy3 gives sequence numbers of triangle
ccc c	ixy and ixy3 are in counterclockwise order
ccc c	xy is not changed
ccc 
ccc c	n = number of points in polygon
ccc 	integer n
ccc 	data n /7/
ccc 
ccc c 	ixy = positions of polygon coordinates in the array xy
ccc 	integer ixy(7)
ccc 	data ixy /1,2,3,4,5,6,7/
ccc 
ccc c	xy = array of coordinates
ccc 	real xy(2,7)
ccc 	
ccc c	ixy3 = positions of triangle coordinates in the array xy
ccc 	integer ixy3(3)
ccc 	data ixy3 /0,0,0/
ccc 
ccc c	idebug = 1 to turn debug on
ccc c	integer idebug
ccc 	data idebug /0/
ccc 
ccc c	coordinates counterclockwise
ccc c	#1 is concave and yes 1
ccc c	#2 is concave and ambiguous 2
ccc c	#3 is concave and no 3
ccc c	#4 is concave and no 2
ccc c	#5 is concave and yes 1
ccc c	#6 is convex (no) 5
ccc c	#7 is colinear 4
ccc 	xy(1,1) = 0.
ccc 	xy(2,1) = 0.
ccc 	xy(1,2) = 1.
ccc 	xy(2,2) = 0.
ccc 	xy(1,3) = 1.
ccc 	xy(2,3) = 1.
ccc 	xy(1,4) = 0.
ccc 	xy(2,4) = 2.
ccc 	xy(1,5) = -1.
ccc 	xy(2,5) = 1.
ccc 	xy(1,6) = 0.
ccc 	xy(2,6) = 1.
ccc 	xy(1,7) = 0.
ccc 	xy(2,7) = 0.5
ccc 
ccc 	write(6,601) n,ixy,ixy3
ccc 	call delsplit (n,ixy,xy,ixy3,idebug)
ccc 	write(6,601) n,ixy,ixy3
ccc 601	format(i3,5x,7i3,5x,3i3)
ccc 
ccc 	call exit(0)
ccc 	end
ccc 
	subroutine delsplit (n,ixy,xy,ixy3,idebug)
c	split a polygon into a Delaunay triangle
c	and a smaller polygon
c	on return (if n>3), n is reduced by one and ixy is shortened
c	ixy3 gives sequence numbers of triangle
c	ixy and ixy3 are in counterclockwise order
c	xy is not changed

c	n = no of points in polygon, before and after
	integer n
c 	ixy = positions of polygon coordinates in the array xy.
	integer ixy(n)
c	ixy3 = position of triangle coordinates in the array xy
	integer ixy3(3)
c	xy = coordinates
	real     xy(2,n)
c	idebug = 1 for debug output

c	internal variables
	integer j1,j2,j3,idel,j

	if (n.eq.3) then
c		move triangle to ixy3 array
		ixy3(1) = ixy(1)
		ixy3(2) = ixy(2)
		ixy3(3) = ixy(3)
		ixy(1) = 0
		ixy(2) = 0
		ixy(3) = 0
		n = 0
	else if (n.gt.3) then
		do 10 j2=1,n
			j1 = j2-1
			if (j1.le.0) j1=j1+n
			j3 = j2+1
			if (j3.gt.n) j3=j3-n
c			get absolute positions in xy array
			ixyj1=ixy(j1)
			ixyj2=ixy(j2)
			ixyj3=ixy(j3)
			call delaunay(n,ixy,xy,ixyj1,ixyj2,ixyj3,
     .			idel,idebug)
			if (idel.eq.1.or.idel.eq.2) then
c				found delaunay triangle

c				print coords of tri eliminated
				ielim1=ixy(j1)
				ielim2=ixy(j2)
				ielim3=ixy(j3)
c				write(6,604)
c     .				xy(1,ielim1),xy(2,ielim1),
c     .				xy(1,ielim2),xy(2,ielim2),
c     .				xy(1,ielim3),xy(2,ielim3)
c604				format(2f7.2)

c				coordinates of triangle
				ixy3(1) = ixy(j1)
				ixy3(2) = ixy(j2)
				ixy3(3) = ixy(j3)

c				new coordinates of polygon
c				ixy(j2) = ixy(j2+1)
c				ixy(j2+1) = ixy(j2+2)
c				...
c				ixy(n-1) = ixy(n)
				do 20 j=j2,n-1
					ixy(j)=ixy(j+1)	
20				continue

c				remove last point of polygon
				ixy(n)=0
				n=n-1

c				removed one triangle, all done
				go to 11
			endif
10		continue

11	continue
	endif

	return
	end

ccc 	program delaunay_test
ccc 
ccc c 	is triangle a delaunay triangle?
ccc c	assume points are in counterclockwise order
ccc 
ccc c	n = number of points in polygon
ccc 	integer n
ccc 
ccc c 	ixy = position of a coordinate in the array xy.
ccc 	integer ixy(7)
ccc 
ccc c	xy = array of coordinates
ccc 	real xy(2,7)
ccc 
ccc c	idel = result
ccc c   	idel=1 if concave and yes,
ccc c   	idel=2 if concave and ambiguous,
ccc c   	idel=3 if concave and no,
ccc c   	idel=4 if colinear,
ccc c   	idel=5 if convex,
ccc 	integer idel
ccc 
ccc c	j1,j2,j3 = 3 consecutive vertices
ccc 	integer j1,j2,j3
ccc 
ccc c	idebug = 1 to turn debug on
ccc c	integer idebug
ccc 
ccc 	data n /7/
ccc 	data ixy /1,2,3,4,5,6,7/
ccc 	data idebug /0/
ccc 
ccc c	coordinates counterclockwise
ccc c	#1 is concave and yes 1
ccc c	#2 is concave and ambiguous 2
ccc c	#3 is concave and no 3
ccc c	#4 is concave and no 2
ccc c	#5 is concave and yes 1
ccc c	#6 is convex (no) 5
ccc c	#7 is colinear 4
ccc 	xy(1,1) = 0.
ccc 	xy(2,1) = 0.
ccc 	xy(1,2) = 1.
ccc 	xy(2,2) = 0.
ccc 	xy(1,3) = 1.
ccc 	xy(2,3) = 1.
ccc 	xy(1,4) = 0.
ccc 	xy(2,4) = 2.
ccc 	xy(1,5) = -1.
ccc 	xy(2,5) = 1.
ccc 	xy(1,6) = 0.
ccc 	xy(2,6) = 1.
ccc 	xy(1,7) = 0.
ccc 	xy(2,7) = 0.5
ccc 
ccc 	do 10 j2=1,7
ccc 
ccc 	j1 = j2-1
ccc 	if (j1.le.0) j1=7
ccc 	j3 = j2+1
ccc 	if (j3.ge.8) j3=1
ccc 
ccc 	call delaunay(n,ixy,xy,j1,j2,j3,idel,idebug)
ccc 	write (6,601) j2,idel
ccc 601	format(2i3)
ccc 10	continue
ccc 
ccc 	call exit(0)
ccc 	end

      SUBROUTINE DELAUNAY(n,ixy,xy,j1,j2,j3,IDEL,idebug)
c given a polygon of n sides (and thus n points), check whether the triangle
c formed by points j1,j2,j3 is concave within that polygon and a Delaunay
c triangle.  Returns:
c   idel=1 if concave and yes,
c   idel=2 if concave and ambiguous,
c   idel=3 if concave and no,
c   idel=4 if colinear,
c   idel=5 if convex,
c Ambiguity arises when a point is exactly on the boundary of the circle drawn
c through j1,j2,j3.

c j1,j2,j3 are positions in xy array of three sequential points

      integer ixy(1)
      real     xy(2,1)
      
      xj1=xy(1,j1)
      xj2=xy(1,j2)
      xj3=xy(1,j3)
      yj1=xy(2,j1)
      yj2=xy(2,j2)
      yj3=xy(2,j3)

      dx23=xj3-xj2
      dx31=xj1-xj3
      dx12=xj2-xj1
      sx23=xj3+xj2
      sx31=xj1+xj3
      sx12=xj2+xj1

      dy23=yj3-yj2
      dy31=yj1-yj3
      dy12=yj2-yj1
      sy23=yj3+yj2
      sy31=yj1+yj3
      sy12=yj2+yj1

      term1=(dx23*sy23)-(dy23*sx23)
      term2=(dx31*sy31)-(dy31*sx31)
      term3=(dx12*sy12)-(dy12*sx12)

      area=-0.25*(term1+term2+term3)

c     convex (triangle not in this polygon)
         if (area.lt.0) idel=5
c     colinear
         if (area.eq.0) idel=4
c     concave (triangle in this polygon)
         if (area.gt.0) idel=3

	if (idebug.eq.1) write(6,6004) j1,j2,j3,xj1,yj1,xj2,yj2,
     x	xj3,yj3,area,idel
6004	format(' j1-3,xy=',3i4,6f7.1,e11.4,i4)

      CALL CIRCUM(xj1,yj1,xj2,yj2,xj3,yj3,X0,Y0,RSQ)

c     continue only for concave angles
         if (idel.gt.3) return

	if (idebug.eq.1) then
		if (rsq.ge.0.) then
			rr=sqrt(rsq)
		else
			rr=-999.
		endif
		write (6,6005) x0,y0,rr
6005		format(' x0,y0,radius=',3e15.4)
	endif

c	infinite radius, points collinear, not Delaunay
	if (rsq.eq.-999) goto 30

c     loop over all points in polygon other than j1,j2,j3
c	for cut doughnut polygons must exclude also the duplicate
c	point at (exactly) the same location
c	kj = location of point in xy array
c	j1,j2,j3 = location of three points in xy array
      do 10 k=1,n
      kj=ixy(k)
	if (xy(1,kj).eq.xy(1,j1).and.xy(2,kj).eq.xy(2,j1)) goto 10
	if (xy(1,kj).eq.xy(1,j2).and.xy(2,kj).eq.xy(2,j2)) goto 10
	if (xy(1,kj).eq.xy(1,j3).and.xy(2,kj).eq.xy(2,j3)) goto 10
      dx=xy(1,kj)-x0
      dy=xy(2,kj)-y0
      r2=(dx*dx)+(dy*dy)

c debug
	if (idebug.eq.1) then
	r=sqrt(r2)
	write (6,6006) k,kj,xy(1,kj),xy(2,kj),r
6006	format(' k,kj,x,y,r=',2i5,2f7.0,e15.4)
	endif

c     if any point kj within circle, j1,j2,j3 not Delaunay
         if (r2.lt.rsq) goto 30
c     if any point kj on circle, ambiguous
         if (r2.eq.rsq) goto 20
   10 continue

c     all points kj outside circle, j1,j2,j3 is Delaunay triangle
      idel=1
      return

c     kj exactly on Delaunay circle
   20 idel=2
      return

c     j1,j2,j3 not Delaunay
   30 return

      end



      SUBROUTINE CIRCUM(x1,y1,x2,y2,x3,y3,X0,Y0,RSQ)
c find the circle which passes through (x1,y1),(x2,y2),(x3,y3)
c output = center (x0,y0) and squared radius (rsq)
c if input points colinear, returns x0 = y0 = rsq = -999

       x0 = -999.0
       y0 = -999.0
      rsq = -999.0

      sx13=0.5*(x1+x3)
      sy13=0.5*(y1+y3)
      dx13=(x3-x1)
      dy13=(y3-y1)

      sx12=0.5*(x1+x2)
      sy12=0.5*(y1+y2)
      dx12=(x2-x1)
      dy12=(y2-y1)

      den = (dx13*dy12) - (dx12*dy13)

      if (den.eq.0) return

      xfac1 =  (sy13*dy13) + (sx13*dx13) 
      xfac2 =  (sy12*dy12) + (sx12*dx12) 
      yfac1 =  (sx13*dx13) + (sy13*dy13)
      yfac2 =  (sx12*dx12) + (sy12*dy12)
      xnum  =  (xfac1*dy12) - (xfac2*dy13)
      ynum  =  (yfac1*dx12) - (yfac2*dx13)
      x0    =  xnum/den
      y0    =  -ynum/den
      dx10  =  x1-x0
      dx20  =  x2-x0
      dx30  =  x3-x0
      dy10  =  y1-y0
      dy20  =  y2-y0
      dy30  =  y3-y0
      rsq1  =  (dx10*dx10)+(dy10*dy10)
      rsq2  =  (dx20*dx20)+(dy20*dy20)
      rsq3  =  (dx30*dx30)+(dy30*dy30)

c in debug mode, check that rsq1 = rsq2 = rsq3
      rsq = rsq1

      return
      end



      SUBROUTINE SORT(ia,n,ib)
c ia and ib are integer vectors of length n.
c on return, integers in ib are original positions of elements of ia.
c on return, ia is sorted in ascending order.
      dimension ia(1)
      dimension ib(1)
      do 10 i=1,n
   10 ib(i)=i
      kk=n
   20 if(kk.le.1) return
      kk=kk/2
      imax=n-kk
   30 nswt=0
      do 40 i=1,imax
      ll=i+kk
      if(ia(i).le.ia(ll))goto 40
      iasave=ia(i)
      ibsave=ib(i)
      ia(i)=ia(ll)
      ib(i)=ib(ll)
      ia(ll)=iasave
      ib(ll)=ibsave
      nswt=nswt+1
   40 continue
      if(nswt.gt.0)goto 30
      goto 20
      end



      SUBROUTINE AREA(ipos,n,xy,A)
c ipos is an integer vector of length n.
c a=4*area is determined from the ipos() entries of the xy(2,.) array.
c xy(2,.) array does NOT contain "duplicate" endpoint.
c
c ipos (ixy) is sorted in correct order around polygon.
c the values of ipos are the values of ncoord that were stored
c in ixy, namely the position of the coordinate in the array xy.
      dimension ipos(1),xy(2,1)
      a=0
      do 10 i=1,n-1
      xp=xy(1,ipos(i+1))+xy(1,ipos(i))
      xm=xy(1,ipos(i+1))-xy(1,ipos(i))
      yp=xy(2,ipos(i+1))+xy(2,ipos(i))
      ym=xy(2,ipos(i+1))-xy(2,ipos(i))
      a=a+(xp*ym)-(yp*xm)
   10 continue
      xp=xy(1,ipos(1))+xy(1,ipos(n))
      xm=xy(1,ipos(1))-xy(1,ipos(n))
      yp=xy(2,ipos(1))+xy(2,ipos(n))
      ym=xy(2,ipos(1))-xy(2,ipos(n))
      a=a+(xp*ym)-(yp*xm)
      return
      end