;+
; PROJECT:
; HESSI
; NAME:
; HSI_CALIB_EVENTLIST__DEMODULATE
;
; PURPOSE:
; This procedure tries to determine the intrinsic source
; time history, i.e. to demodulate the effect of the grids.
;
; CATEGORY:
;
;
; CALLING SEQUENCE:
;
;
; CALLS:
; none
;
; INPUTS:
;
;
; OPTIONAL INPUTS:
; none
;
; OUTPUTS:
;
;
; OPTIONAL OUTPUTS:
;
;
; KEYWORDS:
; none
; COMMON BLOCKS:
; none
;
; SIDE EFFECTS:
; none
;
; RESTRICTIONS:
; none
;
; PROCEDURE:
; This procedure attempts to separate the true source time variation from
; the modulation produced by the grids as the spacecraft rotates. A blend of two
; strategies is used for each detector in the calibrated eventlist.
;   1. For each det_index/harmonic, the observed flux is
; correlated against the time profile for a source at the FOV center allowing
; for differences in phase.  The profile with the highest  correlation is saved for each
; det_index.  Then smoothed profiles, using the smoothing time, are created for the actual
; count flux and the central profile. Then, the real smoothed data are divided by the smoothed model,
; and that profile is saved for each det_index.
; 2.
; A smoothed time profile is created for each det_index,
; using the smoothing time.
;
; The demodulated and smoothed profiles are averaged together over all of the
; set elements of det_index and used for the FLUX_VAR template for each det_index.
;
; MODIFICATION HISTORY:
; 16-aug-2001, richard.schwartz@gsfc.nasa.gov
; 29-nov-2001, Use ptr_valid to select detectors, ras
; 1-jul-2002, ras, prevent nan from division by 0.
; 19-apr-2007, ras, modify to allow harmonics
;-


pro hsi_calib_eventlist_raw::demodulate, cbe, nflux=nflux

caller = self.caller

 ;
 ;Estimate the true flare morphology by summing and filtering
 ;the rates observed with all 9 detectors, limited to selected a2d.

 det_index_mask = bytarr(9)+1b ;caller->get(/det_index_mask)

 seldet=where( ptr_valid(cbe[*,0])*det_index_mask, ndet)

 if caller->Get(/use_flux_var) ne 1 then begin
     for i=0,ndet-1 do (*cbe[seldet[i]]).flux_var = 1.0
     return
 endif

 user_flux_var = Caller->Get(/user_flux_var)
 if n_elements(user_flux_var) gt 1 then begin
     flux_var = user_flux_var
     flux_var = flux_var / avg(flux_var)
     for i=0,ndet-1 do begin

         idet = seldet[i]
         (*cbe[idet]).flux_var = congrid( flux_var, $
         n_elements( (*cbe[idet]).flux_var ), /interp)
         endfor
     return
     endif

 time_unit = caller->get(/time_unit )
 smoothing_time = caller->get(/smoothing_time)



 checkvar, nflux,  128
 time_bin_min = caller->get(/time_bin_min)
 flux_var_template = fltarr(nflux)
 corrdet =replicate( {corr: 0.0, $
 flux_demod:flux_var_template*0.0+1.0, $
 flux_smth:flux_var_template*0.0+1.0}, ndet)
 nphase = 40
 phase = (findgen(nphase)+0.5)/nphase * !pi * 2.0
 corrsave = 0.0
 ysav = 0.0
 for i=0,ndet-1 do begin
     correl = fltarr(nphase)
     idet = seldet[i]
     temp_cos =(*cbe[idet]).livetime * (*cbe[idet]).gridtran* $
     (1+(*cbe[idet]).modamp[0]*cos( (*cbe[idet]).phase_map_ctr[0]))
     temp_sin =(*cbe[idet]).livetime *(*cbe[idet]).gridtran* $
     (1+(*cbe[idet]).modamp[0]*sin( (*cbe[idet]).phase_map_ctr[0]))
     time_bin = time_bin_min * (caller->get(/time_bin_def))[idet]
     nsmooth  = 2.^20/time_bin * smoothing_time + 1 > 3

     corrsave = 0.0
     phasesave = 0
     ysave= temp_cos
     cnt = (*cbe[idet]).count *1.0
     for j=0,nphase-1 do begin
         ymod = temp_cos*cos(phase[j]) - temp_sin*sin(phase[j])

;YMOD is the time profile for the map center, for each of 40 phases (index j)
;CNT is the observed time profile for this binning.
         correl = correlate( ymod, cnt)
         if correl gt corrsave then begin
             ysave = ymod
             corrsave = correl
             phasesave = j
             endif
         endfor
     ;
     ;Save the best correlation profile.
     corrdet[i].corr = corrsave

   ;Demodulate by smoothing both the observed and model profiles.
   ;Divide the smoothed model into the smoothed actual counts.
   ;Interpolate them into nflux bins using congrid.
     corrdet[i].flux_demod =  congrid( f_div( $
     smooth( cnt,nsmooth<(n_elements( cnt ) -1),/edge_trunc),$
     smooth( ysave,nsmooth<(n_elements( cnt ) -1),/edge_trunc)), nflux)
     corrdet[i].flux_demod = f_div(corrdet[i].flux_demod , avg( corrdet[i].flux_demod))

   ;Smooth the observed counts and interpolate into nflux bins.
     corrdet[i].flux_smth  =  congrid( smooth( cnt, $
     nsmooth<(n_elements( cnt )-1),/edge_trunc), nflux)

     corrdet[i].flux_smth = f_div(corrdet[i].flux_smth,  avg( corrdet[i].flux_smth))
     endfor

 ;Form the final smoothed profile by averaging over the smoothed and demodulated profiles
 ;with equal weighting.
 flux_var = avg( corrdet.flux_smth+ corrdet.flux_demod, 1)/2

 ;Normalize the averaged smoothed demodulated profile and put in FLUX_VAR
 ;
 flux_var = flux_var / avg(flux_var)

 ;Interpolate FLUX_VAR into the flux_var tag for each det_index
 for i=0,ndet-1 do begin

     idet = seldet[i]
     (*cbe[idet]).flux_var = congrid( flux_var, $
     n_elements( (*cbe[idet]).flux_var ), /interp)
     endfor

 end