;+
; NAME:
;  oc_search_eph_hd_hbt
;
; PURPOSE: (one line)
;  Search the HD catalog using a list of positions in an ephemeris.
;
; DESCRIPTION:
;  Search the HD catalog using a list of positions in an ephemeris.
;
;  Routine is not very fast.  It calls 'scat' once per ephemeris point,
;  and each call takes ~1 second to run.
;
;  This routine currently returns each star ID a maximum of once, no 
;  matter how many times the ephemeris may pass in front of it.  This is
;  a limited and should be fixed, e.g., by calling LAY's oc_which_uniquev.
;
; CATEGORY:
;  Occultations
;  
; CALLING SEQUENCE:
;  stars = oc_search_eph_hd_hbt( eph, ephindex [, NSTARS=nstars] $
;                                [, MISS=miss] $
;                                [, RADIUS=radius] [, /VERBOSE] )
;
; INPUTS:
;  eph      -- Ephemeris.  Structure, or array of structures.
;
; OPTIONAL INPUT PARAMETERS:
;  None
;  
; OPTIONAL OUTPUT PARAMETERS:
;  ephindex -- Index linking stars to ephemeris points at which they 
;              were found.  Ephindex is of length nstars.  Ephindex[7] = 100
;              indicates that the 7th stars was found at the 100th 
;              ephemeris point searched.
;
; KEYWORD INPUT PARAMETERS:
;  MISS     -- Miss distance between object and star.  Double, radians.
;  /VERBOSE -- Print diagnostic messages to screen.
;
; KEYWORD OUTPUT PARAMETERS:
;  NSTARS   -- Number of unique stars found.
;  RADIUS   -- Search radius = sqrt(miss^2 + arclen^2)
;
; OUTPUTS:
;  None
;  
; COMMON BLOCKS:
;  None
;
; SIDE EFFECTS:
;  None
;
; RESTRICTIONS:
;  None
;
; EXAMPLE:
;  stars = oc_search_eph_hd_hbt(eph, ephindex, nstars, MISS=2*arcmin, $
                                radius=radius, /VERBOSE)
; 
; MODIFICATION HISTORY:
;  Written   9-Dec-2005 by Henry Throop, SwRI.  Based on Leslie Young's
;                    program occlpred/uvis/search_hd.pro
;                    LAY and HBT in parallel modified that code into two
;                    strains of the present code.
;  Modified 15-Dec-2005 by HBT.  Fixed ephindex algorithm.
;  Modified 23-Feb-2006 by HBT.  Improved documentation and formatting.
;  Modified 13-Mar-2006 by HBT.  Added MISS= keyword.
;  Modified 14-Mar-2006 by HBT.  Added RADIUS= keyword.  Pared down by 
;                    removing TYC2 matching and plotting routines, and 
;                    putting these into the _test routine.
;  Modified 14-Apr-2006 by Leslie Young. Made loop index long.
;
;-
function oc_search_eph_hd_hbt, eph, ephindex, nstars, $
  MISS=miss, RADIUS=radius, VERBOSE=verbose

    nstars 	= 0
    stars	= -1
    ephindex 	= -1
    neph 	= n_elements(eph)
    method	= 'subcat'
    deg 	= !dpi/180d         ; radian/deg
    arcmin 	= !dpi/(180d*60d)   ; radian/arcmin
    arcsec 	= !dpi/(180d*3600d) ; radian/arcsec

; Coarse because HD is only tabulated to 1 arcmin

    if not(keyword_set(MISS)) then miss 	= 2d*arcmin     

; Calculate the arclength.  It is not passed in but can be looked up from the
; ephemeris points.

    arclen	= sqrt( (eph[0].ra - eph[1].ra)^2 + $
                        (eph[0].dec - eph[1].dec)^2 )

    radius 	= sqrt(miss^2 + (arclen/2)^2 )

; Generate a circle

    cx 		= radius*cos(findgen(361.)*!pi/180)
    cy 		= radius*sin(findgen(361.)*!pi/180)

; For testing: search only a limited number of ephemeris points

;      neph	= 30		

; Loop over all the ephemeris points

    for i=0L, neph-1 do begin
      if keyword_set(VERBOSE) then print, 'Looking up stars for point ' + $
	    st(i) + '/' + st(neph) + ', RA = ' + st(eph[i].ra)

; Search the HD catalog on this ephemeris point

      stars1	= oc_search_pos_single_hd(eph[i].ra, eph[i].dec, $
	  radius, nstars1, method=method, VERBOSE=verbose)
      if keyword_set(VERBOSE) then print, '  Found ' + st(nstars1) + $
	  ' HD stars'

; If stars were found, save them.

      if nstars1 gt 0 then begin
            
        if nstars eq 0 then begin		; First time through...
          stars = stars1
          nstars = nstars1
          ephindex = replicate(i, nstars1)
        endif else begin			; On next iterations...
          stars = [stars,stars1]
          nstars = nstars + nstars1
          ephindex = [ephindex,replicate(i, nstars1)]
        endelse
      endif
    endfor

    if nstars eq 0 then return, -1
    
    index	= frange(0, neph-1, neph)

; Sort the stars by ID
; In reality, this is the incorrect way to do this, since it misses
; cases where the ephemeris loops back over the same star again
; (e.g., famous HST Saturn occultation).  See LAY's oc_which_uniqev
; for a solution.

    stars_sort		= stars[sort(stars.id)]
    ephindex_sort	= ephindex[sort(stars.id)]

    stars_uniq		= stars_sort[uniq(stars_sort.id)]
    ephindex_uniq	= ephindex[uniq(stars_sort.id)]

    stars		= stars_uniq
    ephindex		= ephindex_uniq

    return, stars

end