;+
; NAME:
;  oc_search_pos_single_stararr
;  
; PURPOSE: (one line)
;  Searches around a given position in a stararray
;
; DESCRIPTION:
;  Searches around a given position in a stararray
;
; CATEGORY:
;  Star catalogs
;
; CALLING SEQUENCE:
;  oc_search_pos_single_stararr, RA, DEC, RADIUS, NSTARS, stararr=stararr
;
; INPUTS:
;  RA     -- Right ascension of center of search position.  Radians, J2000.
;  DEC    -- Declination     of center of search position.  Radians, J2000.
;  RADIUS -- Search radius.  Radians. 
;
;  Inputs are single values, not arrays.
;
; INPUT KEYWORD PARAMETERS:
;  stararr=stararr - the star array to be matched.
;
; OPTIONAL INPUT PARAMETERS:
;  sorted  true is star array is increasing in ra
;  verbose 1 for informational outputs. 2 for debugging 
;
; KEYWORD OUTPUT PARAMETERS:
;  None
;
; OUTPUTS:
;  result -- Array of structures [nstars], with each element being a 
;    fully-populated star structure containing RA, Dec, magnitude,
;    etc.
;  NSTARS -- Number of stars found
;  
; COMMON BLOCKS:
;  None
;
; SIDE EFFECTS:
;  None
;
; RESTRICTIONS:
;
; EXAMPLE:
;  print, (oc_search_pos_single_stararr(0.1, 0.1, 0.01, nstars, stararr=ss)).ID
;
; MODIFICATION HISTORY:
;  Written 2006-Mar-19  by Leslie Young, SwRI
;-

function oc_search_pos_single_stararr, ra, dec, radius, nstar, $
    stararr=stararr, verbose=verbose, sorted=sorted

    hpi = !dpi/2.d
    twopi =  2.d*!dpi

    ; get the ra and dec ranges
    case (1) of 
        radius ge hpi: begin
            decrange = [-hpi, hpi]
        end
        dec+radius gt hpi: begin
            decrange = [(dec-r) < (!dpi-(dec+radius)), hpi]
        end
        dec+radius gt hpi: begin
            decrange = [-hpi, (dec+r) > (-!dpi+(dec-radius))]
        end
        else: begin
            decrange = [dec - radius, dec + radius]
        end
    endcase
    cosdec = min(cos(decrange))
    if cosdec gt 0.d then dra = radius/cosdec
    case (1) of 
        cosdec le 0.d: begin
            rarange = [[0.d, twopi]]
        end
        ra + dra ge twopi and ra - dra le 0.: begin
            rarange = [[0.d, twopi]]
        end
        ra + dra gt twopi: begin
            rarange = [ [ra - dra, twopi], [0., ra + dra - twopi] ]
        end
        ra - dra lt 0.d: begin
            rarange = [ [twopi + ra - dra, twopi], [0., ra + dra] ]
        end
        else: begin
            rarange = [[ra - dra, ra + dra] ]
        end
    endcase


    nstar = 0
    stars = -1
    nrarange = n_elements(rarange)/2
    for irarange = 0, nrarange-1 do begin
        ; postpone the sorted bit
        ibox = where(stararr.ra ge rarange[0,irarange] and $
                     stararr.ra le rarange[1,irarange] and $
                     stararr.dec ge decrange[0] and $
                     stararr.dec le decrange[1], nbox)
        if nbox gt 0 then begin
            rastar = stararr[ibox].ra
            decstar = stararr[ibox].dec
            igd = where(rd2ang(ra,dec,stararr[ibox].ra,stararr[ibox].dec) $
                        le radius, ngd)
            if ngd gt 0 then begin
                indx = ibox[igd]
                if nstar eq 0 then begin
                    nstar = ngd
                    stars = stararr[indx]
                endif else begin
                    nstar = nstar + ngd
                    stars = [star, stararr[indx]]
                endelse
            endif
        endif
    endfor
   
    if keyword_set(verbose) then if verbose gt 1 then begin
        print, 'ibox = ', ibox
        stop
    end
    return, stars
    
end