FUNCTION poly_inside,X,Y,P
;+
;r=poly_inside(X,Y,P) determines whether a point P is
;within the boundary of a polygon with vertices described
;by the coordinate vectors X and Y, or on the boundary, or
;outside the boundary.
;
;H. Rhody
;October 5, 2004
;-

IF n_params() NE 3 THEN MESSAGE,'poly_inside requires 3 parameters'
npts=n_elements(X)
IF npts NE n_elements(Y) THEN MESSAGE,'X and Y must be of equal length'
IF npts LT 3 THEN MESSAGE,'Polygon must have at least three vertices'

;Find the side indicator for the first segment
V=[x[1]-x[0],y[1]-y[0]] ;Segment vector
W=P-[x[0],y[0]]         ;Corner to point
U=[-V[1],V[0]]          ;Normal to segment
Z=W-inprod(V,W)/inprod(V,V)*V ;Orthogonal component of W
IF inprod(z,z) EQ 0 THEN RETURN,0 ;Point is on the boundary
dir0=inprod(U,Z)         ;Direction of normal and ortho
FOR k=1,NPTS-1 DO BEGIN
	V=[x[(k+1) mod npts]-x[k],y[(k+1) mod npts]-y[k]]
	U=[-V[1],V[0]]
	W=P-[x[k],y[k]]
	Z=W-inprod(V,W)/inprod(V,V)*V ;Orthogonal component of W
	IF inprod(Z,Z) EQ 0 THEN RETURN,0 ;Point is on the boundary
	dir1=inprod(U,Z)
	IF (dir0 GT 0) NE (dir1 GT 0) THEN RETURN,-1
ENDFOR

RETURN,1 ;Point is inside
END