c Deane Merrill's point-in-polygon routine 4/11/96
c This file is sample.f
c USAGE:
c     	integer inpoly			! output value (1,0) if (inside/outside)
c	real x,y			! point to be tested
c     	real xinf,yinf			! any point outside polygon
c	integer n 			! points in polygon, no repeated endpt
c	real xy(1:2,1:n)		! polygon boundary
c					! can be clockwise or counterclockwise
c	integer iotmsg			! file for error messages
c	call ptinpoly(x,y,xinf,yinf,n,xy,iotmsg,inpoly)

c EXECUTION: (on parep2)
c	module load lang
c	f77 sample.f
c	a.out

c TIMING:
c The following sample program with n=4, npts=100000 takes 5 secs on a Sparc 10.

c SAMPLE PROGRAM:
	program main
     	integer inpoly			! output value
					! (1,-1) if (inside/outside) polygon
	real x,y			! point to be tested
     	real xinf,yinf			! any point outside polygon
	integer n 			! points in polygon, no repeated endpt
	real xy(1:2,1:200)		! polygon boundary, 200 segments max
					! can be clockwise or counterclockwise
	integer iotmsg			! file for error messages

	data xinf,yinf,n,iotmsg/100.,100.,4,6/
c	polygon is square with corners at +/- 1
	xy(1,1) = -1.
	xy(2,1) = -1.
	xy(1,2) = -1.
	xy(2,2) =  1.
	xy(1,3) =  1.
	xy(2,3) =  1.
	xy(1,4) =  1.
	xy(2,4) = -1.
	
	npts = 100000
	ninside = 0

	do 10 i=1,npts
c	random pt in square with corners at +/- 2
	x=-2. +4.*rand(0)
	y=-2. +4.*rand(0)
	call ptinpoly(x,y,xinf,yinf,n,xy,iotmsg,inpoly)
	ninside = ninside + inpoly
	if (i.le.10) write(iotmsg,100) i,x,y,inpoly
100	format('i, x, y, inpoly=',i5,2f10.3,i3)
10	continue

	write(iotmsg,101) npts, ninside
101	format('npts, ninside=',2i6)

	end


	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
		ptinpoly=0
	else
		ptinpoly=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