c SUBROUTINE putcases: put nc(i) random cases into each polygon i of a map
c 5/23/96 1405

	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

c	last segment location of polygon i, in xys array
	islast = 0

c	number of points kept so for all polygons
	ncases = 0
c	number of random point tries so far for all polygons
	ntrytot = 0

c	loop over polygons
	do 20 i=1,np

c	first and last segment locations of polygon i, in xys array
	isfirst = islast+1
	islast = islast+ns(i)

	if (nc(i).le.0) goto 20

c	calculate area of polygon
	call polyarea(ns(i),xys(1,isfirst),area)

	if (area.le.0.) then
c		write(iotmsg,400) i,area
c400		format('error. polygon ',i5,' has area ',f7.3)	
c		call exit(1)
		goto 20
	endif

c	calculate center of bounding horizontal rectangle
	call polycenter(ns(i),xys(1,isfirst),x0,y0)

c	find minimum-area rectangle around polygon
	call rectangle(ns(i),xys(1,isfirst),x0,y0,angle,x1,x2,
     x	y1,y2,areamin,iotmsg)

c	write(iotmsg,100) i,ns(i),area,x0,y0
c100	format(/'poly,points,area,x0,y0=',2i4, 3f7.3)
c	write(iotmsg,200) angle, x1,x2,y1,y2,areamin
c200	format('angle,x1,x2,y1,y2,areamin=',6f7.3)

c	number of points kept so far for polygon i
	nkeep=0
c	number of tries so far for polygon i
	ntry=0
	maxtry=200000

	do 30 j=1,maxtry
c	plot random point in rectangle
	call ptinrect(angle,x0,y0,x1,x2,y1,y2,x,y,iotmsg)

	call ptinpoly(x,y,xyinf(1),xyinf(2),ns(i),xys(1,isfirst),
     x	iotmsg,inpoly)
c	inpoly = 1/0 if x,y is within/outside polygon at xys(1,isfirst)
c	xyinf(1),xyinf(2) is a point guaranteed to be outside the polygon

	ntry=ntry+1			! no of tries for polygon i
	ntrytot=ntrytot+1		! total number of tries
	nkeep=nkeep+inpoly		! number kept for polygon i
	ncases=ncases+inpoly		! number kept for all polygons

c	coordinates and polygon id's of cases kept
	if (inpoly.gt.0) then
		xyc(1,ncases) = x
		xyc(2,ncases) = y
		ipoly(ncases) = i
	endif

	if (nkeep.ge.nc(i)) goto 31

	if (ntry.ge.maxtry) then
c		write(iotmsg,500) nc(i),i,nkeep,ntry,area,areamin
c500		format(/'error. '
c     x		,i8,' cases requested for polygon ',i5,'.'/
c     x		'found only ',i8,' in ',i8,' tries.'/
c     x		'area of polygon =',f7.3/
c     x		'area of bounding rectangle =',f7.3/)
	endif

30	continue

31	continue

c	write(iotmsg,300) i,areamin,area,ntry,nkeep
c300	format('poly,areamin,area,ntry,nkeep=',i5,2f7.3,2i8)

20	continue

c	write(iotmsg,600) ncases,ntrytot,np
c600	format(/'generated ',i7,' cases in ',i7,' tries in ',i5,' polygons.')

	end

c SUBROUTINE polyarea to calculate area of a polygon

	subroutine polyarea(ns,xys,area)
	integer ns			! segments in polygon
	real xys(1:2,1:ns)		! coords of poly boundary segments
					! no repeated endpoints
 					! can be clockwise or counterclockwise
					! bridges connecting islands OK
					! cuts connecting lakes OK
	real area			! area of polygon

	area = 0.
	do 10 i1=1,ns
	if (i1.eq.ns) then
		i2=1
	else
		i2=i1+1
	endif

	area = area + (xys(1,i1)*xys(2,i2)) - (xys(2,i1)*xys(1,i2))
10	continue
	area = 0.5*abs(area)
	end

c SUBROUTINE polycenter to calculate center of bounding horizontal rectangle

	subroutine polycenter(ns,xys,x0,y0)
	integer ns			! segments in polygon
	real xys(1:2,1:ns)		! coords of poly boundary segments
					! no repeated endpoints
 					! can be clockwise or counterclockwise
					! bridges connecting islands OK
					! cuts connecting lakes OK
	real x0,y0			! center of rectangle

	xmin = 9999.
	xmax = -9999.
	ymin = 9999.
	ymax = -9999.

	do 10 i=1,ns
	if (xys(1,i).gt.xmax) xmax=xys(1,i)
	if (xys(1,i).lt.xmin) xmin=xys(1,i)
	if (xys(2,i).gt.ymax) ymax=xys(2,i)
	if (xys(2,i).lt.ymin) ymin=xys(2,i)
10	continue

	x0 = 0.5 * (xmax+xmin)
	y0 = 0.5 * (ymax+ymin)

	end

c SUBROUTINE rectangle: get minimum-area rectangle around polygon

	subroutine rectangle(ns,xys,x0,y0,angle,x1,x2,y1,y2,areamin,
     x	iotmsg)

	integer ns			! segments in polygon
	real xys(1:2,1:ns)		! coords of poly boundary segments
					! no repeated endpoints
 					! can be clockwise or counterclockwise
					! bridges connecting islands OK
					! cuts connecting lakes OK
	real x0,y0			! center of rectangle
	real a				! orientation of minimum-area bounding
					! rectangle, counterclockwise radians
	real x1,x2,y1,y2		! coords x,y of bounding rectangle
					! in rotated system, relative to x0,y0
	real areamin			! area of smallest bounding rectangle
	integer iotmsg			! file for error messages

				! xys(1)= x0 + x*cos(a) - y*sin(a)
				! xys(2)= y0 + x*sin(a) + y*cos(a)
				! xys(1)-x0 = x*cos(a) - y*sin(a)
				! xys(2)-y0 = x*sin(a) + y*cos(a)
				! x =  (xys(1)-x0)*cos(a) + (xys(2)-y0)*sin(a)
				! y = -(xys(1)-x0)*sin(a) + (xys(2)-y0)*cos(a)
	data pi / 3.14159 /

	areamin=1e10
c	try several values of a = 0,1,...90 degrees
	do 10 j=1,91
	ndeg=1*(j-1)
	a=ndeg*pi/180.

c	find xmax,xmin,ymax,ymin in transformed coords
	xmin = 9999.
	xmax = -9999.
	ymin = 9999.
	ymax = -9999.
	do 20 i=1,ns
	x =   ((xys(1,i)-x0)*cos(a)) + ((xys(2,i)-y0)*sin(a))
	y = - ((xys(1,i)-x0)*sin(a)) + ((xys(2,i)-y0)*cos(a))
	if (x.gt.xmax) xmax=x
	if (x.lt.xmin) xmin=x
	if (y.gt.ymax) ymax=y
	if (y.lt.ymin) ymin=y
c	end loop over points
20	continue

c	for this value of a, calculate area of bounding rectangle
	area = (xmax-xmin)*(ymax-ymin)
c	write(iotmsg,100) j,ndeg, area
c100	format('j,ndeg,area=',i3,i4,f10.5)

c	if smallest area so far, save angle,x1,x2,y1,y2,areamin
	if (area.lt.areamin) then
		jbest=j
		angle=a
		areamin=area
		x1=xmin
		x2=xmax
		y1=ymin
		y2=ymax
	endif

c	end loop over angles
10	continue

c	write(iotmsg,200) jbest,areamin
c200	format(/'jbest,areamin=',i3,f10.5)

	end


c SUBROUTINE ptinrect to plot one random point in rectangle

	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

c	coords in rotated system
	xx = (x2-x1)*rand(0) + x1
	yy = (y2-y1)*rand(0) + y1
c	coords in original system
	x = x0 + xx*cos(angle) - yy*sin(angle)
	y = y0 + xx*sin(angle) + yy*cos(angle)

	end

c SUBROUTINE ptsinpoly for multiple test points in one polygon

	subroutine ptsinpoly(npts,x,y,xinf,yinf,n,xy,iotmsg,inpoly)
	integer npts			! number of points to be tested
	real x(1:npts),y(1:npts)	! coords of points to be tested
     	real xinf,yinf			! any point outside polygon
	integer n 			! points in polygon, no repeated end pt
	real xy(1:2,1:n)		! polygon boundary
 					! can be clockwise or counterclockwise
	integer iotmsg			! file for error messages
     	integer inpoly(1:npts)		! output values (1,0) if (in/out)

	do 20 i=1,npts
	call ptinpoly(x(i),y(i),xinf,yinf,n,xy,iotmsg,inpoly(i))
20	continue

	end

c SUBROUTINE ptinpoly for a single test point

	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)

c	local variables
	integer i,nint
	integer xi,yi			! point of intersection
	real f1,f2			! 0<f1<1 and 0<f2<1 if intersection
c

	nint=0
	do 10 i=1,n
	i1=i
	if (i.eq.n) then
		i2=1
	else
		i2=i1+1
	endif
c	returns 0<f1<1 and 0<f2<1 if (x-xinf,y-yinf) intersects segment i
	call xsect(x,y,xy(1,i1),xy(2,i1),xy(1,i2),xy(2,i2),xinf,yinf
     .	,xi,yi,f1,f2,iotmsg)
	if (0.lt.f1.and.f1.lt.1.and.0.lt.f2.and.f2.lt.1) nint=nint+1

c	write(iotmsg,*) 'n,i,i1,i2,x,y=',n,i,i1,i2,x,y
c     	write(iotmsg,*) 'xy=',xy(1,i1),xy(2,i1),xy(1,i2),xy(2,i2)
c     	write(iotmsg,*) 'xinf,yinf,xi,yi,f1,f2,nint=',
c     x		xinf,yinf,xi,yi,f1,f2,nint

10	continue
c	if nint even, point is outside

	if (mod(nint,2).eq.0) then
		inpoly=0
	else
		inpoly=1
	endif

	return
	end

      subroutine xsect(xa,ya,xb,yb,xc,yc,xd,yd,xe,ye,f1,f2,iotmsg)
c		given a segment a-d which crosses a segment b-c,
c		find the intersection e.
c
c	WARNING: if f1<0 or f1>1 or  f2<0 or f2>1, segments do not intersect
c	and the intersection (xe,ye) is not within the segments
c	
c	example:
c	xa=0,ya=1, xd=0,yd=-1
c	xb=-1,yb=0, xc=1,yc=0
c	xe=0,ye=0, f1=0.5,f2=0.5

c	f1 is the fractional distance (e-a)/(d-a)
c	f2 is the fractional distance (e-b)/(c-b)
c
c	xe=xa+f1*(xd-xa)
c	ye=ya+f1*(yd-ya)
c	xe=xb+f2*(xc-xb)
c	ye=yb+f2*(yc-yb)
c	example:
c	0=0+f1*(0-0)
c	0=1+f1*(-1-1)
c	0=-1+f2*(1+1)
c	0=0+f2*(0-0)

c	xa+f1*(xd-xa)=xb+f2*(xc-xb)
c	ya+f1*(yd-ya)=yb+f2*(yc-yb)
c	example:
c	0+f1*(0-0) = -1+f2*(1+1)
c	1+f1*(-1-1) = 0+f2*(0-0)
c	
c	(xd-xa)*f1 + (xb-xc)*f2 = (xb-xa)
c	(yd-ya)*f1 + (yb-yc)*f2 = (yb-ya)
c	example:
c	(0-0)*f1   + (-1-1)*f2  = (-1-0)
c	(-1-1)*f1  + (0-0)*f2   = (0-1)
c	
c	define:
c	A=(xd-xa)	B=(xb-xc)
c	C=(yd-ya)	D=(yb-yc)
c	det=AD-BC
c	example:
c	A=(0-0)=0	B=(-1-1)=-2
c	C=(-1-1)=-2	D=(0-0)=0
c	det=0*0-(-2)(-2)=-4
c
c	inverse matrix
c	A'=D/det	B'=-B/det
c	C'=-C/det	D'=A/det
c	example:
c	A'=0/(-4)=0	B'=2/(-4)=-0.5
c	C'=2/(-4)=-0.5	D'=0/(-4)=0
c
c	f1 = A'*(xb-xa) + B'*(yb-ya)
c	f2 = C'*(xb-xa) + D'*(yb-ya)
c	example:
c	f1 = 0*(-1-0)    + (-0.5)*(0-1) = 0.5
c	f2 = -0.5*(-1-0) + 0*(0-0)      = 0.5
c	

	real A,B,C,D,det,AI,BI,CI,DI,f1,f2
	real xa,ya,xb,yb,xc,yc,xd,yd,xe,ye
	integer iotmsg

	xe=-999.
	ye=-999.

c	A=(xd-xa)	B=(xb-xc)
c	C=(yd-ya)	D=(yb-yc)
c	det=AD-BC

	A=xd-xa
	B=xb-xc
	C=yd-ya
	D=yb-yc
	det=A*D-B*C

	if (det.ne.0) then

		AI=D/det
		BI=-B/det
		CI=-C/det
		DI=A/det

		f1 = AI*(xb-xa) + BI*(yb-ya)
		f2 = CI*(xb-xa) + DI*(yb-ya)

		xe=xa+f1*(xd-xa)
		ye=ya+f1*(yd-ya)
c	the following should give the same result
		xe=xb+f2*(xc-xb)
		ye=yb+f2*(yc-yb)

	endif

c	write(iotmsg,'(a,4f8.4)') 'xsect: xa..xd=',xa,xb,xc,xd
c	write(iotmsg,'(a,4f8.4)') 'ysect: ya..yd=',ya,yb,yc,yd
c	write(iotmsg,'(a,5e12.4)') 'xsect: a,b,c,d,det=',a,b,c,d,det
c	write(iotmsg,'(a,4e12.4)') 'xsect: ai,bi,ci,di=',ai,bi,ci,di
c	write(iotmsg,'(a,4f8.4)') 'xsect: f1,f2,xe,ye=',f1,f2,xe,ye

	return
	end