function interior,xp,yp,x,y,double=double,warning=warning
;
; PURPOSE:
;   returns 1 if (x,y) is within polygon formed by vertices xp,yp
;   returns 0 otherwise
; INPUTS:
;   xp,yp = vectors of the polygon's vertices
;   x = scalar or vector x coordinate(s) of the point(s) to be tested
;   y = scalar or vector y coordinate(s) of the point(s) to be tested
;   /double will cause to be run in double precision
;   /warning will give warning message if a point is close to a edge
; METHOD:
;   Uses sum of angles subtended by successive vertices at (x,y)
;   Angle sum is 0 when (x,y) is outside polygon and sum is 2*!pi inside
;   Equivalent to Cauchy integral \intgr_0^{2\pi} dz/(x-z)
; RESTRICTIONS:
;   Will give unpredictable results for (x,y) on the perimeter.
; VERSION HISTORY:
;   ejs Nov 2006

if keyword_set(double) then begin
 xp=double(xp) & yp=double(yp)
 x=double(x) & y=double(y)
endif

np= n_elements(xp)

; convert x to vector if necessary
  x=[x] 
  y=[y]
nx=n_elements(x)
xxp=fltarr(np,nx)
yyp=xxp

for j=0,nx-1 do xxp(*,j)=xp
for j=0,nx-1 do yyp(*,j)=yp

x1=xxp
for j=0,nx-1 do x1(*,j)=xxp(*,j)-x(j)
if (size(x1))[0] eq 1 then x2=shift(x1,1) else x2=shift(x1,1,0)
y1=yyp
for j=0,nx-1 do y1(*,j)=yyp(*,j)-y(j)

if (size(y1))[0] eq 1 then y2=shift(y1,1) else y2=shift(y1,1,0)

tg=atan(x2*y1-x1*y2,x1*x2+y1*y2)

t=abs(total(tg,1)/(2*!pi)) + 1.e-5 ; sum over polygon vertices
if keyword_set(warning) then begin
  w=where(abs(t-1.d0) lt 1.d-6,nw)
  if nw gt 0 then message,'Warning,'+strcompress(nw)+' pts close to edge',/info
endif

return,fix(t)

end