pro xrt_ecid,year,month,day,miss,idx,tstamp,miss_pack,tot_pack,zidx,file,obevfile=obevfile,stop_default=stop_default, ql=ql
;+
; PROJECT: 
;       Solar-B / Hinode / XRT
; 
; NAME: 
;       XRT_ECID
;
; CATEGORY:
;
;       XRT Data Verification
; 
; PURPOSE:
;
;       Identify missing files or missing packets in files for a given
;       date. 
;
; CALLING SEQUENCE:
; 
;       XRT_ECID, YEAR, MONTH, DAY, MISS_FILE, MISS_INDEX, $
;       TIME_STAMP, MISS_PACKET, TOTAL_PACKET, MISS_PACKET_INDEX, FILENAME $
;       [,OBEVFILE=OBEVFILE] [,/STOP_DEFAULT]
;
; INPUTS:
;
;       YEAR           - Year (e.g., 2008)
;
;       MONTH          - Month (e.g.,   5)
;
;       DAY            - Day  (e.g.,   12)
;
;       OBEVFILE       - [Optional] (String; Filename) Specify the
;                        OBEV file that covers the given day. This
;                        file is used to identify SAA_ENTRY and
;                        NGT_ENTRY times (i.g., when CTRL_MANU is
;                        issued). 
;
; KEYWORDS:
;       
;       /STOP_DEFAULT  - [Optional] (Boolean) Use the default
;                        negative time-delay duration. By default, it
;                        uses SAA_ENTRY = -0s, SAA_EXIT = +60s, NGT_ENTRY =
;                        -480s and NGT_EXIT = + 480s delays. 
;
;                        Currently this flag is meaningless. 
;
;       /QL            - [Optional] Switch to change the rootdir
;                        pointer to designated QLraw directory. 
;
;
; OUTPUTS:
;       
;       MISS_FILE      - Total number of missing files identified. 
; 
;       MISS_INDEX     - Identifier index that may be useful to find
;                        out where missing files occurred. Also useful
;                        to find out where negative or identical
;                        EC_IDs have occurred (if any). 
;
;       TIME_STAMP     - Time stap (DATE_OBS) for each file examined. 
; 
;       MISS_PACKET    - Total number of missing packets in unit of 64
;                        pixels (one minimal packet). If this
;                        number isn't a whole number, then the
;                        users are encouraged to examine the
;                        MISS_PACKET_INDEX to see if any files are
;                        corrupted. 
; 
;       TOTAL_PACKET   - Total of the packet size expected (64 pixels
;                        per unit)
; 
;       MISS_PACEKT_INDEX - Similar to MISS_INDEX, contains how many
;                           packets (again in unit of 64pixels) are
;                           missing. 
; 
;       FILENAME       - Names of FITS files examined. 
;       
; EXAMPLES:
;
;       There are two ways to use this tool. 
;   
;       Identify missing files and packets without obev file. If XRTCO
;       is not taking data during SAA and NGT, then this will result
;       in spurious missing files. 
;
;       IDL> xrt_ecid, 2007, 3, 29, miss, idx, tstamp, missp, zidx, filename 
;
;
;       Identify missing files and packets using an obev file. This is
;       recommended if the user needs an accurate account for missing
;       files. 
; 
;       IDL> xrt_ecid, 2008, 5, 12, miss, idx, tstamp, missp, zidx, $
;       filename, obevfile='obev_200805110100.evt', /stop_default
; 
;
; NOTES:
;
;       !!CAVEAT!!  
;       Each user must modify "rootdir" variable to the
;       right location for Hinode/XRT data archive directory. 
;
;      
;       (1) MISS_INDEX: The index can have the following values. 
;             0        = No anomaly (or gap)
;             1 ... N  = One or ...N files missing between index[i-1]
;                        and index[i]. 
;            -1        = Two identical EC_ID detected at this index. 
;            -9999     = Negative EC_ID detected at this index. 
; 
;       (2) See more gory details on how this algorithm works in the
;           code. 
;
; HISTORY: 
; 
;       K. Ishibashi   21-May-2008  Original
;       K. Ishibashi   17-Jun-2008  Bug fixes, introducing a better
;                                   version control 
;       K. Ishibashi   19-Jun-2008  Add switch for rootdir pointer
;                                   from level0 and QLraw data dir. 
;       K. Ishibashi   16-Jul-2008  Add SAA_EXIT, NGT_EXIT to make 
;                                   options more explicit. 
;       K. Ishibashi   23-Jul-2008  Add one output parameter (total
;                                   packet numbers) RFE M. Weber. 
;       K. Ishibashi   24-Jul-2008  Modify logic such that the script
;                                   will always identify the last
;                                   previous file. 
;       K. Reeves      08-Sep-2010  fixed bug for ec_id reset in case
;                                   of safehold
;       K .Reeves      09-Dec-2010  changed logic for ECID wrapping to
;                                   account for case where there are 
;                                   missing immages around the time of
;                                   the wrap.
;
;-


; Version control
;
; Ver = x.y.z 
;             where x == main version control
;                   y == significant update
;                   z == minor bug fixes
  ver = '0.5.0'                                      
  ver = '0.5.1'  ;; added logic for negative last_ecid value (manual exposure last image on previous day
  ver = '0.5.2'  ;; fixed bug for ec_id reset in case of safehold
  ver = '0.5.3'  ;; changed logic for ec_id wrap
  
;----  change rootdir as you see fit (this is the hardcoded variable) ----- 
;
;  Be sure to have the "/" at the end of string. 
;
;  Introducing /QL switch. 


  if keyword_set(ql) then begin
     rootdir = '/archive/hinode/xrt/QLraw/' ; SAO
;     rootdir = '/soda/solarb/xrt/QLraw/'           ;ISAS
  endif else begin
     rootdir = '/archive/hinode/xrt/level0/' ; SAO
;  rootdir = '/data/cora4/hinode/xrt/level0/'    ; NWRA
;  rootdir = '/soda/hinode/xrt/level0/'          ; ISAS
  endelse


;-- read obev file if requested -- 
; 
; Parts of this algorithm were taken from xrt_plan.pro written by Narukage. 
;


  if keyword_set(obevfile) then begin

; remove leading and trailing white space from the string obevfile
     obevfile = strtrim(obevfile,2) ;

     if file_test(obevfile) then begin
        obev = {event: '', time: '', tsec:0d0}
        get_lun, unit
        openr,unit,obevfile
        temp = ''               ; so that readf reads in variables as string. 
        
        while ~ EOF(unit) do begin
           readf, unit, temp
           if strmid(temp, 0, 9) eq 'NGT_ENTRY' then obev = [obev, {event: 'NGT_ENTRY', time: STRJOIN(STRSPLIT( STRMID( STRJOIN(STRSPLIT(temp, ' ;', /EXTRACT)) , 18, 19, /REVERSE_OFFSET) , '.', /EXTRACT), ' '), tsec:0d0}]
           if strmid(temp, 0, 8) eq 'NGT_EXIT' then obev = [obev, {event: 'NGT_EXIT', time: STRJOIN(STRSPLIT( STRMID( STRJOIN(STRSPLIT(temp, ' ;', /EXTRACT)) , 18, 19, /REVERSE_OFFSET) , '.', /EXTRACT), ' '), tsec:0d0}]

           if strmid(temp, 0, 9) eq 'SAA_ENTRY' then obev = [obev, {event: 'SAA_ENTRY', time: STRJOIN(STRSPLIT( STRMID( STRJOIN(STRSPLIT(temp, ' ;', /EXTRACT)) , 18, 19, /REVERSE_OFFSET) , '.', /EXTRACT), ' '), tsec:0d0}]
           if strmid(temp, 0, 8) eq 'SAA_EXIT' then obev = [obev, {event: 'SAA_EXIT', time: STRJOIN(STRSPLIT( STRMID( STRJOIN(STRSPLIT(temp, ' ;', /EXTRACT)) , 18, 19, /REVERSE_OFFSET) , '.', /EXTRACT), ' '), tsec:0d0}]
        endwhile
        free_lun, unit          ; this will close the obevfile. No loose end. 
      
; populate the second elapse time into obev.tsec. 
        obev.tsec = anytim(obev.time)

; Time correction for CTRL_MANU
; Assume the following:
;        SAA_ENTRY -- No time delay
;        NGT_ENTRY -- -8min (-480s) time delay. 
;
; As of right now, the default delay values are always used. 
; The code may be changed to accept user-specific values in future. 
        if keyword_set(stop_default) then begin
           SAA_ENTRY = 0d0  
           SAA_EXIT  = 60d0
           NGT_ENTRY = 480d0
           NGT_EXIT  = 480d0
        endif else begin
           print,'obev: forcing to use stop_default...' 
           SAA_ENTRY = 0d0  
           SAA_EXIT  = 60d0
           NGT_ENTRY = 480d0
           NGT_EXIT  = 480d0
        endelse 
        
        for i = 0, n_elements(obev.tsec) - 1 do  begin
           if (obev[i].event eq 'SAA_ENTRY') then obev[i].tsec = obev[i].tsec - SAA_ENTRY
           if (obev[i].event eq 'SAA_EXIT')  then obev[i].tsec = obev[i].tsec + SAA_EXIT
           if (obev[i].event eq 'NGT_ENTRY') then obev[i].tsec = obev[i].tsec - NGT_ENTRY
           if (obev[i].event eq 'NGT_EXIT')  then obev[i].tsec = obev[i].tsec + NGT_EXIT
        endfor
      
; Select events within the specified day
        daystrt = anytim(strtrim(year,1)+'/'+strtrim(month,1)+'/'+strtrim(day,1)+' 00:00:00')
        daystop = anytim(strtrim(year,1)+'/'+strtrim(month,1)+'/'+strtrim(day,1)+' 23:59:59.59')
        h       = where((obev.tsec lt daystrt) or (obev.tsec gt daystop), hcnt)
        
        if (hcnt eq 0) then print,'obev: nothing to be removed (this may be OK)'
        if (hcnt eq n_elements(obev)) then begin
           print,'obev: no valid entryin obevfile... aborting.'
           return
        endif else begin 
           remove,h,obev
        endelse      
     endif else begin
        print,'obev: no valid obevfile found...aborting...'
        return
     endelse
  endif

;---- Locate EC_ID used in the last file taken on the previous day  ----  

  mm    = month
  dd    = day-1
  yearx = year
  if (month lt 10) then mm = '0' + strtrim(month,1)
  if (day lt 11)   then dd = '0' + strtrim(day-1,1)


; Avoid using Solarsoft when possible. 
  if (day eq 1) then begin
     if (month eq 11) then begin
        dd = 31
        mm = 10
     endif
     if ((month eq 2) or (month eq 4) or (month eq 6) or (month eq 8) or (month eq 9)) then begin
        dd = 31
        mm = '0' + strtrim(month - 1,1) 
     endif
     if (month eq 12) then begin
        dd = 30
        mm = month - 1
     endif
     if ((month eq 5) or (month eq 7) or (month eq 10)) then begin
        dd = 30 
        mm = '0' + strtrim(month - 1,1)
     endif
     if (month eq 3) then begin 
        if ((year eq 2008) or (year eq 2012) or (year eq 2016) or (year eq 2020)) then begin
           dd = 29 
           mm = '0' + strtrim(month - 1, 1)
        endif else begin
           dd = 28 
           mm = '0' + strtrim(month - 1, 1) 
        endelse
     endif 
     if (month eq 1) then begin
        dd    = 31
        mm    = 12
        yearx = year - 1
     endif 
  endif

; identify the last FITS file obtained in the prior calendar day
  datadir = strtrim(rootdir,1) + strtrim(yearx,1) + '/' + strtrim(mm,1) + '/' + strtrim(dd,1) + '/*/XRT*fits*'
  file    = file_search(datadir,count=cnt) ;


  if (cnt eq 0) then begin 
     datadir = strtrim(rootdir,1) + strtrim(yearx,1) + '/' + strtrim(mm,1) + '/*/*/XRT*fits*'
     file    = file_search(datadir,count=cnt) ;
     print,'XRT_ECID: no prior day files...switching to contingency. Be alerted.'  
  endif

  if (cnt gt 0) then begin
     NULL = mrdfits(file[cnt-1],0,hdr,/silent)
  endif else begin
     print, 'XRT_ECID:  no valid file found in the prior day....aborting...'
     return
  endelse

; record the last EC_ID number for later use
  last_ecid = sxpar(hdr,'EC_ID')  ;
  last_hdr  = hdr

;---- Search files obtained on the specific date ----- 
;
; [it assumes the same directory structure for data storage as in ISAS.] 
;
  mm = month
  dd = day

  if (month lt 10) then mm= '0' + strtrim(month,1)
  if (day lt 10)   then dd= '0' + strtrim(day,1)

; search for the data files taken in the specified date
  datadir = strtrim(rootdir,1) + strtrim(year,1) + '/' + strtrim(mm,1) + '/' + strtrim(dd,1) + '/*/XRT*fits*'
  file    = file_search(datadir,count=cnt) 


; if no file found (cnt == 0), abort the process.
  if (cnt eq 0) then begin
     print, 'XRT_ECID:  no valid file found for the given date...aborting.'
     return
  endif

; initialize index arrays. 
  n_fi   = n_elements(file)
  idx    = lonarr(n_fi)
  zidx   = fltarr(n_fi)
  tstamp = strarr(n_fi)
  tzidx  = fltarr(n_fi)


; read data files in, record its EC_ID and if its data array contains
; some missing data packets (pixel_value == 0). 
;
; for the packet size, see Chapter 16, page 38 in Solar-B blue
; book. For simplicity, I set the minimum packet unit to be 64
; pixels. 
;
  for i=0, n_fi-1 do begin
     NULL      = mrdfits(file[i],0,hdr,/silent)
     idx[i]    = sxpar(hdr,'EC_ID')
     tstamp[i] = sxpar(hdr,'DATE_OBS')
     hz        = where(NULL eq 0,zp)
     zp        = zp/64.         ;1 packet = 64Kpixel = 64 * 1024 pixels = 65536 pixels max. 
     zidx[i]   = zp
     tzidx[i]  = n_elements(NULL)/64. ; again, 1 packet
  endfor

  ;; detect if images on previous day taken with low level command,
  ;; added KKR 12-15-09
  if (last_ecid lt 0) then begin
     print,  'negative EC_ID detected for previous day: Some EC_ID tag values are corrupted or '
     print,'                     images taken with a low level command (MDP_XRT_EXP_START).'
     last_ecid = idx[0] - 1
  endif

;; logic added to detect XRT SAFEHOLD mode. 
;; changed to check if idx[0] eq 1.  Example on 29 May 2010 has EC_ID
;; reset to 1, not 0

  if ((last_ecid ne 32767) and (idx[0] eq 1)) then begin
     last_ecid = idx[0] - 1     ;
     print,'XRT_ECID:  oddity identified...probably XRT in SAFEHOLD lately?';
     print,'XRT_ECID:  (last_ecid reset to the latest index value)'         ;
  endif


; EC_ID rolls over from 32767 to 0. Check if the index has rolled over
; and then if it did, add 32768 to make it look like a running index. 
; This logic fails if there are more than 32768 files obtained in one
; day. Very unlikely considering one must obtain one image every
; 2.64s on average. I wouldn't say it's inconceivable, though. 
;
; Need to allow for missing images around rollover.  If difference
; between max(idx) and min(idx) is large, probably a rollover occured.
; Will fail if there is a continuous block of 767 missing images that
; bridges the time of the wrap.  Unlikely, but not impossible.
; Added 2010-12-08 KKR

  if (max(idx) - min(idx) gt 32000) then begin
     hh          = where(idx eq max(idx)) 
     idx[hh+1:*] = idx[hh+1:*] + 32768
     print,'EC_ID index rolling over...'
  endif 

; when a low-level command is used to take an image (e.g.,
; MDP_XRT_EXP_START), a negative EC_ID is registered. It
; doesn't seem to matter if the command is issued in OP or RT.
; In any case, those EC_IDs are flagged as "hbad" with the value 
; of -9999. 
  hbad = -1
  if (min(idx) lt 0) then begin
     hbad = where(idx lt 0, mbad) 
     print, strtrim(mbad,1) + ' negative ID detected: Some EC_ID tag values are corrupted or '
     print,'                     images taken with a low level command (MDP_XRT_EXP_START).'
     for j=0,mbad-1 do begin
        idx[hbad[j]] = idx[hbad[j]-1]
     endfor
  endif


; Derive how many gaps in EC_ID are recorded. Use last_ecid here. 
; 
; The idea here is "simple". Suppose idx array looks like this:
;
; ... 12000 12001 12002 12004 12005 12006 12006 -32767 12007 ...
; 
; then the difference "idx - shift(idx,1) is 
; 
; ... 12000 12001 12002 12004 12005 12006 12006 -32767 12007 ...
; minus
; ... 11999 12000 12001 12002 12004 12005 12006  12006 -32767 ...
; ---------------------------------------------------------------
; ...     1     1     1     2     1     1     0 -44773 44774 ...
; 
; Now there are four cases here to consider.
;
; (a) There is no gap, then the difference is 1. 
; (b) There is a gap, then the difference is greater than 1. 
; (c) Identical EC_IDs in sequence, then the difference is 0. 
; (d) There is a negative EC_ID, then the value is gibberish. 
;
; To identify missing files correctly, the cases (c) and (d) need to
; be dealt with before summing the missing file numbers. 
;
; For (c), we flag wherever the difference is 0 and then replace
; the value with 1 (since it is not actually missing). Later we 
; recall the flag and set the index value (idx) to be -1 for
; these. 
;
; For (d) (which we already dealt with in the "IF" loop above), we
; identify the negative EC_ID index prior to taking the difference,
; i.e., changing idx from 
; 
; ... 12000 12001 12002 12004 12005 12006 12006 -32767 12007 ...
; 
; to 
;
; ... 12000 12001 12002 12004 12005 12006 12006  12006 12007 ...
;
; then subtract shift(idx,1) 
; 
; ... 11999 12000 12001 12002 12004 12005 12006  12006 12006 ...
; ---------------------------------------------------------------
; ...     1     1     1     2     1     1     0      0     1 ...
; 
; So the difference of those negative EC_ID indices would be set to 0
; as well. Since we flag these negative ids as hbad, we can
; distinguish between the negative IDs and identical IDs. Those
; negative IDs are set to -9999 later on. 
;
; Phew, I went on a bit of hyperbole here, eh? 
  idx2    = idx - shift(idx,1)
  idx2[0] = idx2[0] + idx[n_fi-1] - last_ecid


; Flagging the identical EC_IDs (idx2 == 0) and then count how many
; gaps there are. If idx2 == 0, then there is no file actually missing
; and therefore we need to re-initialize these to the value of 1 so
; that they won't be counted as missing. 
;
; The variable idx should generally be zero when there's no gap. 

  hm = where(idx2 eq 0,mcnt)
  if (mcnt ne 0) then begin
     print,''
     print,'XRT_ECID: Two identical EC_ID detected (flagged as -1).'
     print,''
     idx2[hm] = 1
     miss     = fix(total(idx2-1))
     idx2[hm] = 0
     idx      = idx2-1
  endif else begin
     miss     = fix(total(idx2-1))
     idx      = idx2-1
  endelse

; Flag the negative EC_IDs as -9999
  if (hbad[0] ne -1) then idx[hbad] = -9999


; Evaluate some of the identified gaps to see if the files really are
; missing. 
;
; [Here is where true nightmare begins.]
;
; It turns out that EC_IDs aren't gap-less all the time. Here
; are some cases where there is a gap in EC_IDs. 
;
; (a) Case I: Switching XOBs. 
;     If CTRL_MANU interrupts an exposure in order to switch to a new
;     XOB, then the exposure would be discarded and there will be a
;     gap in EC_ID. It is also possible that the last exposured image
;     is truly lost (hence a gap in EC_ID), but we may not have a
;     trivial way to identify that case (*). Anyway, in this case, we
;     consider that the file is NOT missing. 
; 
; (b) Case II: Stop, Pause, (re-point?), and Set the Same XOB.
;     The same XOB will be set twice consecutively into QT slot (e.g.,
;     run an XOB, CTRL_MANU, re-point, and then set the same XOB to
;     run again). When this event takes place, the MAIN, SUBR, and
;     SEQN counters are reset.  Rare, but this does happen. 
;
; (c) Case III: Stop and Start (no reset in QT, using the same XOB)
;     XOB is halted when entering SAA and NGT; likewise, XOB is resumed
;     when exiting from SAA and NGT. CTRL_MANU and CTRL_AUTO are
;     issued without setting a new XOB into QT. This event will
;     not reset the counters; instead it will increment the value of
;     some counters. 
;
; There may be other cases where the gap in EC_ID is triggered.
; 
; Handling (a) and (b) are not that hard. All I had to do is to check
; the PROG_VER value as well as some of the MAIN, SUBR, and SEQN
; counter values to identify if the (a) or (b) casese apply. 
; 
; However, the case (c) isn't trivial, as one can imagine. To
; handle this one, one must import the obev file that covers the
; particular date of user's choice. Or else, the algorithm will
; not identify this condition and some files may remain to be
; incorrectly identified as missing. 
;
; Unfortunately, this part of the condition evaluation algorithm is
; not fail-proof. Users, be aware. 
  hh     = where(idx gt 0)
  nomiss = 0

  if (hh[0] ne '-1' and hh[0] ne 0) then begin 
     n_fih = n_elements(hh)

     for i=0, n_fih-1 do begin 
        NULL = mrdfits(file[hh[i]],0,hdr,/silent)
        
        if (idx[hh[i]] le 2) then begin
;
; if two ore more files are missing, then I suspect them to be
; completely missing. 
;
           prgv1 = sxpar(hdr,'PROG_VER')
           manr1 = sxpar(hdr,'MAIN_RPT')
           sbrn1 = sxpar(hdr,'SUBR_NO')
           sbrc1 = sxpar(hdr,'SUBR_CNT')
           sbrr1 = sxpar(hdr,'SUBR_RPT')
           sbrp1 = sxpar(hdr,'SUBR_POS')
           if (hh[i] eq 0) then begin
              prgv0 = sxpar(last_hdr,'PROG_VER')
           endif else begin
              NULL  = mrdfits(file[hh[i]-1],0,hdr,/silent)
              prgv0 = sxpar(hdr,'PROG_VER')
           endelse

; Case I: Switching XOBs. 
           if (prgv0 ne prgv1) then begin
              idx[hh[i]] = idx[hh[i]] - 1
              nomiss     = nomiss + 1 
           endif 

; Case II: Stop, Pause, (re-point?), and Set the Same XOB.
           if ((prgv0 eq prgv1) and ((manr1 eq 0) and (sbrn1 eq 1) and (sbrc1 eq 1) and (sbrr1 eq 1) and (sbrp1 eq 1))) then begin
              idx[hh[i]] = idx[hh[i]] - 1
              nomiss     = nomiss + 1 
           endif


; Case III: Stop and Start.
           if ((prgv0 eq prgv1) and ((manr1 ne 0) or (sbrn1 ne 1) or (sbrc1 ne 1) or (sbrr1 ne 1) or (sbrp1 ne 1))) then begin
              if keyword_set(obevfile) then begin
                 dtime = 65.536 + 15. ;the longest exposure plus an adhoc value. 
;
; If the gap occurred near the CTRL_MANU time (due to SAA or
; NGT_ENTRY), then it is likely that the file is not actually
; missing. See hcm below. 
;
                 hcm   = where((obev.tsec gt anytim(tstamp[hh[i]-1])) and (obev.tsec lt anytim(tstamp[hh[i]-1])+dtime) and (obev.tsec lt anytim(tstamp[hh[i]])), hcmcnt)
                 if (hcmcnt gt 0) then begin
                    idx[hh[i]] = idx[hh[i]] - 1
                    nomiss     = nomiss + 1
                 endif
              endif else begin
                 print,'XRT_ECID: There may be some files that are not actually missing.'
                 print,'        :  Use obevfile option if you want to avoid them.'
              endelse
           endif

        endif
     endfor   
  endif


; Now it's done. 
  miss = miss - nomiss
  miss_pack = total(zidx)
  tot_pack  = total(tzidx)
 
  print, 'miss ', miss
  print, 'nomiss', nomiss

  return
end


; 2007,1,8: some corrupted files encountered. Missing pixels is not
; multiple of 64 here.
;
; 2007,1,19: two identical ec_id detected. 
; 2007,1,20: two identical ec_id detected. 
; 2007,1,27: two identical ec_id detected. 
; 2007,2,3:  two identical ec_id detected. 
; 2007,2,5:  something really wicked happening here. 
; 2007,3,29: ditto
;2008,5,17: all 2 missing files.