;+
; NAME:
;  oc_search_pos_subcat
;
; PURPOSE: (one line)
;  Searches for stars by position, using pre-computed subcatalogs.
;
; DESCRIPTION:
;  This routine searches around a specified (ra, dec, radius) in a 
;  specified star catalog.  It returns a fully-populated star structure
;  for each star found.
;
;  Rather than searching the full catalog, this routine searches subcatalogs
;  which break the catalog into small regions of RA/DEC.  For text catalogs
;  such as HD, this routine is faster than searching the entire catalog.
;
;  If the subcatalogs do not exist, they are created and saved on-the-fly.
;  The first time this routine is called, there is substantial overhead
;  in extracting the subcatalogs.  Subsequent searches are much faster.
;
;  This routine currently works only with the HD catalog.  It could
;  be set up to work with Tycho-2.
;
; CATEGORY:
;  Star catalogs.
;  
; CALLING SEQUENCE: 
;  result = oc_search_pos_subcat(RA, DEC, RADIUS, NSTARS,  $
;                            CATALOG=catalog [, /VERBOSE] )
;
; 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.
;
; OPTIONAL INPUT PARAMETERS:
;  None
;
; KEYWORD INPUT PARAMETERS:
;  CATALOG= -- Name of the catalog to use; e.g., 'HD'.
;    Defaults to 'HD' if not set.
;  /VERBOSE -- If set, print diagnostics to screen.
;
; 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:
;  The proper environment variable must be set for the catalog called.
;   $WCS_CATDIR -- For HD and other WCS-based catalogs
;  The directory $HDOCCUL_TMP must exist and be writeable for the subcatalogs
;   to be stored in.  [2006-Mar-05 - obviated.  LAY]
;  The temporary directory used by tmpfile (nominally /tmp) should
;   exist and be writable for the subcatalogs to be stored in
;   [2006-Mar-5 LAY]
;
; EXAMPLE:
;  print,(oc_search_pos_subcat(0.1, 0.1, 0.01, nstars, catalog='hd')).id
;    [Locates three stars]
; 
; MODIFICATION HISTORY:
;  Written  31-Jan-2005 by Henry Throop, SwRI.
;  Modified 22-Feb-2005 by HBT. Renamed, and generalized to work with other
;                                catalogs.
;  Modified  9-Dec-2005 by HBT. Renamed again.
;  Modified 27-Feb-2006 by HBT. Improved documentation and formatting.
;  Modified 5-Mar-2006 by LAY.  Replaced get_pid with tmpfile
;         Replace ip_down with built-in floor
;         Uncomment the spawn, command_rm
;         Calculate distances of grid centers with rd2ang
;  Modified 14-Apr-2006 by Leslie Young. Made loop index long.
;-

function oc_search_pos_subcat, ra, dec, radius, nstars, $
  VERBOSE=verbose, CATALOG_IN=catalog_in

  d2r		= 2*!dpi / 360d
  r2d		= 1/d2r
  d2as		= 60d*60d

; Maximum number of stars to return

  nstars_max	= 1000000				
  catalog_default= 'hd'

  if not(keyword_set(CATALOG_IN)) then begin
    catalog	= catalog_default
  end $
  else begin
    catalog	= strlowcase(CATALOG_IN)
  end

  if (catalog eq 'ty2') then begin
    print, 'Error: oc_search_pos_fast.pro does not support TY2 catalog'
    stop
  end

;   path_wcs	= getenv('WCS_CATDIR')
  path_subcat	= getenv('HDOCCUL_TMP') + '/'

; grid size, degrees, RA
  dra_deg	= 1d					
; grid size, degrees, Dec
  ddec_deg	= 1d					

; Process the coordinates

  ra_rad	= ra
  dec_rad	= dec
  radius_rad	= radius

  ra_deg	= ra*r2d
  dec_deg	= dec*r2d

  radius_deg	= radius*r2d
  radius_arcsec	= radius_rad * r2d * d2as


; Make a grid of all possible ra, dec grid points

  ra_grid_deg	= frange(0d +dra_deg/2d, 360d -dra_deg/2d, 360d/dra_deg)
  dec_grid_deg	= frange(-90d +ddec_deg/2d, 90d -ddec_deg/2d, 180d/ddec_deg)

; Calc distance from the *center* of each gridpoint, to the given ra+dec

; Original distance calculation by HBT.  This is intended to 
; calculate an n_ra x n_dec
; array, where dist_2d[i_ra, i_dec] = distance from grid center to 
; star.  This is replaced with rd2ang because (1) it fails when 
; ra_grid_deg-ra_deg > 180 (eg when ra_deg = 0 and ra_grid_deg = 360-eps).
; and (2) it fails to account for the cos(dec) factor.
;  dist_2d	= sqrt(coadd( (ra_grid_deg-ra_deg)^2, $
;                              (dec_grid_deg-dec_deg)^2))
  ra_grid = ra_grid_deg*d2r
  dec_grid = dec_grid_deg*d2r
  nra_grid = n_elements(ra_grid_deg)
  ndec_grid = n_elements(dec_grid_deg)
  dist_2d = dblarr(nra_grid, ndec_grid)
  for ira_grid = 0L, nra_grid-1 do begin
      for idec_grid = 0L, ndec_grid-1 do begin
          dist_2d[ira_grid,idec_grid] = $
            rd2ang(ra, dec, ra_grid[ira_grid], dec_grid[idec_grid])*r2d
      end
  end

  ra_deg_2d		= coadd(ra_grid_deg, 0d*dec_grid_deg)
  dec_deg_2d	= coadd(0d*ra_grid_deg, dec_grid_deg)

; To make sure we cover everything, be sure to search to the edge of the box

  in_range	= (dist_2d le (radius_deg + (dra_deg>ddec_deg)*sqrt(2)/2d))

; Figure out which boxes this point covers

  ra_boxes_deg	= ra_deg_2d[where(in_range)]
  dec_boxes_deg	= dec_deg_2d[where(in_range)]

; Get the minimum and maximum values for each gridpoint

  ra_min_deg	= dra_deg * floor(ra_boxes_deg/dra_deg)
  dec_min_deg	= ddec_deg * floor(dec_boxes_deg/ddec_deg)
  ra_max_deg	= ra_min_deg + dra_deg
  dec_max_deg	= dec_min_deg + ddec_deg

  nboxes	= sizex(ra_boxes_deg)

  subcat_str	= strarr(nboxes)

  for i = 0L, nboxes-1 do begin
    subcat_str[i]	= strcompress( string(catalog, ra_min_deg[i], $
                          ra_max_deg[i], dec_min_deg[i], dec_max_deg[i], $
                          format = $
	                 '(A, "_", F5.1, "_", F5.1, "_", F5.1, "_", F5.1)'), $
                          /remove_all)
  end

  subcat_str	= subcat_str[uniq(subcat_str)]

  nboxes	= sizex(subcat_str)

  if keyword_set(VERBOSE) then print, 'Grid size = ' + st(dra_deg) + ' deg'
  if keyword_set(VERBOSE) then print, 'Search radius = ' + $
                                       st(radius_deg) + ' deg'
  if keyword_set(VERBOSE) then print, 'Number of boxes = ' + st(nboxes)

; Make sure that each box exists; if not, create it

  for i = 0L, nboxes-1 do begin
    if not(oc_subcat_exists(subcat_str[i])) then begin
    
      oc_subcat_create, subcat_str[i]
      if keyword_set(VERBOSE) then print, $
             'File created: ' + subcat_str[i] + ' ' + st(i) + '/' + st(nboxes)
    end $
    else begin
      if keyword_set(VERBOSE) then print, $
             'File exists: ' + subcat_str[i]
    end
  end

; Create the temporary catalog file, which has from 1-4 boxes in it
; First create the catalog header

  a		= strarr(2)

;  file_tmp	= 'search_' + st(get_pid()) + '.tmp'
;  path_tmp	= getenv('HDOCCUL_TMP') + '/' + file_tmp
  path_tmp = tmpfile('search_', 'tmp')
  fdecomp, path_tmp,disk,dir,name,qual
  file_tmp = name+'.'+qual
  a[0]		= file_tmp + $		; Catalog name
    			'/d' + $	; Decimal degrees
    			'/j' + $	; J2000
    			'/2' + $	; Two magnitudes
    			'/s' + $	; Spectral type
			'/i' + $	; Ignore everything else 
; 			     (i.e., the distance in arcsec from previous scat)
			''
  a[1]		= 'Created with oc_search_pos_subcat.pro, ' + systime()
  write_txt_file, path_tmp, a

; Now populate it with stars from all of the appropriate boxes

  command = 'cat '

; Pipe concatenation through sort and uniq to remove the occasional
; duplicate star

  for i = 0L, nboxes-1 do command = command + (path_subcat+subcat_str[i] + ' ')
  command 	= command + ' | sort | uniq >> ' + path_tmp 
  spawn, command
  if keyword_set(VERBOSE) then print, 'Created: ' + path_tmp

; Now search the temporary file

; Construct the command

  command	= 'scat ' + $ 
; catalog: temporary file
                  ' -c ' + path_tmp + $			
; search radius, arcsec
		  ' -r ' + string(radius_arcsec,format='(F20.5)') + $	
; degrees (not hms)
                  ' -d' + $				
; tab-delimited columns (also adds 11 lines of header info)
                  ' -t' + $				
; J2000
                  ' -j' + $				
; Max number returned
		  ' -n ' + st(nstars_max) + $		
		  ' '    + string(ra_deg,format='(F19.12)') + ' ' + $
		  string(dec_deg,format='(F19.12)') + ' j2000'

; Execute the command and return output to array 'scat'

  if keyword_set(VERBOSE) then print, 'Executing command: ' + command
  spawn, command, scat

; Remove the temporary file

  command_rm = 'rm ' + path_tmp
  spawn, command_rm

; Now read the output from scat command

  case catalog of
    'hd' : lines_header_table	= 11
  end

  if grep(command, '-t') then lines_header	= lines_header_table $
    else		      lines_header	= 0
  line_headings	= lines_header-2		; index of the headings line
  nstars	= n_elements(scat)-lines_header

; If none found, then return

  if (nstars eq 0) then return, -1

; If stars found, then prepare to read them in

  scat_chop	= scat[lines_header:*]

; Parse and store the scat output Scan will not output which has fixed column
; locations; therefore, IDL's FORMAT= routines do not work.  We therefore use
; tab-delimited output, and split on tabs, which prevents problems with ID's
; containing spaces.

  char_tab	= string(9B)
  ids		= strarr(nstars)

for i = 0L, nstars-1 do begin
    ids[i]	= (ht_str_split(scat_chop[i], char=char_tab))[0]
  end

; 'scat' uses an abbreviated catalog.  We now take the ID's, and use these to
; look up the full original catalog line.

  case strlowcase(catalog) of
    'hd' : begin
      ids	= long(ids)		; For HD: convert strings to integers
      stars	= oc_getstar_hd(ids, count=count, verbose=verbose)
    end
  end

; Return the entire populated structures.
; Remove any duplicates.  We don't expect there would be duplicate stars here,
; but this is just to prevent it.

  ids		= stars[*].id
  return, stars[uniq(ids)]

  return, stars
end