;+
; NAME:
;   GTIREAD
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   Read a Good Time Interval (GTI) from a FITS input file.
;
; CALLING SEQUENCE:
;   GTI = GTIREAD(INPUT, EXTENSION, [COUNT=, INTERSECT=, UNION=])
;
; DESCRIPTION: 
;
;   The function GTIREAD reads a set of good time intervals from one
;   or more input FITS files.  By default the FITS extension number is
;   1, however the user can specify a different extension with the
;   EXTENSION keyword.
;
;   The GTIs must be stored in the standard OGIP FITS format for good
;   time intervals.  An additional extension keyword named EMPTYGTI is
;   allowed.  When this keyword is present in the header, the GTI is
;   assumed to be empty and no good intervals are generated.
;
;   As the degenerate case, it is acceptable for the user to pass an
;   existing GTI to GTIREAD via the INPUT parameter.  If this is the
;   case then the same GTI is returned.  Thus, a procedure can
;   overload a GTI-like keyword to allow the user to pass either the
;   data or a filename.
;
;   If multiple input files are specified, then the user has a choice
;   how to join them together.  The INTERSECT keyword will compute the
;   intersection of all files (ie, only those intervals which are
;   shared by all files); the UNION keyword will compute the union of
;   all files (ie, those intervals which appear at least once in any
;   file).
;
; INPUTS:
;
;   INLIST - either an existing valid good time interval, as a
;            numerical array, or a string array containing a list of
;            FITS file names.
;
;            If an array of numbers are supplied, then they are
;            returned immediately to the user.
;
;            If an array of strings are supplied, then each file is
;            opened in turn and joined together as described above.
;            The "batch" style of GTI is allowed (using the
;            "@gtilist.txt" format to indicate a list of GTI files is
;            stored in gtilist.txt).
;
;   EXTENSION - a optional scalar FITS extension number.
;               Default: 1
;
; KEYWORDS:
;
;   COUNT - upon return, the number of resulting good intervals.  A
;           value of zero indicates no good intervals were found.
;
;   INTERSECT - if set, then the GTIs of multiple files are joined
;               together using their intersection.
;
;   UNION - if set, then the GTIs of multiple files are joined
;           together using their union.
;
;   STATUS - upon return, this keyword is set to 1 if at least one GTI
;            file was successfully read, or 0 if none were successful.
;
;   QUIET - if set, then no error messages will be printed.
;
; RETURNS:
;
;   A new GTI array containing the resulting intervals.  The array is
;   2xCOUNT where COUNT is the number of resulting intervals.
;   GTI(*,i) represents the start and stop times of interval number i.
;   The intervals are non-overlapping and time-ordered.
;
;   If COUNT is zero then the returned array is a scalar value of
;   zero, indicating no good intervals were found.
;
; SEE ALSO:
;
;   GTIMERGE, GTITRIM
;
; MODIFICATION HISTORY:
;   Written, CM, 1997-2001
;   Documented, CM, Apr 2001
;
;  $Id: gtiread.pro,v 1.4 2002/04/12 09:58:28 craigm Exp $
;
;-
; Copyright (C) 1997-2001, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
function gtiread, inlist, extension, count=ngti, $
                  intersect=intersect, union=union, $
                  timezero=timezero, mjdref=mjdref, toffset=toffset, $
                  quiet=quiet, status=status

  if n_elements(extension) EQ 0 then extension = 1
  errmsg = ''
  ngti = 0L
  status = 0

  sz = size(inlist)
  ;; Numbers that are passed in are returned directly
  if sz(sz(0)+1) NE 7 then begin
      if n_elements(inlist) EQ 1 then begin
          count = 0L
          return, 0
      endif
      ngti = n_elements(inlist)/2
      status = 1
      return, reform(inlist, 2, ngti)
  endif

  ;; Default is INTERSECT
  intersect = keyword_set(intersect)
  union     = keyword_set(union)
  if (union EQ 0) AND (intersect EQ 0) then intersect = 1
  if (union EQ 1) AND (intersect EQ 1) then begin
      if NOT keyword_set(quiet) then begin
          message, 'WARNING: cannot perform both UNION and INTERSECT.', /info
          message, '         INTERSECT assumed.', /info
      endif
      union = 0
  endif

  ;; Look for "batch" files, which list other files
  batchscan, inlist

  first_t = 1
  for i = 0L, n_elements(inlist)-1 do begin
      filename = inlist(i)
      fxbopen, un, filename, extension, header, errmsg=errmsg
      if errmsg NE '' then begin
          if NOT keyword_set(quiet) then $
            message, 'WARNING: could not open GTI file '+filename, /info
          goto, NEXT_FILE
      endif
      
      ingti = fxpar(header, 'NAXIS2')
      if ingti LE 0 then goto, CLOSE_FILE
      ggstart = 0D
      emptygti = fxpar(header, 'EMPTYGTI')
      if emptygti then goto, EMPTY_GTI

      errmsg = ''
      fxbread, un, ggstart, 'START', errmsg=errmsg
      if errmsg EQ '' then fxbread, un, ggstop,  'STOP', errmsg=errmsg
      if errmsg NE '' then begin
          if NOT keyword_set(quiet) then $
            message, 'WARNING: could not read GTI file '+filename, /info
          goto, CLOSE_FILE
      endif
      
      !err = 0
      itimezeri = double(fxpar(header, 'TIMEZERI'))
      itimezerf = double(fxpar(header, 'TIMEZERF'))
      if !err LT 0 then itimezero = double(fxpar(header, 'TIMEZERO')) $
      else              itimezero = itimezeri + itimezerf
      !err = 0

      imjdrefi = double(fxpar(header, 'MJDREFI'))
      imjdreff = double(fxpar(header, 'MJDREFF'))
      if !err LT 0 then imjdref = double(fxpar(header, 'MJDREF')) $
        else            imjdref = imjdrefi + imjdreff
      !err = 0
      
      if first_t then begin
          first_t = 0
          timezero = itimezero
          mjdref   = imjdref
      endif else if keyword_set(toffset) then begin

          ;; If TOFFSET keyword is set, then offset the values to a
          ;; common TIMEZERO value.  Often this is not desired, if
          ;; everything is in spacecraft time to begin with.
          ggstart = ggstart + (itimezero - timezero)
          ggstop  = ggstop  + (itimezero - timezero)
      endif
          
      ingti = n_elements(ggstart)
      igti = make_array(2, ingti, value=ggstart(0)*0)
      igti = reform(igti, 2, ingti, /overwrite)

      igti(0,*) = ggstart
      igti(1,*) = ggstop
      ;; Special case of an empty GTI
      if n_elements(ggstart) EQ 1 AND n_elements(ggstop) EQ 1 then begin
          if ggstart(0) EQ 0 AND ggstop(0) EQ 0 then begin
              EMPTY_GTI:
              igti = ggstart(0)*0
              ingti = 0L
          endif
      endif

      if ngti EQ 0 then begin
          gti  = igti
          ngti = ingti
      endif else begin
          gti = gtimerge(gti, igti, count=ngti, $
                         intersect=intersect, union=union)
      endelse

      CLOSE_FILE:
      fxbclose, un
      NEXT_FILE:
  endfor

  if first_t EQ 0 then status = 1

  if ngti GT 0 then begin
      gti = reform(gti, 2, ngti, /overwrite)
      return, gti
  endif else begin
      return, 0L
  endelse

end