;+
; PROJECT:
;	HESSI
; NAME:
;	HSI_ANNSEC_MODUL_PATTERN
;
; PURPOSE:
;	This procedure calculates various annular sector map building blocks.
;
; CATEGORY:
;	HESSI, IMAGE
;
; CALLING SEQUENCE:
;	hsi_annsec_MODUL_PATTERN, self, det_index,  mpat, $
;	  harmonic=harmonic
;
; CALLS:
;	HESSI OBJECTS
;
; INPUTS:
;      self - object containing calibrated eventlist object.
;
;	Det_index     - detector index
;
; OUTPUTS:
;	Mpat - annular sector map structure. partially complete.
;
;
; OPTIONAL OUTPUTS:
;	none
;
; KEYWORDS:
;	HARMONIC - one, default, is first harmonic, 2nd harmonic is 2, 3rd is 3.
;
; COMMON BLOCKS:
;	none
;
; SIDE EFFECTS:
;	none
;
; RESTRICTIONS:
;	none
;
; PROCEDURE:
;	This procedure builds the universal elements of the annular
;	sectors.  These depend only on the FOV dimensions and scale.
;	The phase corrections are computed by hsi_annsec_phase.
;
; MODIFICATION HISTORY:
;	Aug 22, 2001 richard.schwartz@gsfc.nasa.gov
;	Based on the obsolete hsi_annsec_map.pro
; 15-feb-2018, RAS, standard format (ctrl shft F)

;-




pro hsi_annsec_modul_pattern, self, det_index, mpat,   $
  harmonic=harmonic

  twopi = 2.0 * !pi
  checkvar,harmonic,3

  par = self->get(/image_dim,/xyoffset,/pixel_size,/factor_by,$
    /pixel_scale,/modpat_skip,/r0_offset,/xaxis_fov)

  image_dim  = par.image_dim
  xyoffset   = par.xyoffset
  pixel_size = par.pixel_size
  factor_by  = par.factor_by > 1.0
  pixel_scale= par.pixel_scale > 1
  modpat_skip  = par.modpat_skip  > 1
  r0_offset  = par.r0_offset
  xaxis_fov  = 1 ;par.xaxis_fov
  if not xaxis_fov then self->set, /xaxis_fov ;must set xaxis_fov, no longer an option, 4-jul-2001.



  rmap_dim   = hsi_rmap_dim( par )
  ;   pixel_scale must be a factor of both elements of rmap_dim
  ;   pixel_scale is used to set rmap_dim in hsi_rmap_dim

  ;***r0_asec    = sqrt( total(xyoffset^2) )

  if xaxis_fov then r0_asec   = r0_offset

  dphi       = min(pixel_size*pixel_scale) / r0_asec
  total_phi  = ceil(twopi/2.0/dphi) + 48 ;cover at least 1/2 rotation plus some extra for convenience
  ;*** cbe_ptr    = self->getdata(class='hsi_calib_eventlist'))[det_index,harmonic]
  ;*** nbin       = n_elements(*cbe_ptr)
  gp         = hsi_grid_parameters()
  pitch      = gp[det_index].pitch

  ;*** phi0_rad = (atan(xyoffset[1],xyoffset[0]) + twopi) mod twopi
  ;*** if xaxis_fov then phi0_rad = 0.0

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


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

  rmin = r0_asec - num_fov_r/2* (pixel_size[0]*pixel_scale)
  rmax = rmin + num_fov_r * (pixel_size[0]*pixel_scale)

  area_max = !pi * rmax^2
  area_min = !pi * rmin^2
  darea    = (area_max-area_min)/num_fov_r
  darea    = darea < ( !pi*r0_asec^2/(num_fov_r-1 - num_fov_r/2.0 + 0.5))
  area      = area_min+findgen(num_fov_r+1)*darea

  annsec_center = $
    ((findgen(num_fov_r)-num_fov_r/2.+0.5)*darea/!pi + r0_asec^2)^0.5

  ;***

  ;
  ;Compute the area of a pixel, r * dr * dphi
  ;
  area_pixel = avg(interpol(annsec_center[1:*]-annsec_center, $
    findgen(num_fov_r-1)+.5,findgen(num_fov_r))*annsec_center*dphi)


  ;;;annsec_phi = findgen(total_phi + rmap_dim[1]) * dphi + phi_map_first
  annsec_phi = findgen(total_phi + rmap_dim[1]) * dphi ;+ phi_map_first

  wave_number = twopi * (1.+findgen(harmonic)) / pitch

  ;phase should correspond to the centers of fields of view.  The
  ;first field of view starts at rmap_dim[1]/2


  ;The phase at the center of each rotation column. The first index of
  ;each map is for the radius, the second for the rotation around the annulus.
  phase_mctr    = wave_number ## (r0_asec * $
    cos( annsec_phi[ ((rmap_dim[1]+1)/2-1) : ((rmap_dim[1]+1)/2+total_phi-1) ] )) mod twopi
  ;cos( dphi/2.+annsec_phi[ ((rmap_dim[1]+1)/2-1) : ((rmap_dim[1]+1)/2+total_phi-1) ] )) mod twopi

  ;*** phz_ptr  = ptr_new( replicate( {hsi_annsec_phase_str}, nbin ))

  phz_param = {hsi_annsec_phz_param, modpat_skip:modpat_skip,$
    time_bin_def:0L,$
    gp_orient:0.0,$
    roll:ptr_new(0) }
  nc = num_fov_r
  nphi = total_phi + num_fov_phi
  mpat =  ptr_new({hsi_annsec_modul_pattern_str,$
    det_index, $
    harmonic, $
    image_dim,$
    pixel_size,$
    factor_by, $
    pixel_scale, $
    ;The next two parameters are included for convenience. Universal
    ;modulation patterns don't depend on these.
    modpat_skip, $  ;property of phase, not modulation patterns
    xyoffset, $		;property of pointing, not modulation patterns
    rmap_dim, $
    xaxis_fov, $
    r0_offset, $
    total_phi, $
    dphi, $
    ptr_new(annsec_center), $
    ptr_new(annsec_phi), $
    ;ptr_new(reform(phase_mctr[0,*])), $
    ptr_new(phase_mctr), $ harmonic phase info here, total_phi+1 x (harmonic)
    area_pixel, $
    phz_param, $
    ptr_new(0), $ 			;phase_ptr,
    ptr_new(0), $			;weight_map_ptr
    ptr_new(0), $			;cmap_aver_ptr
    ptr_new(0), $			;smap_aver_ptr
    ptr_new(0), $			;cmapb_ptr
    ptr_new(0), $ 			;smapb_ptr
    ptr_new( reform(cos( (cos(annsec_phi)##annsec_center)[*]#wave_number),$
    /over,nc,nphi,harmonic)),$
    ptr_new( reform(sin( (cos(annsec_phi)##annsec_center)[*]#wave_number),$
    /over,nc,nphi,harmonic)) })

end