################### begin triangulate.s ########################

triangulate.test1 _ function() {
#
# Delaunay triangulation of a map

#input:
#	ns [1:np] = # of segments in each polygon (np = # of polygons)

#	xys [1:sum(ns), 1:2] = coordinates of polygon boundaries
#	(no repeated endpts) (must be counterclockwise)
#	external boundary is polygon 1, not to be triangulated

#output:
#	xytri [1:6,1:ntri] = coordinates of triangles in Delaunay
#	triangulation

# set up polygons
# polygon 1: external boundary: 13 points in clockwise order
# polygon 2: 7 points
# polygon 3: 10 points

# number of segments in each polygon
ns _ c(13,7,10)

#
# polygon 1 (external boundary, 13 points, clockwise)
p1 _ c(0,.5, 0,1, -1,1, 0,2, 2,2, 2,1, 3,1, 3,0, 2,0, 2,-1, 1,-1, 1,0, 0,0)
# polygon 2 (7 points, counterclockwise)
p2 _ c(0,0, 1,0, 1,1, 0,2, -1,1, 0,1, 0,.5)
# polygon 3 (10 points, counterclockwise)
p3 _ c(0,2, 1,1, 1,0, 1,-1, 2,-1, 2,0, 3,0, 3,1, 2,1, 2,2)

#
# xys is a matrix with two columns
xys _ matrix(c(p1,p2,p3),ncol=2,byrow=T)

#
# given ns and xys, returns a matrix xytri (6 x ntri)
xytri _ triangulate(ns,xys)

# each row contains the coordinates of a triangle
}

triangulate _ function(ns,xys){
#input:
#	ns [1:np] = # of segments in each polygon (np = # of polygons)

#	xys [1:sum(ns), 1:2] = coordinates of polygon boundaries
#	(no repeated endpts) (must be counterclockwise)
#	external boundary is polygon 1, not to be triangulated

#output:
#	xytri [1:6,1:ntri] = coordinates of triangles in Delaunay
#	triangulation
#

np _ length(ns)
nstot _ sum(ns)
iotmsg _ 0
# expected number of triangles = sum(ns-2) for all polygons except #1
ntri _ sum((ns-2)[-1])
xytri _ rep(0,6*ntri)

z _ .Fortran("triangulate",
	as.integer(np),
	as.integer(nstot),
	as.integer(ns),
	as.single(xys),
	as.integer(ntri),
	xx=as.single(xytri),
	as.single(iotmsg)
	)

# convert z$xx to a matrix 6 x ntri and put it out
matrix(z$xx,ncol=6,byrow=T)

}

###################  end  triangulate.s ########################