;+
; PROJECT:
;   HESSI
; NAME:
;   HSI_DECIM_CORRECTION
;
; PURPOSE:
;   This function returns a decimation correction table for all energies, front segments, and
;   weight/channel states
;
; CATEGORY:
;   HESSI SPECTRA
;
; CALLING SEQUENCE:
;
;   correction = hsi_decim_correction( energy2d, weight, channel, coef)
;
; CALLS:
;   none
;
; INPUTS:
;       Energy2d - energy in 2xN form. In keV.  Obtain from spec_obj->getaxis(/energy, /edges_2)
;     Weight   - Weight state for decimation.  Use this value to renormalize for pha channels
;      less than Channel
;     Channel  - PHA channel below which data is decimated by weight
;      weight and channel are obtained from spec_obj->decim_time()
;          i.e. decim_time = spec_obj->decim_time_range()
;     Coef     - obtained from spec_obj->Get_offset_gain()
;
; OPTIONAL INPUTS:
;   none
;
; OUTPUTS:
;       none explicit, only through commons;
;
; OPTIONAL OUTPUTS:
;   none
;
; KEYWORDS:
;   none
; COMMON BLOCKS:
;   none
;
; SIDE EFFECTS:
;   none
;
; RESTRICTIONS:
;   none
;
; PROCEDURE:
;   none
;
; MODIFICATION HISTORY:
;   23-jun-2003, richard.schwartz@gsfc.nasa.gov
;   19-nov-2003, richard schwartz, protect against extreme situations
;     in building decimation correction
;   5-may-2004, added patch to ensure all channels below energy boundaries
;   do have full decimation weighting
;   5-may-2005, ras, make changes to support rear decimation correction
;   25-jul-2005, ras, handle case of degenerate wt/energy better
;   1-Feb-2007, Kim Tolbert, fix catch. previously fell through to run routine again
;   16-feb-2007,ras change gt to ge when looking for fully weighted bins.
;-




function hsi_decim_correction, energy2d, decim_str, coef, $
    front=front, interp_spectrum=interp_spectrum, $
    apar=apar, pflux=pflux, etab=etab, frac=frac, ffff=ffff,$
    error_message = error_message

;Establish an error handler
;
error_message = ''

Catch, error_status
If Error_status NE 0 then begin
    Catch,/cancel
    ;Print, 'Error index: ',error_status
    ;Print, 'Error message:',!Error_state.msg
    ;error_message =  'HSI_DECIMATION_CORRECTION ERROR: Disabling decimation correction ;
    ;Print, error_message
    message, !Error_state.msg + '  Aborting decimation correction.'
    endif

channel = decim_str.channel
weight  = decim_str.weight
atten   = decim_str.atten
ch_wt_atten_mask = [1000L, 10L, 1L]
ch_wt_atten = [[channel],[weight],[atten]]#ch_wt_atten_mask
sel_decim_uniq = uniq( ch_wt_atten, sort( ch_wt_atten) )
decim_uniq =  decim_str[ sel_decim_uniq]

ndecim_uniq = n_elements( decim_uniq ) ;need a correction spectrum for each front seg for each ndecim_uniq
default, front, 1
front = byte(front)
rear = 1-front
default, apar, [1.,2.,1.,4.,400.,4.]
all_atten = get_uniq(atten) ;atten[ uniq( atten, sort(atten))]
natten = n_elements( all_atten )
sel9 = indgen(9)+rear*9
;If front then begin
    cf = float(coef)
    etab=cf[sel9,1]#(channel+1) + cf[sel9,0]#(channel*0.+1.) > 0.0
    el0 = n_elements(energy2d[0,*]) eq 1? 0*etab : value_locate( energy2d[0,*], etab )
    el = el0 >0
    ;create a table for all
;    below_etab = where(  le $
;        rebin( etab, 9, n_elements(energy2d[0,*])), nbelow)


;Before building a count rate spectrum for interpolation, see if we really need one

    if n_elements(interp_spectrum) eq 0 and $
        (min(energy2d) lt max(etab)) then begin
       ;Find the limits on the channels that contain etab

         e2d = minmax(energy2d[*,el0]) + [-3.,13.] >1.01
         ;In case we're dealing with uncorrectable large energy bands
         ovrscl = long([e2d[1]-e2d[0]]/100. + 1)>1.<100.
         eres =(ovrscl/2.)
         nres = fix((((e2d[1]-e2d[0])/eres)>1)[0])
         e2d = get_edges( e2d[0]+(e2d[1]-e2d[0])*findgen(nres+1)/nres, /edges_2)
         ne2d = n_elements( e2d ) / 2
         default, pflux, f_vth_bpow(e2d, apar)
         pflux = pflux * get_edges( e2d, /width)
         asrm = fltarr( ne2d, ne2d, 9, natten)
         for j=0, natten-1 do begin
          hessi_build_srm, get_edges(e2d,/edges_1),[bytarr(9)+front,bytarr(9)+rear], $
              temp_srm, /pha_on_row, all_simp=1,atten_state = all_atten[j], /sep_det
              asrm[0,0,0, j] = temp_srm
              endfor
         interp_spectrum = fltarr( ne2d, 9, natten )
         for i=0,8 do for j=0,natten-1 do interp_spectrum[0,i,j] = asrm[*,*,i,j] # pflux

         endif

       ufrac = fltarr(9,ndecim_uniq)
       for j=0, ndecim_uniq-1 do begin
         sel_etab = (where( decim_uniq[j].channel eq channel))[0]
         els = el[*,sel_etab]
         sel_atten = (where( decim_uniq[j].atten eq all_atten))[0]
         etabs = etab[*,sel_etab]
         for i=0,8 do begin
          temp = etabs[i] gt energy2d[0,els[i]] and etabs[i] lt energy2d[1,els[i]] ? $
          interp2integ( [energy2d[0,els[i]], etabs[i], energy2d[1,els[i]] ], $
          e2d, interp_spectrum[*, i, sel_atten] ) : 0
          ufrac[i,j] = f_div(temp[0] , total(temp))
          endfor
         endfor
         ;Expand frac
         ;else    frac = (etab - energy2d[0,el])/(energy2d[1,el]-energy2d[0,el])
    ;
    nel = (size(/dim,energy2d))[1]



    weight_arr = weight ## (intarr(9)+1.)
    nfrac = n_elements( weight_arr )
    frac = weight_arr * 0.0
    these_decim = ndecim_uniq eq 1 ? intarr(n_elements(ch_wt_atten)) : $
        value_locate( ch_wt_atten[sel_decim_uniq], ch_wt_atten)

    for i=0,n_elements(weight)-1 do frac[0,i] = ufrac[*, these_decim[i]]
    ffff = 1/( frac * 1/weight_arr + 1 * (1-frac))
    correction = fltarr( nel, nfrac ) + 1.
    correction = reform(/over, correction, nel*9L, nfrac/9 )
    ;This patch for now until we can come up with a better solution
    ;It will find all channels, for all decim states, for all a2d that
    ;are completely below the boundary energy

    for i=0,nfrac/9 -1 do begin
       use_full_weight= where( transpose( reproduce(etab[*,i], nel)) ge $;ras, 16-feb-2007, change gt to ge
         reproduce(reform(energy2d[1,*]), 9), nfull)
       if nfull ge 1 then correction[use_full_weight, i] = weight[i]

        endfor


    select = where( correction gt correction[1:*], nselect)
    if nselect ge 1 then begin
        ij = get_ij( select, nel )
        select = nel*ij[1,*] + ij[0,*]+1
        correction[select] = ffff[ij[1,*]]
        endif

    correction = reform(/over, correction, nel, 9, nfrac/9  )
;    for i=0,n_elements(el)-1 do begin
;       if el[i] gt 0 then correction[0:el[i]-1,i] =weight_arr[i]
;       if el0[i] ge 0 then correction[el[i],i]= ffff[i] >1.0
;       endfor




    ;endif

return, float( correction)
end