;---------------------------------------------------------------------------
;+
; NAME:
;     avsig_1
; PURPOSE:
;     Average and dispersion of an array, zeros can be not included,
; CALLING SEQUENCE:
;     xbar=avsig_1(x, sigma, no_zeros=no_zeros, sig_mean=sig_mean, dimension=dimension, /fractional)
; INPUT:
;     x = an array
; OUTPUT:
;     xbar = mean, total(x)/n_elements(x)
;     sigma = standard deviation, sqrt(total((x-xbar)^2/(nx-1)))
; KEYWORDS:
;     no_zeros= if set, strip out zeros, note that this option does
;               not work if the dimension keyword is used. At least not as of
;               2-13-95...
;     sig_mean = sigma/sqrt(nx), the standard deviation of the mean of the array
;     dimension = the dimension of the array to find the mean in,
;                 passed into the total command, it must be a scalar.
;     fractional = the fractional error is passed out as sigma,
;                  don't use this if zero is a valid value of xbar...
; HISTORY:
;     12-9-94, jmm
;     2-13-95, jmm added dimension keyword, switched from ok_zeros to no_zeros
;     5-sep-1996, jmm, switched to double precision
;-
FUNCTION Avsig_1, x0, sigma, no_zeros=no_zeros, sig_mean=sig_mean, $
                  dimension=dimension, fractional=fractional

   x = double(x0)
   IF(KEYWORD_SET(dimension)) THEN BEGIN
;get the size of the array, if it's only 1d, goto the else part
      size_x = size(x)
      IF(size_x(0) LE 1) THEN GOTO, its_a_vector
;Ok, now be sure that the given dimension exists
      d0 = long(dimension(0))
      IF(d0 GT size_x(0)) THEN BEGIN
         message, 'NOT ENOUGH DIMENSIONS IN ARRAY, USING WHOLE ARRAY', /info
         GOTO, its_a_vector
      ENDIF
;N-elements
      nx = float(size_x(d0))
      xbar = total(x, d0)/nx ;the mean is easy
      sigma = total(x^2, d0) ;gotta use two sums

      IF(nx GT 1.0) THEN sigma = sqrt((sigma-nx*xbar^2)/(nx-1.0)) $
        ELSE sigma = make_array(size = size(xbar)) + 0.0
      sig_mean = sigma/sqrt(nx)
   ENDIF ELSE BEGIN
      its_a_vector: IF(KEYWORD_SET(no_zeros)) THEN BEGIN
         ok = where(x GT 0.0, nx)
         nx = float(nx)
         IF(nx GT 0) THEN BEGIN
            xbar = total(x(ok))/nx
            IF(nx GT 1) THEN sigma = sqrt(total((x(ok)-xbar)^2)/(nx-1.0)) $
              ELSE sigma = 0.0
            sig_mean = sigma/sqrt(nx)
         ENDIF ELSE BEGIN
            xbar = 0.0
            sigma = 0.0
            sig_mean = 0.0
         ENDELSE
      ENDIF ELSE BEGIN
         nx = float(N_ELEMENTS(x))
         xbar = total(x)/nx
         IF(nx GT 1) THEN sigma = sqrt(total((x-xbar)^2)/(nx-1.0)) $
           ELSE sigma = 0.0
         sig_mean = sigma/sqrt(nx)
      ENDELSE
   ENDELSE

   IF(KEYWORD_SET(fractional)) THEN BEGIN
      xxx = where(xbar NE 0.0)
      fraction = xbar & fraction(*) = 0.0
      fraction(xxx) = sigma(xxx)/xbar(xxx)
      sigma = fraction
   ENDIF

   RETURN, xbar
END

;+
;NAME:
; hsi_obs_summ_allrates
;PROJECT:
; HESSI
;CATEGORY:
; Hessi catalog generation
;PURPOSE:
; Calculates count rates from hessi eventlist data, to be
; used in the HESSI Observing Summary. The idea is that the
; eventlist should only have to be created once per file
;CALLING SEQUENCE:
; hsi_obs_summ_allrates, infile, time_intv = time_intv, $
;                        energy_edges = energy_edges, $
;                        Quiet=quiet, $
;                        simulated_data=simulated_data, $
;                        obs_time_interval=obs_time_interval
;INPUT:
; infile = A fits file to get data from
; data = a data structure of type  {Hsi_obssummratedata}
;OUTPUT:
; All in a common block hsi_obs_summ_rates,
; pkt = the packet array for the given eventlist
; pkttbl = the packet array lookup table
; pktt = the packet collect time in anytim format, sec from 1-jan-1979
; rates = the count rate, for each detector for each energy band
; rates_mvar = the modulation variance front detectors, colls 7 and 8
; ut_ref_rates = reference time for the rates
; dt_rates = the interval time
; ntimes_rates = the number of time intervals
; av_lvtim = av live time, for all detectors
;KEYWORDS:
; time_intv = the nominal duration, in seconds,
;             of a single time interval, note that
;             there is no info about the duration of
;             individual time intervals, the default is 4.0 seconds
; energy_edges = the energy bin edges
; quiet = if set, run quietly, this is the default
; obs_time_interval = if passed in, the obs_summary is calculated for this
;                  time_range, and the summary start and end times
;                  will be set to this range
; simulated_data = set to 1 for simulated data
;HISTORY:
; Hacked from hsi_obs_summ_fill, 21-nov-2000, jmm, jimm@ssl.berkeley.edu
; Added particle rate calculation, 30-nov-2000, jmm
; Combined rear segments, 1-dec-2000, jmm
; uses lightcurve object, 19-jan-2001, jmm
; Demands file input adds live time and fastrates, endfileinfo
; keyword, 23-may-2001, jmm
; changed common to an anonymous structure, added _extra,
; 6-jun-2001,jmm
; to avoid problems with lightcurve object, defines ut_ref_rates
; to be the start of the file, or obs_time_interval, not at the 
; start of the eventlist, jmm, 31-oct-2001
; Added det_index and seg_index_mask keywords, total_rates is now per
; detector jmm, 7-nov-2001
;-
PRO Hsi_obs_summ_allrates, infile, $ 
       time_intv = time_intv, $
       eventlist_obj=eventlist_obj, $
       Energy_edges=energy_edges, $
       Quiet=quiet, Xxx=xxx, $
       simulated_data=simulated_data, $
       Var_nbin=var_nbin, $
       Var_erange=var_erange, $
       obs_time_interval=obs_time_interval, $
       endfileinfo=endfileinfo, $
       ref_time=ref_time, $
       file_time_intv=file_time_intv, $
       seg_index_mask=seg_index_mask, $
       det_index=det_index, $
       _extra=_extra

;Hsi_obs_summ_rates holds the count rates
   COMMON hsi_obs_summ_rates, hsi_rates_struct
   
   print, 'SIMULATED_DATA:', KEYWORD_SET(simulated_data)
   IF(KEYWORD_SET(quiet)) THEN notquiet = 0 ELSE notquiet = 1
   t00 = !Stime
   IF(notquiet) THEN print, 'Allrates start:  ', t00
   t00 = anytim(t00)
;Initialize
   pkt = -1 & pktt=-1 & pkttbl=-1 & rates=-1 & tot_rates=-1
   rates_mvar = -1 & part_rates=-1 & av_livtim = -1
   ut_ref_rates = 0 & dt_rates=0 & ntimes_rates=0
;keywords
   IF(keyword_set(det_index)) THEN BEGIN
      segi = bytarr(18)
      segi[[det_index, det_index+9]] = 1
   ENDIF ELSE IF(keyword_set(seg_index_mask)) THEN BEGIN
      segi = seg_index_mask
   ENDIF ELSE segi = bytarr(18)+1
   IF(NOT keyword_set(ref_time)) THEN ref_time = anytim('4-jul-2000 0:00')
   IF(KEYWORD_SET(var_nbin)) THEN vnbin = var_nbin ELSE vnbin = 128
   IF(KEYWORD_SET(var_erange)) THEN ch0v = var_erange $ 
   ELSE ch0v = [6.0, 25.0]
   IF(KEYWORD_SET(energy_edges)) THEN enbins = energy_edges $
   ELSE enbins = [3.0, 6.0, 12.0, 25.0, 50.0, 100.0, 300.0, 1000.0, $
                  2500.0, 20000.0]
   nenbins = N_ELEMENTS(enbins)-1
   rates_ebands = enbins

   IF(KEYWORD_SET(time_intv)) THEN dt_rates = time_intv[0] $
   ELSE dt_rates = 4.0d0
   IF(keyword_set(file_time_intv)) THEN dt_files = file_time_intv $
   ELSE dt_files = 5.0d0*dt_rates
;Can you find the file?
   aa = loc_file(infile[0], loc=loc_aa, count=bb)
   IF(bb EQ 0) THEN BEGIN       ;No file
      message, /info, 'Telemetry File not found:'
      print, infile[0]
      RETURN
   ENDIF
   infile = aa[0]
   
;Get the packet array,
   help, infile
   IF(notquiet) THEN message, /info, 'packet array'
   pkt_obj = obj_new('hsi_packet')
   pkt = pkt_obj -> getdata(/all, filename=infile[0])
   pkttbl = pkt_obj -> getdata(/lookup_table)
   pktt = hsi_sctime2any(pkt.collect_time)
   npkt = N_ELEMENTS(pkt)
   IF(datatype(pkt) NE 'STC') THEN BEGIN
      IF(notquiet) THEN message, /info, 'No Packets'
      RETURN
   ENDIF
      
;Ok, now we have packets, Set up the time range,
   IF(KEYWORD_SET(obs_time_interval)) THEN BEGIN
      print, 'OBS TIME INTERVAL IS SET'
      oti = anytim(obs_time_interval)
   ENDIF ELSE BEGIN
      oti = pkt_obj -> get(/file_time_range)
;All files start an integral number of dt from the start of the day
      oti[0] =  hsi_set_file_time(oti[0], time_intv=dt_files, $
                                  ref_time=ref_time)
   ENDELSE
   dt_total = oti[1]-oti[0]
   tu = hsi_adjust_time_unit(dt_total, 1, quiet=quiet) ;10-nov-2001
   ntimes_rates = long(dt_total/dt_rates)
   IF((dt_total MOD dt_rates) NE 0.0) THEN ntimes_rates = ntimes_rates+1
   oti[1] = oti[0]+ntimes_rates*dt_rates
   ut_ref_rates = oti[0]
;Ok, now you may get the eventlist
   ok100 = where(pkttbl.app_id EQ 100)
   IF(ok100[0] NE -1) THEN BEGIN
      eventlist_obj = obj_new('hsi_eventlist')
      eventlist_obj -> set, source = pkt_obj
      pscore = eventlist_obj -> getdata(/all, time_range=[0, 0], $
                                        simulated_data=simulated_data, $
                                        time_unit=tu)
      t01 = !Stime
      IF(notquiet) THEN BEGIN
         print, 'Allrates has eventlist:  ', t01
         t01 = anytim(t01)
         Elapsed_time = t01-t00
         per_pkt = (elapsed_time/npkt)*1000.0
         print, 'Elapsed time: ', elapsed_time, '  s'
         print, 'Per 1000 Packets: ', per_pkt, '  s/1000pkt'
      ENDIF

;Set up the time range, be careful...
      ut_ref = eventlist_obj -> get(/ut_ref)
      rel_trange = oti-ut_ref   ;can be negative
;Next get the light curve, use all segments
      stop
      oltc1 = obj_new('hsi_lightcurve')
      oltc1 -> set, source = eventlist_obj
      rates1 = oltc1 -> getdata(time_range=rel_trange, $
                                ltc_energy_band=enbins, $
                                this_seg_index=bindgen(18), $
                                seg_index_mask=bytarr(18)+1, $
                                ltc_time_resolution=dt_rates, sum_flag=0, $
                                energy_band=[0, 0], $
                                simulated_data=simulated_data)
;you need to synch the absolute time range up with the oti time range
      ntimes_rates1 = N_ELEMENTS(rates1[*, 0, 0])
      rates = fltarr(ntimes_rates, nenbins, 18)
      tr1 = oltc1 -> get(/absolute_time_range)
      tr1 = ref_time+dt_rates*round((tr1-ref_time)/dt_rates)
      IF(total(abs(oti-tr1)) EQ 0.0 AND ntimes_rates EQ ntimes_rates1) $
         THEN rates = temporary(rates1) ELSE BEGIN
;to do the synching correctly, embed rates1 in a larger array
         txmin = min([oti, tr1], max=txmax)
         nx = long((txmax-txmin)/dt_rates)
         ratesx = fltarr(nx, nenbins, 18)
;x1=1st ss of rates, in ratesx
;y1=1st ss of rates1, in ratesx
         y1 = long((tr1[0]-txmin)/dt_rates)
         ratesx[y1:y1+ntimes_rates1-1, *, *] = rates1
         x1 = long((oti[0]-txmin)/dt_rates)
         rates = ratesx[x1:x1+ntimes_rates-1, *, *]
         delvarx, rates1, ratesx
      ENDELSE
;normalize rates by time
      rates = rates/dt_rates
      
      IF(notquiet) THEN BEGIN
         print, 'ut_ref:', anytim(ut_ref_rates, /ccsds)
         print, 'eventlist ut_ref:', anytim(ut_ref, /ccsds)
      ENDIF
      t01 = !Stime
      IF(notquiet) THEN BEGIN
         print, 'Allrates has rates:  ', t01
         t01 = anytim(t01)
         Elapsed_time = t01-t00
         per_pkt = (elapsed_time/npkt)*1000.0
         print, 'Elapsed time: ', elapsed_time, '  s'
         print, 'Per 1000 Packets: ', per_pkt, '  s/1000pkt'
      ENDIF
      IF(notquiet) THEN message, /info, 'Arrays for mod_variance'
;Do the variance, one number for the lowest three bands,
;for the last two detectors,  default of 6 to 25keV
;For only detectors 8 & 9, front segments
      det = [7, 8]
      rates1 = bytarr(2, ntimes_rates)
;Time resolution??
      dt1 = dt_rates/vnbin
      n1 = ntimes_rates*vnbin
;Another lightcurve
      oltc2 = obj_new('hsi_lightcurve')
      oltc2 -> set, source = eventlist_obj
      ltc2 = oltc2 -> getdata(time_range=rel_trange, ltc_energy_band=ch0v, $
                              ltc_time_resolution=dt1, sum_flag=0, $
                              this_seg_index=det, energy_band=[0, 0], $
                              simulated_data=simulated_data, $
                              seg_index_mask=[bytarr(7), 1, 1, bytarr(9)])
      IF(ltc2[0] NE -1) THEN BEGIN
         nltc2 = N_ELEMENTS(ltc2[*, 0, 0])
         IF(nltc2 GT n1) THEN ltc2 = ltc2[0:n1-1, *, *] ELSE BEGIN
            temp = ltc2
            ltc2 = fltarr(n1, 1, 2)
            ltc2[0:nltc2-1, *, *] = temp
            delvarx, temp
         ENDELSE
;Sum for the variances         
         FOR i=0, 1 DO BEGIN
            lc0 = reform(reform(ltc2[*, 0, i]), vnbin, ntimes_rates)
            avcts = avsig_1(lc0, dimension=1, sigma)
            ok = where(avcts GT 0 AND sigma GT 0)
            IF(ok[0] NE -1) THEN rates1[det[i]-7, ok] = $
                  byte(round(10.0*(sigma[ok]^2/avcts[ok])))
         ENDFOR
      ENDIF
;you need to synch the absolute time range up with the oti time range
      ntimes_rates1 = N_ELEMENTS(rates1[0, *])
      rates_mvar = fltarr(2, ntimes_rates)
      tr1 = oltc2 -> get(/absolute_time_range)
      tr1 = ref_time+dt_rates*round((tr1-ref_time)/dt_rates)
      IF(total(abs(oti-tr1)) EQ 0.0 AND ntimes_rates EQ ntimes_rates1) $
         THEN rates_mvar = temporary(rates1) ELSE BEGIN
;to do the synching correctly, embed rates1 in a larger array
         txmin = min([oti, tr1], max=txmax)
         nx = long((txmax-txmin)/dt_rates)
         ratesx = fltarr(2, nx)
;x1=1st ss of rates, in ratesx
;y1=1st ss of rates1, in ratesx
         y1 = long((tr1[0]-txmin)/dt_rates)
         ratesx[*, y1:y1+ntimes_rates1-1] = rates1
         x1 = long((oti[0]-txmin)/dt_rates)
         rates_mvar = ratesx[*, x1:x1+ntimes_rates-1]
         delvarx, rates1, ratesx
      ENDELSE
      t01 = !Stime
      IF(notquiet) THEN BEGIN
         print, 'Allrates has mod_variance:  ', t01
         t01 = anytim(t01)
         Elapsed_time = t01-t00
         per_pkt = (elapsed_time/npkt)*1000.0
         print, 'Elapsed time: ', elapsed_time, '  s'
         print, 'Per 1000 Packets: ', per_pkt, '  s/1000pkt'
      ENDIF
   ENDIF ELSE BEGIN
      IF(notquiet) THEN message, /info, 'No Event Packets (app_id=100)'
   ENDELSE
   
;Particle rates
   ok = where(pkttbl.app_id EQ 102)
   IF(ok[0] EQ -1) THEN BEGIN
      message, /info, 'No Monitor rate packets'
   ENDIF ELSE BEGIN             ;got packets
      mpak = pkt[ok]
      IF(keyword_set(endfileinfo)) THEN BEGIN
         IF(datatype(endfileinfo.monitor_packet) EQ 'STC') THEN $
            mpak=[endfileinfo.monitor_packet, mpak]
      ENDIF
;Use hsi_monitor_rate_read directly
      mr_ut_ref = 0
      hsi_monitor_rate_read, packet_array=mpak, mr, nmr, dtmr, $
         ut_ref=ut_ref_mr, time_array=mr_time_array
;Use 
      help, mr, /str
      nmr = N_ELEMENTS(mr.particle_lo)
;Interpolate to the interval times
      tim_arr_rates = ut_ref_rates+dt_rates*dindgen(ntimes_rates)+dt_rates/2.0
      part_rates = fltarr(2, ntimes_rates)
      FOR i=0, 1 DO BEGIN
         IF i EQ 0 THEN temp = reform(mr.particle_lo, nmr) $
         ELSE temp = reform(mr.particle_hi, nmr)
         part_rates[i, *] = interpol(temp, mr_time_array, tim_arr_rates)
      ENDFOR
;Need livetime, too
      nseg = n_elements(mr[0].live_time)
      lvtim = fltarr(ntimes_rates, nseg)
      
      FOR i=0, nseg-1 DO BEGIN
         lvtim[*, i] = interpol(mr.live_time[i]+1.0/256.0, $ 
                              mr_time_array, tim_arr_rates)
         lvtim[*, i] = livtm[*, i] < 1.0
         lvtim[*, i] = livtm[*, i] > 1.0e-20 ;min value to avoid overflows
;Divide rates by lvtim
         IF(rates[0] NE -1) THEN $
            FOR k=0, nenbins-1 DO rates[*, k, i]=rates[*, k, i]/lvtim[*, i]
      ENDFOR
      av_livtim = total(livtim, 2)/float(nseg)
   ENDELSE
  
;total rates, divide by the number of front segments used
   IF(rates[0] NE -1) THEN BEGIN
      used = where(segi EQ 1, nused)
      IF(nused EQ 0) THEN BEGIN
         message, /info, 'No Detectors?'
         ndet = 0
      ENDIF ELSE BEGIN
         det_index_l = where(segi[0:8] EQ 1, ndet_l)
         IF(det_index_l[0] NE -1) THEN ndet = ndet_l ELSE ndet = nused
         IF(nused GT 1) THEN BEGIN
            tot_rates = total(rates, 3)/ndet
         ENDIF ELSE tot_rates = rates
      ENDELSE
   ENDIF
   
   t01 = !Stime
   IF(notquiet) THEN print, 'Allrates End:  ', t01
   t01 = anytim(t01)
   Elapsed_time = t01-t00
   per_pkt = (elapsed_time/npkt)*1000.0
   IF(notquiet) THEN BEGIN
      print, 'Elapsed time: ', elapsed_time, '  s'
      print, 'Per 1000 Packets: ', per_pkt, '  s/1000pkt'
      help, /memory
   ENDIF
   
   hsi_rates_struct = {pkt:pkt, pkttbl:pkttbl, pktt:pktt, $ 
                       rates:rates, tot_rates:tot_rates, $
                       rates_mvar:rates_mvar, ut_ref_rates:ut_ref_rates, $
                       dt_rates:dt_rates, ntimes_rates:ntimes_rates, $
                       part_rates:part_rates, rates_ebands:rates_ebands, $ 
                       av_livtim:av_livtim, seg_index_mask:segi, $
                       ndet:ndet}
   RETURN
END