;+
; PROJECT:
;   HESSI
; NAME:
;   HSI_ANNSEC_BPROJ_WEIGHT_MAP
;
; PURPOSE:
;   This procedure computes the flatfielding maps for back projection.
;
; CATEGORY:
;   HESSI, UTIL, IMAGE
;
; CALLING SEQUENCE:
;   hsi_annsec_bproj_weight_map, cbe_obj, this_mpat, $
;     force=force, this_harmonic=this_harmonic, this_det_index=this_det_index
;
; CALLS:
;   HSI_ANNSEC_MAP
;
; INPUTS:
;       Cbe_obj - object containing calibrated eventlist.
;     This_mpat - modulation pattern structure for this detector and this harmonic
;     THIS_DET_INDEX - must be included to specify detector index
;
;
; OPTIONAL INPUTS:
;     FORCE - flatfield map is loaded into this_mpat's structure's flatfield_map_ptr field.
;
; OUTPUTS:
;       Result is loaded into (*this_mpat).weight_map_ptr as an array consisting of
;     <R> and <R^2> where R is the probability of photon transmission for a given
;     pixel and rotation bin.  <> denotes an average over the rotation modulation.
;
; OPTIONAL OUTPUTS:
;   none
;
; KEYWORDS:
;   none
; COMMON BLOCKS:
;   none
;
; SIDE EFFECTS:
;   none
;
; RESTRICTIONS:
;   none
;
; PROCEDURE:
;   Implements the flatfielding factors given in
;   http://hesperia.gsfc.nasa.gov/~schmahl/bproj_correction/corrbproj_gh.html
;   Equation is broken into components.  Unsquared terms obtained by calling hsi_annsec_bproj
;   with suitable weighting.  Squared terms computed by obtaining coefficient products, collecting
;   like terms, summing over coefficients, and then applying summed coefficients to
;   the squared terms and summing the maps.  By collecting terms for like coefficients, the
;   operation of summing weighted sine and cosine maps is minimized.
;
; MODIFICATION HISTORY:
;   Version 1, richard.schwartz@gsfc.nasa.gov, renamed from hsi_annsec_bproj_weight.
;   26-may-2000.
;
; Version 5, richard.schwartz@gsfc.nasa.gov, implemented interpolation of modpat from
; two base cmaps and two base smaps.  27-jul-2000.
; Version 5.1, richard.schwartz@gsfc.nasa.gov, fixed bug in computing terms.  This error
;   was evidenced by the flatfielding function denominator going negative for coarse
;   collimators close to the telescope axis due to an error in one of the cross terms.
; Version 5.2, richard.schwartz@gsfc.nasa.gov - new formulation by GH dropping gridtran term
; from inside the averages.- 19-dec-2000.
; 11-jun-2001 - remove livetime from weight map.  livetime applied at module profile only.
;      richard.schwartz@gsfc.nasa.gov
; 7-jul-2002, ras, changed bsort to ibsort, much faster.
; 23-oct-2002, ras, change loop indices to Long
; 23-apr-2004, ras, substitute a running average for the flat_rate, correction suggested by ghurford

; use a running mean over this angle (360 degreees / smoothing_angle) instead of flat_rate in flatfield map
; If smoothing_angle is lt 2  it has no effect on the flatfield numerator or denominator.
; 26-dec-2004, ras, made a little more robust and removed development message about local average
; 24-apr-2007, ras, only use first harmonic, correct modamp if cbe.modamp has 3 terms
; 22-dec-2011, ras, warning added about too few samples, for num valid samples (>0 livetime) <6
; 23-jun-2017, ras, updated computation of ptot to use hsi_annsec_bproj_1det
; 12-feb-2018, RAS, changed format blocking
; ;if we're calling this here, flatfield must be
; 20-feb-2020, ras,  now using hsi_annsec_bproj_1det in the
; flatfield branch. Should have been changed in Jun 2017 just as in the changed noted above
;-
pro hsi_annsec_bproj_weight_map, cbe_obj, this_mpat,  $
  force=force, $
  smoothing_bins=smoothing_bins, $
  use_rate=use_rate, $
  this_det_index=this_det_index



  det_index = fcheck(this_det_index,0)

  twopi = 2.0 * !pi
  if size(/tname,cbe_obj) eq 'POINTER' then cbe = cbe_obj else $
    cbe = cbe_obj->getdata(class='hsi_calib_eventlist')
  id = det_index
  if n_elements(cbe) gt 1 then cbe = cbe[id]
  structname = tag_names(*cbe,/stru)
  kill_cbe_ptr = 0
  if stregex(structname,/boolean, 'X') then begin ;need to recast harmonic cbe and harmonic mpat as non-harmonic
    ;for the weight map
    newstr = stregex(structname,/boolean,'STACK') ? {hsi_calib_event_stack} : {hsi_calib_event}
    temp   = cbe
    cbe = replicate( newstr, n_elements(*temp))
    struct_assign, *temp, cbe

    cbe.gridtran =  cbe.gridtran * hsi_gridtran_correction(cbe.modamp)
    cbe = ptr_new(cbe)
    kill_cbe_ptr = 1
    cmap = 	*this_mpat.cmap_ptr
    smap =  *this_mpat.smap_ptr
    *this_mpat.cmap_ptr = cmap[*,*,0]
    *this_mpat.smap_ptr = smap[*,*,0]
  endif

  i1=0
  nzro = where( (*cbe).livetime gt 0.0, npat)
  if npat lt 5 then begin
    print,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    print,'******************************************
    print,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    print,'******************************************
    print,'Warning. Too few livetimes for this Grid:',det_index+1
    print,'It may crash in this routine because there are too few samples
    print,'with non-zero livetime.  This may occur because you
    print,'have problems with SAS solution or the energy range
    print,'is out of the valid range for this grid (detector).
    print,'For the best imaging results, remove this grid
    print,'from the image reconstruction.
    print,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    print,'******************************************
    print,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    print,'******************************************

    print,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    print,'******************************************
  endif
  npat0 = n_elements(*cbe)


  ;vrate = (*cbe).gridtran*float((*cbe).livetime)*(*cbe).modamp ## (fltarr(2)+1)
  ;vrate = float((*cbe).livetime)*(*cbe).modamp ## (fltarr(2)+1)   ;19-dec-00 new formulation by gh implem. by ras.
  vrate = (*cbe)[nzro].modamp ## (fltarr(2)+1)         ;19-dec-00 new formulation by gh implem. by ras.
  rmap_dim  = this_mpat.rmap_dim
  npixel   = rmap_dim[0]*rmap_dim[1] - 1 ;used in calculating final index

  phase_ptr = (*(this_mpat.phase_ptr))[nzro]
  coef         = phase_ptr.phz_coef[*,*,0]
  out = fltarr(long(rmap_dim[0])*rmap_dim[1])

  ;load coef with[ [c1, c2],- [s1, s2]]
  for i=0L,1 do coef[*,i,*] = coef[*,i,*] * vrate * (-1)^i

  mpi = phase_ptr.map_index
  cprd = replicate( {cprd, cs_product: 0.0, cs_type:bytarr(2), map_index:lonarr(2)}, 10, npat)

  ;make coefficients of this equation
  ; (c[0]*CM[0]+s[0]*SM[0] + c[1]*CM[1] + s[1]*SM[1])^2
  ; There will be ten independent terms
  ; c[0]^2, c[1]^2, s[0]^2, s[1]^2, 2c[0]c[1], 2s[0]s[1], 2c[0]s[0], 2c[0]s[1],2c[1]s[0],2c[1],s[1]

  ;First record the coefficient products for the cosine and sine coefficients on themselves
  for i=0L,1 do for j=0b,1b do begin
    k = j*2+i
    cprd[k,*].cs_product = cprd[k,*].cs_product + coef[I, J, *]^2
    cprd[k,*].cs_type     = cprd[k,*].cs_type + (bytarr(npat)+1b)## (bytarr(2) + J)
    cprd[k,*].map_index = cprd[k,*].map_index + mpi[i,*] ## (lonarr(2)+1L)

  endfor

  for j=0b,1b do begin
    k = 4 + j
    cprd[k,*].cs_product = cprd[k,*].cs_product + 2.0* coef[0, J, *] * coef[1, J,*]
    cprd[k,*].cs_type     = cprd[k,*].cs_type + (bytarr(npat)+1b)## (bytarr(2) + J)
    cprd[k,*].map_index = cprd[k,*].map_index + [mpi[0,*], mpi[1,*]]

  endfor

  for j=0b,1b do begin

    k = j*2 +  6
    ;cprd[k,*].cs_product = cprd[k,*].cs_product + 2.0* coef[0, J, *] * coef[0, 1-J, *]
    cprd[k,*].cs_product = cprd[k,*].cs_product + 2.0* coef[J, J, *] * coef[J, 1-J, *]

    cprd[k,*].cs_type     =  cprd[k,*].cs_type  + (bytarr(npat)+1b)## [J,1-J]

    ;cprd[k,*].map_index = cprd[k,*].map_index + [mpi[0,*], mpi[0,*]]
    cprd[k,*].map_index = cprd[k,*].map_index + [mpi[J,*], mpi[J,*]]

    k = k +1
    cprd[k,*].cs_product = cprd[k,*].cs_product + 2.0* coef[0, J, *] * coef[1, 1-J, *]
    cprd[k,*].cs_type       = cprd[k,*].cs_type +  (bytarr(npat)+1b)## [J,1-J]
    cprd[k,*].map_index = cprd[k,*].map_index + [mpi[0,*], mpi[1,*]]
  endfor

  select_rev = where( cprd.cs_type[0] gt cprd.cs_type[1])
  cprd[select_rev].cs_type = cprd[select_rev].cs_type[[1,0]]
  cprd[select_rev].map_index = cprd[select_rev].map_index[[1,0]]
  ord = sort( cprd.map_index[1])
  cprd0 = cprd
  cprd = cprd[ord]
  ord = ibsort( cprd.map_index[0])
  cprd = cprd[ord]


  ;First process  c[0]^2, c[1]^2, s[0]^2, s[1]^2 terms

  for cs_type = 0,1 do begin
    sel =  where( cprd.cs_type[0] eq cs_type and cprd.cs_type[1] eq cs_type and $
      cprd.map_index[0] eq cprd.map_index[1], nsel)
    if nsel ge 1 then begin
      if cs_type eq 0 then $
        temp = (*this_mpat.cmap_ptr)^2 else temp = (*this_mpat.smap_ptr)^2

      if nsel ge 2 then $
        sdiff = [-1,where((cprd[sel[1:*]].map_index[0] - cprd[sel].map_index[0] ) $
        ne 0, nsdiff),nsel-1] $
      else begin
        nsdiff = 0
        sdiff   = [-1,0]
      endelse
      for i=0L, nsdiff do begin
        this_map_index = cprd[sel[sdiff[i]+1]].map_index[0]
        coefficient = total( cprd[sel[sdiff[i]+1:sdiff[i+1]]].cs_product)
        index1 = this_map_index * rmap_dim[0]
        index2 = index1 + npixel
        out = out +  coefficient * temp[index1:index2]

      endfor
    endif
  endfor

  ;Nex t process 2c[0]c[1], 2s[0]s[1]

  for cs_type = 0,1 do begin
    sel =  where( cprd.cs_type[0] eq cs_type and cprd.cs_type[1] eq cs_type and $
      cprd.map_index[0] ne cprd.map_index[1], nsel)


    if nsel ge 1 then begin
      if cs_type eq 0 then $
        temp = this_mpat.cmap_ptr  else temp = this_mpat.smap_ptr

      if nsel ge 2 then $
        sdiff = [-1,where((cprd[sel[1:*]].map_index[0] - cprd[sel].map_index[0] ) ne 0, $
        nsdiff),nsel-1] else begin
        nsdiff = 0
        sdiff   = [-1,0]
      endelse
      for i=0L, nsdiff do begin
        this_map_index = cprd[sel[sdiff[i]+1]].map_index[0]
        next_map_index = cprd[sel[sdiff[i]+1]].map_index[1]
        coefficient = total( cprd[sel[sdiff[i]+1:sdiff[i+1]]].cs_product)
        index1_0 = this_map_index * rmap_dim[0]
        index2_0 = index1_0 + npixel
        index1_1 = next_map_index * rmap_dim[0]
        index2_1 = index1_1 + npixel
        out = out +  coefficient*(*temp)[index1_0:index2_0]*(*temp)[index1_1:index2_1]

      endfor
    endif
  endfor


  ;Finally process the remaining cross terms: 2c[0]s[0], 2c[0]s[1],2c[1]s[0],2c[1],s[1]


  sel =  where( cprd.cs_type[0] eq 0 and cprd.cs_type[1] eq 1, nsel)


  if nsel ge 1 then begin

    ctemp = this_mpat.cmap_ptr
    stemp = this_mpat.smap_ptr

    if nsel ge 2 then $
      sdiff = [-1,where((cprd[sel[1:*]].map_index[0]  ne  cprd[sel].map_index[0] )  or $
      (cprd[sel[1:*]].map_index[1]  ne  cprd[sel].map_index[1] ), nsdiff),nsel-1] $
    else begin
      nsdiff = 0
      sdiff   = [-1,0]
    endelse
    for i=0L, nsdiff do begin
      this_map_index = cprd[sel[sdiff[i]+1]].map_index[0]
      next_map_index = cprd[sel[sdiff[i]+1]].map_index[1]
      coefficient = total( cprd[sel[sdiff[i]+1:sdiff[i+1]]].cs_product)
      index1_0 = this_map_index * rmap_dim[0]
      index2_0 = index1_0 + npixel
      index1_1 = next_map_index * rmap_dim[0]
      index2_1 = index1_1 + npixel
      out = out +  coefficient*(*ctemp)[index1_0:index2_0]*(*stemp)[index1_1:index2_1]

    endfor
  endif

  ;vrate = (*cbe).gridtran*float((*cbe).livetime)
  ;vrate = float((*cbe).livetime) ;19-dec-00 new formulation by gh implem. by ras.


  flat_rate = fltarr(npat0)
  flat_rate[nzro] = 1.
  ptot = hsi_annsec_bproj_1det( cbe, map_ptr=this_mpat, this_det_index=this_det_index, $
    data_ptr = flat_rate,/max_harmonic, flatfield=0)

  ;out = out + 2 * hsi_annsec_bproj( cbe, map_ptr=this_mpat, this_det_index=this_det_index, $
  ;this_harmonic=this_harmonic, data_ptr = flat_rate) -  total((flat_rate^2))
  ;Ras, simplified above expression
  out= out + 2 * ptot - npat

  ;here we substitute a running average for the flat_rate, correction suggested by ghurford 23-apr-2004, ras
  default, smoothing_bins, 0
  default, use_smoothing, 0
  if smoothing_bins ge 2 then begin
    ; message,/cont,'USING LOCAL SMOOTHING FOR FLAT RATE'
    inarray = f_div( (*cbe)[nzro].count, (*cbe)[nzro].livetime )

    smoothed_rate =hsi_gsmooth( /edge_truncate, inarray , smoothing_bins )
    new_flat_rate = fltarr(npat0)

    new_flat_rate[nzro] = F_DIV(smoothed_rate ,avg(smoothed_rate)) ;protects against division by 0
;    new_ptot = hsi_annsec_bproj( cbe, map_ptr=this_mpat, this_det_index=this_det_index, $
;      data_ptr = new_flat_rate,/max_harmonic ) ;if we're calling this here, flatfield must not
      ;set and this needs the same code as above with the new_flat_rate instead of the simple flat_rate
    ;changed 20-apr-2020, RAS, 
    new_ptot = hsi_annsec_bproj_1det( cbe, map_ptr=this_mpat, this_det_index=this_det_index, $
        data_ptr = new_flat_rate,/max_harmonic, flatfield=0)
    this_mpat.weight_map_ptr = ptr_new(reform( [[new_ptot[*]],[out],[ptot[*]]]/(npat), $
      rmap_dim[0], rmap_dim[1],3))
  endif else this_mpat.weight_map_ptr = ptr_new(reform( [[ptot[*]],[out]]/(npat), $
    rmap_dim[0], rmap_dim[1],2))

  if kill_cbe_ptr then begin ;cleanup transformation from 3 harm to 1 and back
    ptr_free, cbe
    *this_mpat.cmap_ptr = cmap
    *this_mpat.smap_ptr = smap
  endif

end