;+
; NAME:
;   GTIMERGE
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   Merge two Good Time Interval (GTIs) arrays into a single array
;
; CALLING SEQUENCE:
;   NEWGTI = GTIMERGE(GTI1, GTI2, COUNT=COUNT, [/INTERSECT, /UNION,
;                     /INVERT1, /INVERT2, TTOLERANCE=])
;
; DESCRIPTION: 
;
;   The function GTIMERGE accepts two existing valid Good Time
;   Interval (GTI) arrays and merges them into a single array.  Either
;   the intersection or the union of the two GTIs are returned.
;
;   The intersection refers to the set of intervals which lie in both
;   intervals.  The union refers to the set of intervals which lie in
;   at least one but not necessarily both.  Here is an example of both
;   kinds of operations.  Let us start with two intervals here:
;
;        0              50              100                  170   GTI1
;   <----|==============|----------------|====================|------>
;
;                 30                         120                   GTI2
;   <--------------|==========================|---------------------->
;
;   These intervals would be represented by GTI1=[[0,50],[100,170]]
;   and GTI2=[[30,120]]. The intersection of the two sets of intervals
;   are the points which lie in both, ie [[30,50],[100,120]]:
;
;                 30   50               100  120              INTERSECT
;   <--------------|====|----------------|====|---------------------->
;
;   The union is the combination of both intervals, ie [[0,170]]:
;
;        0                                                   170  UNION
;   <----|====================================================|------>
;
;   It is also possible to treat either one of the input arrays as
;   "bad" intervals using the INVERT1 and/or INVERT2 keywords.  When
;   an interval is inverted, then the output is composed only of areas
;   *outside* the specified intervals.
;
;   It should be noted that this function is not constrained to
;   operation only on time arrays.  It should work on any
;   one-dimensional quantity with intervals.
;
;
; PERFORMANCE:  Combining many intervals
;
;   Users who wish to combine many intervals in sequence will find a
;   performance degradation.  The problem is that each GTIMERGE
;   operation is order N^2 execution time where N is the number of
;   intervals.  Thus, if N mostly distinct GTIs are merged, then the
;   running time will be order N^3.  This is unacceptable, but there
;   is a workaround.
;
;   Users can accumulate "sub" GTIs by merging subsets of the full
;   number of intervals to be merged, and then occasionally merging
;   into the final output GTI.  As an example, here first is a simple
;   merging of 1000 different GTIs:
;
;     totgti = 0L                   ;; Empty GTI
;     FOR i = 0, 999 DO BEGIN
;         gti = ...                
;         totgti = gtimerge(totgti, gti, /union)
;     ENDFOR
;
;   This computation may take a long time.  Instead the merging can be
;   broken into chunks.
;
;     totgti = 0L
;     chgti  = 0L      ;; "Chunk" GTI
;     FOR i = 0, 999 DO BEGIN
;         gti = ...
;         chgti = gtimerge(chgti, gti, /union)
;         if (n_elements(chgti) GT 100) OR (i EQ 999) then begin
;             ;; Merge "chunk" gti into final one, and reset
;             totgti = gtimerge(totgti, chgti, /union)
;             chgti = 0L
;         endif
;     ENDFOR
;
;   Note that the final merge is guaranteed because of the (i EQ 999)
;   comparison.
;
; INPUTS:
;
;   GTI1, GTI2 - the two input GTI arrays.
;
;         Each array is a 2xNINTERVAL array where NINTERVAL is the
;         number of intervals, which can be different for each array.
;         GTI(*,i) represents the start and stop times of interval
;         number i.  The intervals must be non-overlapping and
;         time-ordered (use GTITRIM to achieve this).
;
;         A scalar value of zero indicates that the GTI is empty, ie,
;         there are no good intervals.
;
; KEYWORDS:
;
;   INTERSECT - if set, then the resulting GTI contains only those
;               intervals that are in both input sets.
;
;   UNION - if set, then the resulting GTI contains those intervals
;           that are in either input set.
;
;   COUNT - upon return, the number of resulting intervals.  A value
;           of zero indicates no good time intervals.
;
;   INVERT1 - if set, then GTI1 is considered to be inverted, ie, a
;             set of "bad" intervals rather than good.
;
;   INVERT2 - if set, then GTI2 is considered to be inverted, ie, a
;             set of "bad" intervals rather than good.
;
;   TTOLERANCE - a scalar value indicating the tolerance for
;                determining whether values are equal.  This number
;                can be important for intervals that do not match
;                precisely.
;                Default: Machine precision
;
; RETURNS:
;
;   A new GTI array containing the merged 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:
;
;   GTITRIM, GTIENLARGE
;
; MODIFICATION HISTORY:
;   Written, CM, 1997-2001
;   Documented, CM, Apr 2001
;   Handle case of zero-time GTIs, CM, 02 Aug 2001
;   Handle "fractured" GTIs correctly, though worriedly, CM, 15 Oct
;     2001
;   Handle case where both inputs are empty, but /INVERT1 and/or
;     /INVERT2 are set, CM, 08 Aug 2006
;
;  $Id: gtimerge.pro,v 1.7 2006/10/22 09:49:53 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.
;-
;      print, '  ', o1(uu(j):uu(j+1)-1), format='(A0,100(D4.0," ",:),$)'
;      print, ' --> ', oo1
;      print, '  ', o2(uu(j):uu(j+1)-1), format='(A0,100(D4.0," ",:),$)'
;      print, ' --> ', oo2

      ;; yes1 and yes2 are updated here, if there is a change in gti state
      if oo1 GE 0 then yes1 = (oo1 EQ 1)
      if oo2 GE 0 then yes2 = (oo2 EQ 1)

      ;; Here are the conditions for starting a new GTI entry
      if (NOT within AND (yes1 OR  yes2) AND keyword_set(union)) OR $
        ( NOT within AND (yes1 AND yes2) AND keyword_set(intersect)) then $
        begin
          tstart = tt(i)
          within = 1
      endif
      if dbl1 AND yes1 then yes1 = 0
      if dbl2 AND yes2 then yes1 = 0

      ;; And the conditions for ending it.
      if (within AND (NOT yes1 AND NOT yes2) AND keyword_set(union)) OR $
        ( within AND (NOT yes1 OR  NOT yes2) AND keyword_set(intersect)) then $
        begin
          if numgti EQ 0 then $
            newgti = [ tstart, tt(i) ] $
          else $
            newgti = [ newgti, tstart, tt(i) ]
          numgti   = numgti + 1
          within   = 0

          if dbl1 OR dbl2 then begin
              ;; Now restart the GTI if desired
              tstart = tt(i)
              within = 1
              dbl1 = 0 & dbl2 = 0
          endif
      endif

      LOOPCONT:
  endfor

  ;; I think this clause is activated only if the INVERT1 or INVERT2
  ;; keywords are activated.  This should be tested.  Boundary
  ;; conditions can be wierd.
  if within then begin
      if numgti EQ 0 then $
        newgti = [ tstart, max(tt) ] $
      else $
        newgti = [ newgti, tstart, max(tt) ]
      numgti   = numgti + 1
      within   = 0
  endif
  
  ;; Empty GTI is a special case.  Here I have opted to *not* take the
  ;; lame HEASARC route.  Instead, I define a single element GTI array
  ;; to mean no good times.
  if numgti EQ 0 then begin
      EMPTY_GTI:
      count = 0L
      return, 0L
  endif

  ;; Reconstitute the new GTI array
  count    = numgti
  newgti   = reform(newgti, 2, numgti, /overwrite)

  ;; Now do some clean-up, removing zero-time GTIs and GTIs that
  ;; adjoin.  This really should not be necessary, since by definition
  ;; the above technique should not create buggy GTIs.
  if numgti GT 0 then begin
      ocount = count
      newgti = gtitrim(newgti, maxgap=ttol, count=count)
  endif

  return, newgti
end