;+
;FUNCTION:
;  w = enclosed(x,y [cx,cy],NCIRCS=NCIRCS,COUNT=COUNT)
;PURPOSE:
; Returns the indices of a set of x,y points that are inside a contour.
;INPUT:  x,y:    data set of points.  (x and y must be the same dimension)
;        cx,cy:  vector of x,y pairs that describe a closed contour. 
;        if cx,cy are not provided then the cursor is used to obtain it.
;OUTPUT:
;    W:  Array of indices of x (& y) that are within the contour cx,cy.
;    NCIRCS: Same dimension as x (& y); integer array giving the number
;        of times each point is encircled.
;    COUNT:  Size of array W
;-

function enclosed,x,y,cx,cy  $
  ,xlog=xlog, ylog=ylog $
  ,count = count,ncircs=inside,limits=lim
n = n_elements(x)

if keyword_set(lim) then   plot,x,y,_extra= lim

if n_params() lt 4 or keyword_set(cx) eq 0 then begin
   getxy,cx,cy,/continue,psym=0
   plots,cx[0],cy[0],/continue,psym=0
   xlog = !x.type
   ylog = !y.type
endif

cx1 = keyword_set(xlog) ? alog(cx) : cx
cy1 = keyword_set(ylog) ? alog(cy) : cy

;time = systime(1)

inside = make_array(/long,dim=dimen(x))

xrange = minmax(cx)
yrange = minmax(cy)
w = where((x le xrange[1]) and (x ge xrange[0]) and (y le yrange[1]) and (y ge yrange[0]),count)
if count eq 0 then return,-1

wx= keyword_set(xlog) ? alog(x[w]) : x[w]
wy= keyword_set(ylog) ? alog(y[w]) : y[w]

cross = 0
nc = n_elements(cx1)
for i0=0l, nc-1 do begin
  i1 = (i0+1) mod nc
  dx = (cx1[i1]-cx1[i0])
  if dx ne 0 then begin
    ym = double(cy1[i1]-cy1[i0])/dx*(wx-cx1[i0]) + cy1[i0]
    cross = cross +  (fix(cx1[i0] lt wx) - fix(cx1[i1] le wx)) * (ym lt wy)
  endif
endfor

bad = where(finite(wx + wy) eq 0,c)
if c ne 0 then cross[bad] = 0

inside[w] = cross


; Old method (slower by factor of 2 or 3)
;nc = n_elements(cx1)
;i = nc-1
;lastphi = atan(cy1[i] - wy,cx1[i] -wx)
;tdphi = 0.d
;for i=0l,nc-1 do begin
;   phi = atan(cy1[i] - wy,cx1[i] -wx)
;   dphi = phi - lastphi
;   lastphi = phi
;   tdphi = tdphi + dphi +2*!dpi*( fix(dphi lt -!dpi) - fix(dphi gt !dpi) )
;endfor
;bad = where(finite(tdphi) eq 0,c)
;if c ne 0 then tdphi[bad] = 0
;inside[w] = round(tdphi/2/!dpi)

;print,systime(1)-time,' Seconds'

return, where(inside ne 0,count)
end


;Very old method:
;inside = make_array(/long,dim=dimen(x))
;for i=0l,n-1 do begin
;   phi = atan(cy-y(i),cx-x(i))
;   phi = phi - shift(phi,1)
;   phi = phi +2*!dpi*( fix(phi lt -!dpi) - fix(phi gt !dpi) )
;   tphi = total(phi)
;   inside(i) = round(tphi/2/!dpi)
;;print,i,tphi,inside(i)
;endfor