;+
;NAME:
;	hxtimg_accum
;PURPOSE:
; 	plots hxt data, prompts for start and end times for data 
;	accumulation and background, and returns an array with the
;       accumulated counts, for each interval and channel.
;CATEGORY:
;       Yohkoh HXT
;CALLING SEQUENCE:
;       Hxtimg_accum, index, data, index1, hxi_index, counts, sigma, bcounts, $
;                     bsigma, img_trange = img_trange, bck_trange=bck_trange, $
;                     dt_resolution=dt_resolution, alt_tsel=alt_tsel, $
;                     accum_cnts=accum_cnts, ch_accum=ch_accum, cnts_str=cnts_str,$
;                     bkgd_str=bkgd_str, hdafile=hdafile, channel=channel, $
;                     noback=noback, chisqr=chisqr, st_times = st_times, $
;                     int_times = int_times
;INPUT:	
;	index=data info from hda file
;	data=counts in each subcollimator for each channel
;KEYWORDS: (INPUT)
;       img_trange = a time range for the images
;       bck_range = a time_range for the background accumulation
;       dt_resolution = number of 1/2 second (flare mode) or 4 second
;                      (quiet mode) frames per image, default is to
;                      prompt the user. If set to -1, then the whole chosen
;                      interval is used. NOTE THAT THIS KEYWORD IS IGNORED
;                      IF MULTIPLE INTERVALS HAVE BEEN SELECTED EITHER
;                      IN PLOT_LCUR, OR IF THEY HAVE BEEN PASSED IN.
;       alt_tsel = use alt_tim_sel to obtain interval, rather than plot_lcur
;                  Useful if you don't have x-windows
;       accum_cnts = if set, the value of average counts/SC in a given
;                    that you want for each set of images, the routine
;                    will accumulate counts until the number of counts
;                    in the given channel are greater than or equal to
;                    this value. If set with a slash, ie. the value=1
;                    then the default is 200 (a good value).
;       accum_cnts_ch = the channel you want to use for the accum_cnts
;       structures = if set, the data will be passed out in structures,
;                    duplicating the use of MK_HXT_SUM for use in HXTPRO.
;                    The structure hxi_index need not be passed in if
;                    the structures keyword is set. THIS IS IMPORTANT:
;                    DO NOT USE ANY OF THE OPUTPUT VARIABLES IF YOU SET
;                    THIS KEYWORD, I.E. THE CALLING SEQUENCE ONLY INCLUDES
;                    INDEX AND DATA AS PARAMETERS, THE REST MUST BE KEYWORDS,
;                    LIKE THIS: Hxtimg_accum, index, data, /struct, cnts_str=cnts_str, $
;                                  bkgd_str=bkgd_str, other keywords....
;       hdafile= the name of the hda file used, the structures want this
;                i have no clue why, if not set we make one up.
;       channel= a channel number for the plots, default is to prompt the user
;       noback= if set, the background is set to zero. Requested by HSH 21-jun-95
;       /chisqr= if set, sigma(counts) is a weighted sum of total(sqrt(counts))
;                 and sqrt(total(counts)), where the total is taken over the
;                 accumulation intervals. This may be more useful for cases
;                 with low signal-to-noise. This used to be the default,
;                 but not anymore, 6-mar-96, the default is to use sqrt(total(counts))
;       st_times= an array of accumulation start times, in any yohkoh format
;       int_times= interval times corresponding to st_times, may be a scalar,
;                  if this is not set, and st_times is set, the routine will
;                  use the differences between st_times. This must be an integer
;                  number of the appropriate resolution in seconds, if the
;                  value of int_times is less than a frame value.
;KEYWORDS: (OUTPUT)
;       cnts_str = the counts in a structure, the structures keyword must
;                  be set for this to be calculated, for use by HXTPRO
;       bkgd_str = the background in a structure, the structures keyword must
;                  be set for this to be calculated.
;OUTPUT:
;       index1= index structure for the first frame of each image interval
;	hxi_index.chan= channel used
;		=0 L0 (13-23 keV)
;		=1 M1 (23-33 keV)
;		=2 M2 (33-53 keV)
;		=3 HI (53-93 keV)
;	hxi_index.actim= image accumulation time in 10ths of secs
;	hxi_index.interval= actim
;	hxi_index.bkg_time= start time for background accumulation
;                           (millisec OF day)
;	hxi_index.bkg_day= day of start of background accumulation
;                          (days since 1-Jan-79)
;	hxi_index.bkg_actim= background accumulation time in 10ths of secs
;	hxi_index.bkg_interval= actim 
;       counts= the counts accumulated over the chosen intervals
;       sigma= sqrt(counts)
;       bcounts= the counts accumulated over the background interval
;       bsigma= sqrt(bcounts)
;DESCRIPTION:
;	Uses the function plot_lcur to select the data to be imaged,
;	and the background.
;MODIFICATION HISTORY:
;	Written Oct-93 by JMM,
;       Changed to give full resolution, 19-jul-94, JMM
;       9/21/94, account for 4.0 second time lag in hda files, jmm
;       1/9/95, changed calculation of sigma to be the same as
;               in Tom Metcalf's hxtaccumulate.pro, to be consistent
;               in pixon calculations.
;       Added the ability to do non-flare mode data, 2/12/95, jmm
;       Speed things up using new-style TOTAL commands, 2/13/95, jmm
;       13-mar-95 JPW : img_trange can now be an array of ranges (2,n)
;       16-mar-95 JPW : can now use multiset option in plot_lcur
;       Stripped out accumulation stuff, replaced it with calls to
;       CHOOSE_INTERVAL, 5-apr-95, jmm
;       added cnts_str, bkgd_str keyword outputs, 6-apr-95, jmm
;       Added the call to setgra, which gives the subscripts for
;       the subcollimator switch, from HXTSBT_MULT, 6-apr-95, jmm.
;       This accomodates the structures...
;       I don't think anyone will notice.
;       Subtract bcakground before choosing the imaging intervals,
;       You can see where it's less than zero, and ACCUM_CNTS now
;       refers to Excess counts, 9-jun-95
;       Added /noback option, 22-jun-95
;       Fixed bug in calculation of SC uncertainties for single
;       interval accumulations, 1-Aug-95, Jmm
;       Changed  default sigmas back to poisson statistics, added 
;       /chisqr keyword, 6-mar-96, jmm
;       Added st_times, and int_times keywords, 3-jan-1997, jmm
;       Call to hxt_chsums is changed, 6-apr-1999, jmm
;-
PRO  Hxtimg_accum, index0, data, index1, hxi_index, counts, sigma, bcounts, $
                   bsigma, Img_trange=img_trange, bck_trange=bck_trange, $
                   dt_resolution=dt_resolution, alt_tsel=alt_tsel, $
                   accum_cnts=accum_cnts, ch_accum=ch_accum, structures=structures, $
                   Cnts_str=cnts_str, Bkgd_str=bkgd_str, Hdafile=hdafile, $
                   channel=channel, noback=noback, chisqr=chisqr, $
                   st_times=st_times, int_times=int_times

;The structure hxi_index may not be passed in, if the structures keyword is set
   IF(N_ELEMENTS(hxi_index) GT 0) THEN qhxi_index = 1 ELSE qhxi_index = 0
   IF(qhxi_index) THEN serr = hxi_index(0).syserr ELSE serr = 0.0
   IF(KEYWORD_SET(chisqr)) THEN qchisqr = 1 ELSE qchisqr = 0
   
   IF(KEYWORD_SET(structures) AND n_params() GT 2) THEN BEGIN
      message, /info, 'THE STRUCTURES KEYWORD HAS BEEN SET'
      message, /Info, 'DO NOT USE ANY NON-KEYWORDS IN THE FUNCTION CALL OTHER THAN INDEX AND DATA'
      RETURN
   ENDIF
   
;sum SC's to plot and you'll need a time array for each data point
;Also deal with any Cal mode data.
   zz = hxt_chsums(index0, data, /half_sec, tim_arr = tim_arr1, $
                   tim_sec = tim_arr, dt = dx4, index_out = index, $
                   data_out = data_x)
   time0 = tim_arr1(0)
   nc0 = 4                      ;no. of channels
   nc01 = nc0-1
   ndset = N_ELEMENTS(index)	;no. of data sets
   nds1 = ndset-1
;For full resolution, 4 minor frames for each data set
   nds4 = ndset*4 & nds41 = nds4-1
;reform data to get into the proper form
   data_x = reform(data_x, nc0, 64, nds4, /overwrite)

;If i need to choose intervals by plotting
   no_time_selected = NOT (KEYWORD_SET(img_trange) OR KEYWORD_SET(st_times))
   no_bck_selected = NOT (KEYWORD_SET(bck_trange) OR KEYWORD_SET(noback))
   IF(no_time_selected OR no_bck_selected) THEN BEGIN
      IF(N_ELEMENTS(channel) GT 0) THEN ich = fix(channel(0)) ELSE BEGIN
         print, ' INPUT CHANNEL NUMBER FOR THE INTERVAL>>'
         print, ' 0 FOR LO CHANNEL'
         print, ' 1 FOR M1 CHANNEL'
         print, ' 2 FOR M2 CHANNEL'
         print, ' 3 FOR HI CHANNEL'
         read, ich
         ich = fix(ich)
      ENDELSE
   ENDIF ELSE ich = 0
   y0 = reform(zz(ich, *))/(dx4*64.0) ;y0 is the curve used to find the intervals
   IF(KEYWORD_SET(noback)) THEN BEGIN
      print, 'BACKGROUND SET TO ZERO, NOT RECOMMENDED'
      bcounts = fltarr(64, 4)
      bsigma = bcounts
      ib1 = 0 & ib2 = 3         ; i may need this stuff, for the structures
      ib14 = ib1/4 & ib24 = ib2/4
      bactim = total(dx4(ib1:ib2))
   ENDIF ELSE BEGIN      
;choose the background interval
      print, 'BACKGROUND SELECTION::'
      choose_interval, tim_arr1, ib1, ib2, data = y0, plot_title = 'HXT COUNTS/SECOND/SC', $
        time_range = bck_trange, dt_resolution = -1, alt_tsel = alt_tsel
      ib1 = ib1(0) & ib2 = ib2(0)

; Check time interval for non-flare mode data
      ib14 = ib1/4 & ib24 = ib2/4
      not_flare = where(gt_dp_mode(index(ib14:ib24)) NE 9)
      IF(not_flare(0) NE -1) THEN BEGIN
         message, /info, 'WARNING, NON-FLARE MODE BACKGROUND CHOSEN, IS ZERO FOR HI, M1, AND M2 CHANNELS'
      ENDIF
   
; time interval,
      bactim = total(dx4(ib1:ib2))
; Counts array
      tcounts = data_x(*, *, ib1:ib2)
;Speed things up, using total commands on 1 dimension
      nibs = ib2-ib1+1
      IF(nibs GT 1) THEN BEGIN
; new stuff, see explanation below
; bsigma = sqrt(bcounts) is not true anymore
; error1 is total(sqrt(counts)), bsigma is sqrt(total(counts))
         bcounts = transpose(total(tcounts, 3))
         bsigma = transpose(sqrt(total(tcounts, 3)+ $
                                 total((serr*tcounts)^2, 3)))
         IF(qchisqr) THEN BEGIN ; jmm, 6-mar-96
            error1 = transpose(sqrt(total(sqrt(tcounts), 3)^2 $
                                    +total((serr*tcounts)^2, 3)))
            mn = avsig_1(tcounts, sdev, dim = 3) ;mean and st deviation
            mn = mn > 0
            xxx = where(sdev GT 0.0)
            mn(xxx) = ((1.0-sqrt(mn(xxx))/sdev(xxx)) > 0.0) < 1.0
            wt = transpose(mn)
            bsigma = error1*wt + bsigma*(1.0-wt)
         ENDIF ELSE bsigma = bsigma
      ENDIF ELSE BEGIN          ;no need to total
         bcounts = transpose(tcounts)
         bsigma = transpose(sqrt(tcounts+(serr*tcounts)^2)) ;oops, jmm,1-aug-95
      ENDELSE
   ENDELSE

; Ok, now you choose the image interval, subtract the background first
; you need a sum for the background
   IF(KEYWORD_SET(noback)) THEN bcps = fltarr(4) ELSE BEGIN
      IF(ib2 NE ib1) THEN bcps = total(zz(*, ib1:ib2), 2) ELSE bcps = zz(*, ib1)
   ENDELSE
   bcps = bcps/(64.0*bactim)    ;average background rate per SC

   print, 'IMAGE INTERVAL SELECTION::'
   IF(KEYWORD_SET(accum_cnts)) THEN BEGIN
      IF(accum_cnts EQ 1) THEN acc_cnts = 200.0 ELSE acc_cnts = accum_cnts(0)
      IF(N_ELEMENTS(ch_accum) GT 0) THEN ich0 = ch_accum(0) ELSE ich0 = 0
      y0 = reform(zz(ich0, *)/64.0)
      y0 = y0-bcps(ich0)*dx4
      choose_interval, tim_arr1, ss_st, ss_en, data = y0, accum_cnts = acc_cnts, $
        time_range = img_trange, plot_title = 'BACKGROUND SUBTRACTED HXT COUNTS/SC', $
        st_times = st_times, int_times = int_times
   ENDIF ELSE BEGIN
      y0 = y0-bcps(ich)
      choose_interval, tim_arr1, ss_st, ss_en, data = y0,  $
        plot_title = 'BACKGROUND SUBTRACTED HXT COUNTS/SEC/SC', $
        time_range = img_trange, dt_res = dt_resolution, alt_tsel = alt_tsel, $
        st_times = st_times, int_times = int_times
   ENDELSE
; Check time interval for non-flare mode data
   ix14 = min(ss_st)/4 & ix24 = max(ss_en)/4
   not_flare = where(gt_dp_mode(index(ix14:ix24)) NE 9)
   IF(not_flare(0) NE -1) THEN BEGIN
      message, /info, 'WARNING, NON-FLARE MODE DATA IN CHOSEN INTERVAL IS ZERO FOR HI, M1, AND M2 CHANNELS'
   ENDIF

   n_intv = N_ELEMENTS(ss_st)
;You'll not have an index structure for each minor frame, so for
;each minor frame, find the appropriate index, and reset .time
;and .day
   ss_st_index = ss_st/4 & index1 = index(ss_st_index)
   index1.gen.time = tim_arr1(ss_st).time
   index1.gen.day = tim_arr1(ss_en).day

; now accumulate data, dti will be the time interval length, in seconds
   counts = fltarr(64, nc0, n_intv)
   IF(qchisqr) THEN BEGIN ;newer stuff, jmm, 6-mar-96
      wt = counts               ;new stuff, 1/9/95, jmm
      error1 = counts
   ENDIF
   sigma = counts

   dti = fltarr(n_intv)
   FOR j = 0, n_intv-1 DO BEGIN
      dti(j) = total(dx4(ss_st(j):ss_en(j)))
      tcounts = data_x(*, *, ss_st(j):ss_en(j))
      dnn = ss_en(j)-ss_st(j)+1
      IF(dnn GT 1) THEN BEGIN
         counts(*, *, j) = transpose(total(tcounts, 3))
         sigma(*, *, j) = transpose(sqrt(total(tcounts, 3) $
                                          +total((serr*tcounts)^2, 3)))
         IF(qchisqr) THEN BEGIN ;newer stuff, jmm, 6-mar-96
            error1(*, *, j) = transpose(sqrt(total(sqrt(tcounts), 3)^2 $
                                             +total((serr*tcounts)^2, 3)))
            mn = avsig_1(tcounts, sdev, dim = 3) ;mean and st dev
            mn = mn > 0
            xxx = where(sdev GT 0.0)
            mn(xxx) = ((1.0-sqrt(mn(xxx))/sdev(xxx)) > 0.0) < 1.0
            wt(*, *, j) = transpose(mn)
         ENDIF ELSE BEGIN
            error1 = 0.0
            wt = 0.0
         ENDELSE
      ENDIF ELSE BEGIN          ;no need to total
         counts(*, *, j) = transpose(tcounts)
         sigma(*, *, j) = transpose(sqrt(tcounts+(serr*tcounts)^2))
         error1 = 0.0
         wt = 0.0
      ENDELSE
   ENDFOR

   ; now calculate the uncertainty=sqrt(counts)

   ; The error would be total(sqrt(counts)) if each observation
   ; were independent.  The error would be sqrt(total(counts))
   ; if each observation were the same.  Reality is in between.
   ; Use the standard deviation compared to the error on a single
   ; observation to weight each type of error.

   ; error1 is total(sqrt(counts)), sigma is sqrt(total(counts))

   sigma = error1*wt + sigma*(1.0-wt)  
   
; Replicate and fill hxi_index if necessary,
   IF(qhxi_index) THEN BEGIN
      hxi_index.bkg_time = tim_arr1(ib1).time
      hxi_index.bkg_day = gt_day(tim_arr1(ib1))
      hxi_index.bkg_actim = fix(10.0*bactim)
      hxi_index.bkg_interval = fix(10.0*(tim_arr(ib2)-tim_arr(ib1)+dx4(ib2)))      
      hxi_index = replicate(hxi_index(0), nc0, n_intv)
      FOR i = 0, nc01 DO hxi_index(i, *).chan = i
      FOR j = 0, n_intv-1 DO BEGIN
         hxi_index(*, j).actim = fix(10.0*dti(j))
         hxi_index(*, j).interval = fix(10.0*(tim_arr(ss_en(j))-tim_arr(ss_st(j))+dx4(ss_en(j))))
      ENDFOR
   ENDIF
   data_x = 0.0

; hxt_setgra gives the subscripts for the subcollimator switch
   hxt_setgra, sn2gra, gra2sn
   bcounts = bcounts(gra2sn, *)
   counts = counts(gra2sn, *, *)
   bsigma = bsigma(gra2sn, *)
   sigma = sigma(gra2sn, *, *)

;For structure, this should be moderately easy
   IF(KEYWORD_SET(structures)) THEN BEGIN
      ss_en_index = ss_en/4
      set_hxt_struct, cnts_str  ;set counts structure
      cnts_str.gen.id = ' HXTIMG_ACCUM Output' ;and start filling
      cnts_str.inf_raw.comment = 'No Comment'
      IF(KEYWORD_SET(hdafile)) THEN cnts_str.inf_raw.rf_file = hdafile $
        ELSE cnts_str.inf_raw.rf_file = 'hda'+ex2fid(anytim2ex(index(0)))
      
      cnts_str = replicate(cnts_str(0), n_intv)
;the times look formatted. This looks like a mess....
      t7 = anytim2ex(tim_arr1(ss_st))
      t7_str = strcompress(string(t7), /remove_all)
      lt_10 = where(t7 LT 10)
      IF(lt_10(0) NE -1) THEN t7_str(lt_10) = ' '+t7_str(lt_10)
      lt_100 = where(t7(3,*) LT 100)
      IF(lt_100(0) NE -1) THEN t7_str(3, lt_100) = ' '+t7_str(3, lt_100)
      t71 = anytim2ex(anytim2ints(tim_arr1(ss_st), offset=dti))
      t71_str = strcompress(string(t71), /remove_all)
      lt_10 = where(t71 LT 10)
      IF(lt_10(0) NE -1) THEN t71_str(lt_10) = ' '+t71_str(lt_10)
      lt_100 = where(t71(3,*) LT 100)
      IF(lt_100(0) NE -1) THEN t71_str(3, lt_100) = ' '+t71_str(3, lt_100)
      time_str = t7_str(6, *)+'/'+t7_str(5, *)+'/'+t7_str(4, *) $
        +' '+t7_str(0, *)+':'+t7_str(1, *)+':'+t7_str(2, *)+'.'+t7_str(3, *) $
        +' - '+t71_str(0, *)+':'+t71_str(1, *)+':'+t71_str(2, *)+'.'+t71_str(3, *)
      FOR j = 0, n_intv-1 DO BEGIN
         cnts_str(j).inf_raw.range(0, *) = ss_st_index(j) ;the subscript in index for the start
         cnts_str(j).inf_raw.range(1, *) = ss_st(j)-4*ss_st_index(j) ;start minor frame
         cnts_str(j).inf_raw.range(2, *) = ss_en_index(j) ;the subscript in index for the end
         cnts_str(j).inf_raw.range(3, *) = ss_en(j)-4*ss_en_index(j) ;end minor frame
         cnts_str(j).inf_raw.time(*) = time_str(j)
         cnts_str(j).data.acc_tim(*) = dti(j)
         cnts_str(j).data.cts = counts(*, *, j)
         cnts_str(j).data.sigma = sqrt(counts(*, *, j))
      ENDFOR

      set_hxt_struct, bkgd_str  ;set counts structure
      bkgd_str.gen.id = ' HXTIMG_ACCUM Output' ;and start filling
      bkgd_str.inf_raw.comment = 'No Comment'
      IF(KEYWORD_SET(hdafile)) THEN bkgd_str.inf_raw.rf_file = hdafile $
        ELSE bkgd_str.inf_raw.rf_file = 'hda'+ex2fid(anytim2ex(index(0)))

;the times look formatted.
      t7 = anytim2ex(tim_arr1(ib1))
      t7_str = strcompress(string(t7), /remove_all)
      lt_10 = where(t7 LT 10)
      IF(lt_10(0) NE -1) THEN t7_str(lt_10) = ' '+t7_str(lt_10)
      IF(t7(3) LT 100) THEN t7_str(3) = ' '+t7_str(3)
      t71 = anytim2ex(anytim2ints(tim_arr1(ib1), offset = bactim))
      t71_str = strcompress(string(t71), /remove_all)
      lt_10 = where(t71 LT 10)
      IF(lt_10(0) NE -1) THEN t71_str(lt_10) = ' '+t71_str(lt_10)
      IF(t71(3) LT 100) THEN t71_str(3) = ' '+t71_str(3)
      bkgd_str.inf_raw.time(*) = t7_str(6)+'/'+t7_str(5)+'/'+t7_str(4) $
        +' '+t7_str(0)+':'+t7_str(1)+':'+t7_str(2)+'.'+t7_str(3) $
        +' - '+t71_str(0)+':'+t71_str(1)+':'+t71_str(2)+'.'+t71_str(3)
      bkgd_str.inf_raw.range(0, *) = ib14 ;the subscript in index for the start
      bkgd_str.inf_raw.range(1, *) = ib1-4*ib14 ;start minor frame
      bkgd_str.inf_raw.range(2, *) = ib24 ;the subscript in index for the end
      bkgd_str.inf_raw.range(3, *) = ib2-4*ib24 ;end minor frame
      bkgd_str.data.acc_tim(*) = bactim
      bkgd_str.data.cts = bcounts
      bkgd_str.data.sigma = sqrt(bcounts(*, *))
   ENDIF
         
   RETURN
END