;+
; NAME:
;  oc_getstar_tyc2_test.pro
;
; PURPOSE: (one line)
;  Validates the oc_getstar_tyc2 routine by searching for Tycho-2 stars.
;
; DESCRIPTION:
;  Searches for a large number of stars, using made-up ID's.  
;  Some of these will be valid; most will be invalid.  The 
;  results are printed, along with 'benchmark' results 
;  to compare to to validate proper functionality.
;
;  Then, prints a series of diagnostics to verify that proper motion
;  is applied properly.
; 
; CATEGORY:
;  Star catalogs
;  
; CALLING SEQUENCE:
;  oc_getstar_tyc2_test
;
; INPUTS:
;  None
;
; OPTIONAL INPUT PARAMETERS:
;  None
;  
; KEYWORD INPUT PARAMETERS:
;  None
;
; KEYWORD OUTPUT PARAMETERS:
;  None
;
; OUTPUTS:
;  None
;  
; COMMON BLOCKS:
;  None
;
; SIDE EFFECTS:
;  None
;
; RESTRICTIONS:
;  Tycho-2 catalog must be available at $WCS_CATDIR/tyc2/catalog.dat .
;
; EXAMPLE:
;  oc_getstar_tyc2_test
;   [Makes a plot, and prints results along with benchmark results]
; 
; MODIFICATION HISTORY:
;  Written   2-Nov-2005 by Henry Throop, SwRI.
;  Modified 28-Feb-2006 by HBT. Improved documentation and formatting.
;  Modified 24-Feb-2006 by HBT. Added testing of proper motion fields.
;-
pro oc_getstar_tyc2_test

; First, make up a list of TYC2 id's.  Some of these will be valid and some
; will not be; thus, the TYC2 searching is fully tested for error handling.

  sector_1	= '11'
  sector_2	= '15'
  num_search	= 500

  ids		= strarr(num_search)
  for i = 0, num_search-1 do ids[i] = $
    (i le num_search/2 ? sector_1 : sector_2) + '-' + st(i)

; Look up their coords based on ID

  stars		= oc_getstar_tyc2(ids, NSTARS=nstars)

; Get their RA/Dec and plot them

  ra		= stars.ra
  dec		= stars.dec

  plot, ra, dec, psym=2, xtitle = 'RA [rad]', ytitle = 'Dec [rad]', $
    xrange=[0.15, 0.29], yrange=[0, 0.09]

; Print statistics

  print, 'Stars searched: ' + st(num_search)
  print, 'Stars found: ' + st(nstars)
  print, 'Total RA: ' + st(total(ra > 0))
  print, 'Total Dec: ' + st(total(dec > 0))

  print
  print, '  Benchmark results HBT 2-Nov-2005:'
  print, '  Stars searched: ' + st(500)
  print, '  Stars found: ' + st(74)
  print, '  Total RA: ' + st(103.84328)
  print, '  Total Dec: ' + st(22.115512)

;;; Now, test the proper motions and ET conversions
; Check J2000 coordinates.  Epoch 2000, Equinox 2000.

  id		= ['11-49-1']
  yr 		= cspice_jyear()
  d2r             = !dpi / 180.d  ; multiply by to convert degrees to radians
  r2d             = 180.d / !dpi  ; multiply by to convert radians to degrees
  d2as            = 60d*60d       ; multiply by to convert degrees to accsec
  r2as            = r2d * d2as    ; multiply by to convert radians to arcsec
  as2r            = 1/r2as        ; multiply by to convert arcsec to radians

  star = oc_getstar_tyc2(id)
  print
  print, 'Performing motion/precession tests on TYC2 ', id
  print, 'Position from star.ra'
  print, 'RA J2000: ' + radtohms_str(star.ra)
  print, 'DE J2000: ' + radtodms_str(star.dec)
  print

  print, 'Now applying PM to move from catalog ET to ET=2000 yr'
  star = oc_star_apply_pm(star, (2000-2000)*yr)
  print, 'RA J2000: ' + radtohms_str(star.ra)
  print, 'DE J2000: ' + radtodms_str(star.dec)
  print
 
  print, 'Vizier: TYC2 11-49-1: RA J2000 = 00 40 51.989'
  print, 'Vizier: TYC2 11-49-1: DE J2000 =+02 28 15.67'
  print

; Check the proper motion from 2000 to 1950, using oc_star_apply, with no
; equinox change

  print, 'Now applying proper motion to 1950, with no precession'
  star = oc_getstar_tyc2(id)
  star = oc_star_apply_pm(star, (1950-2000)*yr)
  print, 'RA 2000: ' + radtohms_str(star.ra) + ' +- ' + $
    radtohms_str(star.raerr)
  print, 'DE 2000: ' + radtodms_str(star.dec) + ' +- ' + $
    radtodms_str(star.decerr)
  print, 'RA 2000: ' + radtohms_str(star.ra) + $
    ' Epoch 1950, Equinox 2000, PM manually calculated'
  print, 'DE 2000: ' + radtodms_str(star.dec) + $
    ' Epoch 1950, Equinox 2000, PM manually calculated'
  print

; Check the precession alone, with no proper motion
  print, 'Now applying precession to B1950, with no proper motion'
  star = oc_getstar_tyc2(id)
  bprecess, star.ra*r2d, star.dec*r2d, ra_1950_deg, dec_1950_deg
  print, 'RA 1950: ' + radtohms_str(ra_1950_deg*d2r) + $
    ' +- ' + radtohms_str(star.raerr)
  print, 'DE 1950: ' + radtodms_str(dec_1950_deg*d2r) + $
    ' +- ' + radtodms_str(star.decerr)


; Check the proper motion from 2000 to 1950, using manual application, with no
; equinox change

  star = oc_getstar_tyc2(id)
  star.ra	= star.ra + (1950-2000)*yr*star.radot/cos(star.dec)
  star.dec= star.dec+ (1950-2000)*yr*star.decdot
  print

; Check B1950 coordinates.  Epoch 1950, Equinox 1950.  Can be compared to
; Vizier.  Q: does the order of application matter?  I assume first PM, then
; precess?  It does substantially matter which way we do it in.  PM must be
; defined in some coord frame -- if it's J2000, need to use J2000. *** This
; should be a new field -- the equinox of PM!

  print, 'Now applying proper motion and precession to 1950'
  star = oc_getstar_tyc2(id)
  star = oc_star_apply_pm(star, (1950-2000)*yr)
  bprecess, star.ra*r2d, star.dec*r2d, ra_1950_deg, dec_1950_deg

  print, 'Proper motion per year: ' + radtohms_str(star.radot*yr)
  print, 'Proper motion per year: ' + radtodms_str(star.decdot*yr)
  print, 'RA 1950: ' + radtohms_str(ra_1950_deg*d2r) + ' +- ' + $
    radtohms_str(star.raerr)
  print, 'DE 1950: ' + radtodms_str(dec_1950_deg*d2r) + ' +- ' + $
    radtodms_str(star.decerr)
  print
 
  print, 'Now printing Vizier catalog positions, Epoch 1950, Equinox 1950'
  print, 'Vizier: TYC2 11-49-1: RA 1950 = 00 38 17.770'
  print, 'Vizier: TYC2 11-49-1: DE 1950 =+02 11 48.37'

; Print position errors after this

; Print PM rates

  print, 'pmRA: ' + st(star.radot*yr*r2as*1000)  + ' mas/yr'
  print, 'pmDE: ' + st(star.decdot*yr*r2as*1000) + ' mas/yr'

  stop

end