; This IDL procedure should (very nearly) reproduce
; Figure 2 of McCullough et al 2006, ApJ, 648, 1228.
; (also http://xxx.lanl.gov/abs/astro-ph/0605414)

pro xo_dr1

  ref = '240p31'

  ; Reminder: IDL counts from 0, whereas many FITS implementations begin at 1.
  ; For example, 634 here should be 635 if one uses ds9 to display the FITS files.

  id = 634	; XO-1

  ; HJD0, Period, RA (degrees), DEC (degrees) [2000]

  ephem_xo1 = [2453808.9170d0,3.941534d0,240.549334502d0,28.1695641682d0]

  ; keyword file

  key  = readfits(ref+'key.fits',/silent)

  ; Extract heliocentric julian date for nominal center of field of view.
  ; hjd = Heliocentric JD - 2450000.0

  hjd  = key[9,*]

  print,'Julian date range (-2450000.):',min(hjd),max(hjd)

  ; calf = calibrated and flagged magnitudes.

  calf = readfits(ref+'calf.fits',/silent)

  ; Compute phase.
  ph    = ((hjd - (ephem_xo1[0]-245.d4))/ephem_xo1[1]) mod 1.d0
  k     = where(ph lt -0.5d0)
  ph[k] = ph[k] + 1.0d0

  ; Convert from differential magnitudes to flux units.

  flux = 10.^(-0.4*calf[id,*])

  ; Make the plot.

  plot,ph*ephem_xo1[1],flux,psy=3,xr=[-2.,2.],yr=[0.96,1.04],xtit='Time from mid-transit [days]',ytit='Relative Flux',charsi=2.

return
end