;+
; PROJECT:
;   HESSI
; NAME:
;   HSI_annsec_phase
;
; PURPOSE:
;   This function returns the map structure containing the annular sector
;   universal modulation patterns.
;
; CATEGORY:
;   HESSI, IMAGE
;
; CALLING SEQUENCE:
;   mpat = hsi_annsec_phase( self, map_ptr, det_index, $
;     harmonic=harmonic)
;
; CALLS:
;   HESSI OBJECTS, HSI_ANNSEC_BPROJ_WEIGHT, DATATYPE
;
; INPUTS:
;   self - object containing calibrated eventlist object.  Should also contain modpat control object.
;   Map_Ptr - Pointer array which stores the universal modulation pattern map structure.
;
;   Det_index - detector index
;
; OPTIONAL KEYWORD INPUTS:
;   HARMONIC - harmonic index (starts at 1 for first harmonic).
; OUTPUTS:
;
;
; OPTIONAL OUTPUTS:
;   none
;
; KEYWORDS:
;   none
; COMMON BLOCKS:
;   none
;
; SIDE EFFECTS:
;   none
;
; RESTRICTIONS:
;   none
;
; PROCEDURE:
;   none
;
; MODIFICATION HISTORY:
;   4-July-2001 - renamed from hsi_annsec_map.pro, richard.schwartz@gsfc.nasa.gov
;   7-sep-2004 - ras, force roll_angle to be positive
;	19-apr-2007- ras, modify for harmonics
;	21-sep-2010, ras, support modul_pattern saved
; 15-feb-2018, RAS, standard format (ctrl shft F)
;-

pro hsi_annsec_phase, self, mpat_ptr,  error=error

  error = 0

  twopi = 2.0 * !pi
  det_index = (*mpat_ptr).det_index
  gp         = hsi_grid_parameters()
  time_bin_def = (self->get(/time_bin_def))[det_index]
  harmonic  = (*mpat_ptr).harmonic
  total_phi = (*mpat_ptr).total_phi
  rmap_dim  = (*mpat_ptr).rmap_dim
  dphi      = (*mpat_ptr).dphi
  modpat_skip = self->get(/modpat_skip)
  (*mpat_ptr).modpat_skip = modpat_skip
  (*mpat_ptr).xyoffset = self->get(/xyoffset)


  ;
  ;   Here we use the calibrated eventlist to build the universal modulation pattern.
  ;   Of course, the actual counts are irrelevant.  Only the aspect solution and FOV matter.
  ;   This should change later when the calib_eventlist is broken apart into its components.

  ;cbe_ptr    =  self->getdata(class='hsi_calib_eventlist')       ;kim 02-apr-2002
  class='hsi_calib_eventlist'
  oc = self->get(/obj,class=class)
  cbe_ptr    =  oc->framework::getdata()       ;ras, 21-sep-2010

  if total(ptr_valid(cbe_ptr)) eq 0 then begin       ;kim 02-apr-2002
    message,/cont,'No Calibrated eventlist passed.'
    error=1
    return
  endif

  cbe_ptr = cbe_ptr[det_index];,harmonic]     ;ras 19-apr-2007
  ;
  ;   If an empty cbe_ptr is passed, we can do nothing here.  So return.
  if not ptr_valid( cbe_ptr ) then begin
    message,/cont,'No Calibrated eventlist passed.'
    error=1
    return
  endif

  nbin       = n_elements(*cbe_ptr)
  ;Force positive roll_angle, follow on code depends on positive values
  roll_angle = ((*cbe_ptr).roll_angle mod (2.*!pi) )+(2.*!pi)
  droll_at_bin      = roll_angle  ; - (xaxis_fov * phi0_rad)

  ;Check to see if phase needs evaluation
  phz_param= (*mpat_ptr).phz_param
  do_phase = 1 ; leave the question of this update to the object
  ;do_phase = (NOT keyword_set( *(*mpat_ptr).phase_ptr)) or $
  ;   (NOT time_bin_def eq phz_param.time_bin_def) or $
  ;   (NOT modpat_skip  eq phz_param.modpat_skip) or $
  ;   (NOT gp[det_index].orient eq phz_param.gp_orient) or $
  ;   (NOT same_data( droll_at_bin, *phz_param.roll) )

  if NOT do_phase then return

  roll0             = droll_at_bin[0]


  phi0_rad = 0.0 ;provided we always have xaxis_fov set.

  num_fov_r         =  long(rmap_dim[0])
  num_fov_phi       =  long(rmap_dim[1])

  phi0     = phi0_rad - (!pi/2 - gp[det_index].orient)  + roll0  $
    - dphi * ((rmap_dim[1]-1)/2. + (rmap_dim[1] mod 2)/2.0)

  phz  = replicate( {hsi_annsec_phase_str}, nbin )
  ;
  droll_at_bin      = ((droll_at_bin -roll0 + phi0) mod !pi + twopi) mod !pi  ;ras 24-oct-2000
  ipi0 = fix(((phi0 + twopi) mod twopi ) / !pi) ;;; correction for starting segment of circle.

  scoef = intarr(2,nbin)  ;hold the sign associated with rotation_angle mod pi
  wcoef = fltarr(2,nbin)  ;weighting coefficients - two maps linearly weighted for one bin.
  dphase= fltarr(2,nbin)
  ;The number of maps available covers slightly more than pi rotation.
  ;It could be pi, but it would make it more difficult to select the correct map for each bin
  ;because of the endpoint rollover.
  ;

  modpat_skip_4det = modpat_skip * time_bin_def ;;;;;;;;;;;;;;;;;;;;;
  ;
  ; Compute ang_to_map - This gives the map index for an index into the roll_angle

  next_pi = where( droll_at_bin gt droll_at_bin[1:*], num_pi)
  if num_pi eq 0 then next_pi = nbin-1

  map_index = droll_at_bin/(dphi*modpat_skip_4det)
  phz.map_index[0] = floor(map_index)
  phz.map_index[1] = phz.map_index[0] + 1
  ;
  ; The wcoef values are used to interpolate the products of the matrix multiplications for
  ; the blended modulation patterns.  Each sine and cosine map is best expressed as the
  ; weighted sum of its two surrounding maps.  This is always better, and at low computational
  ; cost.  For modpat_skip gt 1, this is particularly important.
  ;

  wcoef[0,*]     = 1.0 - (map_index -  phz.map_index[0])
  wcoef[1,*]     = 1.0 - wcoef[0,*]

  phz.map_index = phz.map_index * modpat_skip_4det < (total_phi-1)


  ;      endfor

  ;
  ;        Computing the phase differences between the centers of the annular sector maps and the phase
  ;   at the calibrated eventlist map center.  With these phase difference we can build the response
  ;   modulation pattern from the sin and cosine maps.
  ;
  ;        if isum eq 0 then phz[select].dphase = ((*cbe_ptr)[select].phase_map_ctr mod twopi) - $
  ;     (phase_mctr[phz[select].map_index]*phz[select].scoef) else $
  ;
  ;   For the each group of rotation angles from n*!pi to (n+1)*!pi, use this expression.  Breaking it apart:
  ;
  ;   dphase = phz[select].dphase             ;Put the results in here
  ;   cbe_phase_map_ctr = ((*cbe_ptr)[select].phase_map_ctr mod twopi)    ;Cal Ev List Map Center
  ;   phase  = phase_mctr[phz[select].map_index]              ;Phase at center of Sine and Cosine Maps.
  ;   coef   = phz[select].scoef)                              ;1 or -1 depending on fix(roll_angle/!pi)
  ;   dphase = cbe_phase_map_ctr - phase * coef
  ;       phz[select].dphase = dphase
  ;
  ;

  pi_change = where(  (phz[1:*].map_index[0]-phz.map_index[0]) lt 0, num_pi_change)
  phase_mctr = *(*mpat_ptr).phase_mctr_ptr

  if num_pi_change eq 0 then pi_change = [-1,nbin-1] else pi_change = [-1,pi_change, nbin-1]
  for hh=0,harmonic-1 do begin ;19-apr-2007, ras
    for ipi = 0,num_pi_change do begin
      select = lindgen( pi_change[ipi+1]-(pi_change[ipi]+1)+1 ) + (pi_change[ipi]+1)
      ;;;scoef[*,select]      =  (-1.)^(ipi)
      scoef[*,select]      =  (-1.)^(ipi+ipi0)


      for icoef=0,1 do $
        dphase[icoef,select] = ((*cbe_ptr)[select].phase_map_ctr[hh] mod twopi) - $
        (phase_mctr[phz[select].map_index[icoef],hh] * scoef[icoef,select])

    endfor
    sincos = ['cos','sin']
    for i=0,1 do for j=0,1 do $
      case j of
      0: phz.phz_coef[i,j,hh] = reform(call_function( sincos[j],dphase[i,*])*wcoef[i,*] )
      1: phz.phz_coef[i,j,hh] = reform(call_function( sincos[j],dphase[i,*])*wcoef[i,*]*scoef[i,*])
    endcase
  endfor ;close harmonic
  ;Now identify all the uniq values of map_index

  hist_map_index = histogram( phz.map_index, min=0, bin=1, max= total_phi-1, reverse_indices= rev_indx)
  slct = where( hist_map_index, nslct )
  ;Need to use temp for isum array because isum has 2 elements for each structure.
  temp = phz.isum
  for i=0,nslct-1 do temp[ rev_indx[rev_indx[slct[i]] : rev_indx[slct[i]+1]-1]] = $
    lindgen( hist_map_index[slct[i]] )
  phz.isum = temp


  *(*mpat_ptr).phase_ptr = phz

  ;
  ;   When initializing the phase pointer, also clear anything out of the
  ;   weight_map_ptr
  ;
  ptr_free, (*mpat_ptr).weight_map_ptr
  phz_param = (*mpat_ptr).phz_param
  phz_param.modpat_skip = modpat_skip
  phz_param.gp_orient      = gp[det_index].orient
  phz_param.time_bin_def= time_bin_def
  *phz_param.roll       = roll_angle
  (*mpat_ptr).phz_param = phz_param

end