FUNCTION inside, x, y, px, py

; Purpose: see if point is inside polygon
; Category: maths
; Input: x, y - [vector of] points
;        px,py - points defining polygon (will be closed automatically)
; Author: "Bård Krane" <bard.krane@fys.uio.no>
; Mods: wmc - make it work with x, y as vectors
; More-info: posted to comp.lang.idl-pvwave on Wed, 01 Apr 1998 12:26:38 +0200
; See-also: http://www.ecse.rpi.edu/Homepages/wrf/geom/pnpoly.html for another method, possibly better
; This is better than my routine "is_inside"
; Note: reduce test from 1e-8 to 1e-4, since we are usually in single precision. In fact, "0.01"
;       would do as well.

    sx = size(px)
    sy = size(py)
    IF (sx[0] EQ 1) THEN NX=sx[1] ELSE message,'px not a vector'
    IF (sy[0] EQ 1) THEN NY=sy[1] ELSE message,'py not a vector'
    IF (NX EQ NY) THEN N = NX ELSE message,'Incompatible dimensions'

    tmp_px =  [px, px[0]]                           ; Close Polygon in x
    tmp_py =  [py, py[0]]                           ; Close Polygon in y
    
    i  = indgen(N)                                  ; indices 0...N-1
    ip = indgen(N) + 1                              ; indices 1...N
   
    nn=n_elements(x) 
    X1 = tmp_px(i)#replicate(1,nn)  - replicate(1,n)#reform([x],nn)
    Y1 = tmp_py(i)#replicate(1,nn)  - replicate(1,n)#reform([y],nn)
    X2 = tmp_px(ip)#replicate(1,nn) - replicate(1,n)#reform([x],nn)
    Y2 = tmp_py(ip)#replicate(1,nn) - replicate(1,n)#reform([y],nn)
 
    dp = X2*X1 + Y1*Y2                               ; Dot-product
    cp = X1*Y2 - Y1*X2                               ; Cross-product
    theta = atan(cp,dp)

    ret=replicate(0,n_elements(x))
    i=where(abs(total(theta,1)) GT 0.01,count)
    if (count gt 0) then ret(i)=1
    if (n_elements(ret) eq 1) then ret=ret[0]
    return,ret

END