ptinrect.test1 _ function(angle,x0,y0,x1,x2,y1,y2){
# ptinrect.test1
#square with corners at +/- 1;
angle _ 0
x0 _ 0
y0 _ 0
x1 _ -1
x2 _ 1
y1 _ -1
y2 _ 1
ptinrect(angle,x0,y0,x1,x2,y1,y2)
}

ptinrect _ function(angle,x0,y0,x1,x2,y1,y2){
  z _ .Fortran("ptinrect",
	as.single(angle),
	as.single(x0),
	as.single(y0),
	as.single(x1),
	as.single(x2),
	as.single(y1),
	as.single(y2),
	x = as.single(0),
	y = as.single(0),
	iotmsg = as.integer(0)
	)

list(z$x, z$y)

#	interface to:
#	subroutine ptinrect(angle,x0,y0,x1,x2,y1,y2,x,y,iotmsg)
#	real angle			! angle of rectangle orientation
#					! radians, measured counterclockwise
#	real x0,y0			! center of rectangle
#	real x1,x2,y1,y2		! coords x,y of bounding rectangle
#					! in rotated system, relative to x0,y0
#	real x,y			! coords of random point, in original
#					! coordinate system
#	integer iotmsg			! file for error messages (not used)
}

ptsinpoly.test1 _ function(x,y,xinf,yinf,xy,inpoly){
# ptsinpoly.test1
#square with corners at +/- 1;  (xinf,yinf) lies outside:
#three points, two in and one out:
xinf_100
yinf_100
xy_c(-1,-1, -1,1, 1,1, 1,-1)
x_c(.5,.8,1.2)
y_c(.4,.7,1.3)
inpoly_ptsinpoly(x,y,xinf,yinf,xy)
#length.x_length(x)
#sum.inpoly_sum(inpoly)
}

ptsinpoly.test2 _ function(x,y,xinf,yinf,xy,inpoly,length.x,sum.inpoly){
# ptsinpoly.test2
#square with corners at +/- 1;  (xinf,yinf) lies outside:
#10 random points in a square with corners at +/- 2:
xinf_100
yinf_100
xy_c(-1,-1, -1,1, 1,1, 1,-1)
npts_10
a_rep(1,npts)
x_-2+2*runif(a)
y_-2+2*runif(a)
inpoly_ptsinpoly(x,y,xinf,yinf,xy)
#length.x_length(x)
#sum.inpoly_sum(inpoly)
}

ptsinpoly.test3 _ function(x,y,xinf,yinf,xy,inpoly,length.x,sum.inpoly){
# ptsinpoly.test3
#square with corners at +/- 1;  (xinf,yinf) lies outside:
# 200000 random points in a square with corners at +/- 2:
xinf_100
yinf_100
xy_c(-1,-1, -1,1, 1,1, 1,-1)
npts_200000
a_rep(1,npts)
x_-2+2*runif(a)
y_-2+2*runif(a)
inpoly_ptsinpoly(x,y,xinf,yinf,xy)
#length.x_length(x)
#sum.inpoly_sum(inpoly)
}

ptsinpoly _ function(x,y,xinf,yinf,xy){
# ptsinpoly - put random cases into each polygon of a map
  z _ .Fortran("ptsinpoly",
	   n.check = as.integer(length(x)),
	   as.single(x),
	   as.single(y),
	   as.single(xinf),
	   as.single(yinf),
	   n.vert = as.integer(length(xy)/2),
	   as.single(xy),
	   iotmsg = as.integer(0),
	   inpoly = as.integer(rep(0,length(x)))
	   )
  z$inpoly
}

ptinpoly _ function(x,y,xinf,yinf,n,xy){
  z _ .Fortran("ptinpoly",
	   as.single(x),
	   as.single(y),
	   as.single(xinf),
	   as.single(yinf),
	   n = as.integer(length(xy)/2),
	   as.single(xy),
	   iotmsg = as.integer(0),
	   inpoly = as.integer(0)
	   )
  z$inpoly

#
# interface to:
#	subroutine ptinpoly(x,y,xinf,yinf,n,xy,iotmsg,inpoly)
#c	returns 1 if x,y is within polygon xy(1:2,1:n)
#c	returns 0 if x,y is outside polygon xy(1:2,1:n)
#c	(points are counterclockwise, with no repeated endpoint)
#c	xinf,yinf is a point guaranteed to be outside the polygon
#
#	integer n 			! no of points in polygon
#	real x,y			! test point
#	real xinf,yinf			! a point guaranteed outside
#	real xy(1:2,1:n)	! polygon boundary, no repeated endpoint
#	integer iotmsg 			! logical unit for messages
#	integer inpoly			! output: (1/0) if (inside/outside)
}
################### begin putcases.s ########################
putcases.test1 _ function(ns,xys,xyinf,nc) {
# put cases into the polygons 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) (can be clockwise or counterclockwise)
#	do not include external boundary unless it is the only polygon

#	xyinf [1:2] = xy coords of a point outside the map

#	nc[1:np] = number of cases reqested for each

#output:
#	xyc [1:sum(nc), 1:2] = coordinates of points

# set up polygons
# polygons 1-9: squares between +/- 1.5 (36 points total)
# polygon 10: oblique hexagon: 6 points total
# segment counter
is _ 0

# number of segments in each polygon
ns_c(rep(4,9),6)

# xys is a matrix with two columns
# number of rows is sum (ns)
a _ rep(0,(2*sum(ns)))
xys _ matrix(a,ncol=2,byrow=T)
# polygon row
for (ix in 1:3)
	{
	# polygon column
	for (iy in 1:3)
		{
		is _ is+4
		# segment in polygon
		xys[is-3,1] _ -1.5 + (1.0*(ix-1)) + 0.
		xys[is-3,2] _ -1.5 + (1.0*(iy-1)) + 0.
		xys[is-2,1] _ -1.5 + (1.0*(ix-1)) + 0.
		xys[is-2,2] _ -1.5 + (1.0*(iy-1)) + 1.
		xys[is-1,1] _ -1.5 + (1.0*(ix-1)) + 1.
		xys[is-1,2] _ -1.5 + (1.0*(iy-1)) + 1.
		xys[is-0,1] _ -1.5 + (1.0*(ix-1)) + 1.
		xys[is-0,2] _ -1.5 + (1.0*(iy-1)) + 0.
		}
	}
# tenth polygon is an oblique hexagon
xys[37,1] _ 4.
xys[37,2] _ 4.
xys[38,1] _ 4.
xys[38,2] _ 4.05
xys[39,1] _ 6.
xys[39,2] _ 5.05
xys[40,1] _ 8.
xys[40,2] _ 5.05
xys[41,1] _ 8.
xys[41,2] _ 5.
xys[42,1] _ 6.
xys[42,2] _ 5.

# xyinf = coordinates of a point guaranteed outside map
xyinf _ c(100,90)

# number of cases requested in each polygon = 0,20,...90
nc_3*(0:9)

putcases(ns,xys,xyinf,nc)

# returns a matrix with 2 columns and sum(nc) rows
# each row contains the coordinates of a point
}

putcases _ function(ns,xys,xyinf,nc){

nstot _ dim(xys)[1] # number of rows in xys
ww_rep(0,2*nstot)
iy_2*(1:nstot)
ix_iy-1
ww[ix] _ xys[,1] # x values of polygon boundaries
ww[iy] _ xys[,2] # y values of polygon boundaries

  z _ .Fortran("putcases",
           np = as.integer(length(ns)),
	   as.single(xyinf),
	   nstot = as.integer(dim(xys)[1]),
	   as.integer(ns),
	   as.single(ww),
	   nctot = as.integer(sum(nc)),
	   nc = as.integer(nc),
	   xyc = as.single(rep(0,2*sum(nc))),
	   ipoly = as.integer(rep(0,sum(nc))),
	   iotmsg = as.integer(0)
	   )

iy_2*(1:sum(z$nc))
ix_iy-1
xc_z$xyc[ix]
yc_z$xyc[iy]
cbind(xc,yc)

#
# splus interface to

#	subroutine putcases(np,xyinf,nstot,ns,xys,nctot,nc,xyc,
#     x	ipoly,iotmsg)
#
#
#	integer np			! number of polygons in map
#     	real xyinf(1:2)			! coords of any point outside 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)
#					! etc
# 					! can be clockwise or counterclockwise
#					! no external boundary polygon
#					! bridges connecting islands OK
#					! cuts connecting lakes OK
#	integer nctot			! total cases to put in map
#	integer nc(1:np)		! cases to put in each polygon
#	real xyc(1:2,1:nctot)		! output: coords of cases
#	integer ipoly(1:nctot)		! polygon id of each case
#	integer iotmsg			! file for error messages
}
###################  end  putcases.s ########################