;+
; PROJECT:
;	HESSI
; NAME:
;	HSI_ANNSEC_MAP
;
; PURPOSE:
;	This function returns the map structure containing the annular sector
;	universal modulation patterns.
;
; CATEGORY:
;	HESSI, IMAGE
;
; CALLING SEQUENCE:
;	mpat = hsi_annsec_map( cbe_obj, map_ptr, det_index, $
;		harmonic=harmonic)
;
; CALLS:
;	HESSI OBJECTS, HSI_ANNSEC_BPROJ_WEIGHT, DATATYPE
;
; INPUTS:
;	Cbe_obj - 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 0 for first harmonic).
; OUTPUTS:
;       Returns a {hsi_annsec_map_a_str}
;
; OPTIONAL OUTPUTS:
;	none
;
; KEYWORDS:
;	none
; COMMON BLOCKS:
;	none
;
; SIDE EFFECTS:
;	none
;
; RESTRICTIONS:
;	This routine is no longer in the active path. It should not be compiled, and
;	is replaced by hsi_annsec_modul_pattern and hsi_annsec_phase.
;
; PROCEDURE:
;	none
;
; MODIFICATION HISTORY:
;	Version 1, richard.schwartz@gsfc.nasa.gov
;	Version 2, ras, 24-apr-2000.
;	Version 2.1, ras, 26-apr-2000. Fix error in assigning map_indices for roll angles gt pi.
;	Version 3, ras, 8-may-2000. The map_index in phz_ptr is now direct.  For each element of the cbe
;	the associated annular sector map index is given directly by *phz_ptr.map_index.
;	Version 3.1, ras, change map_range to modpat_skip, 10-may-2000.
;	Version 3.2, ras, fixed bug in structure assignment statement caused by fewer elements
;	beyond first set of PI.
;	Version 3.3, ras, identify the map_index for all time intervals at once.
;	Version 3.4, ras, 25-may-2000.  Protect against empty calib_eventlist.
;	Version 3.5, ras, 12-jun-2000.  Add protection for mpat_defined.
;	Version 5, richard.schwartz@gsfc.nasa.gov 17-jul-2000.
;	- switching to interpolating on base maps. A modulation pattern is
;	composed of the weighted sum of four base maps. The phases are adjusted in the two sin/cos pairs to
;	match the phase of the real modulation pattern.  The weighting is assigned to the separate pairs based on
;	the distance of each (in rotation angle) from the rotation angle at the center of the event bin.
;	Modpat_skip is now the overall modpat_skip x time_bin_def for this coll/harm.
;	Version 5.1 richard schwartz 24-oct-2000. Add twopi and take mod !pi again to keep droll positive
;	Aug 21, 2001, ras, remove reference to xaxis_fov, now that is always set to 1.
;	This routine is no longer used with the implementation of the
;	universal annsec routines.
;-

function hsi_annsec_map, cbe_obj, map_ptr, det_index, harmonic=harmonic, $
bproj_weight=bproj_weight, period=period

twopi = 2.0 * !pi
checkvar,harmonic,0


if datatype(map_ptr) ne 'PTR' then map_ptr = ptrarr(9,3)

if ptr_valid(map_ptr[det_index,harmonic]) then begin
    mpat = *map_ptr[det_index,harmonic]
    mpat_defined = 1 * (size(/tname, mpat) eq 'STRUCT')
    endif else mpat_defined = 0

if mpat_defined and datatype(cbe_obj) eq 'OBJ' then begin
    cbe_control=cbe_obj->get(/xyoffset,/pixel_size,/image_dim,/pixel_scale,/factor_by,/modpat_skip,$
    /r0_offset)
    mpat_defined = mpat_defined * $
    (total(abs(mpat.xyoffset    - cbe_control.xyoffset))    eq 0 ) * $
    (total(abs(mpat.pixel_size  - cbe_control.pixel_size))  eq 0 ) * $
    (total(abs(mpat.image_dim   - cbe_control.image_dim))   eq 0 ) * $
    (total(abs(mpat.modpat_skip - cbe_control.modpat_skip)) eq 0 ) * $
    (total(abs(mpat.factor_by   - cbe_control.factor_by))   eq 0 ) * $
    (total(abs(mpat.pixel_scale - cbe_control.pixel_scale)) eq 0 )
    if not mpat_defined then ptr_free, map_ptr[det_index,harmonic]
    endif

if not mpat_defined then begin

    ;
    ;   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    =  (cbe_obj->getdata(class='hsi_calib_eventlist'))[det_index,harmonic]

    ;
    ;   If an empty cbe_ptr is passed, we can do nothing here.  So return.
    if not ptr_valid( cbe_ptr ) then return, 0

    xyoffset   = cbe_obj->get(/xyoffset)

    nbin       = n_elements(*cbe_ptr)
    droll_at_bin      = (*cbe_ptr).roll_angle  ; - (xaxis_fov * phi0_rad)

    roll0             = droll_at_bin[0]

    ;droll_at_bin      = (droll_at_bin - roll0) mod !pi ;now roll is differential from offset, and relative to !pi


    ;;;hsi_annsec_map_info, cbe_obj, det_index, roll0, mpat, phase_mctr, harmonic=harmonic
    hsi_annsec_map_info, cbe_obj, det_index, roll0, mpat, phase_mctr, phi0, harmonic=harmonic

    ;now roll is differential from offset, and relative to !pi
    ;;;droll_at_bin      = ((droll_at_bin - roll0) mod !pi + twopi) mod !pi 	;ras 24-oct-2000
    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.

    phz_ptr = mpat.phase_ptr
    dphi  = mpat.dphi
    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.
    ;
    total_phi = mpat.total_phi
    modpat_skip = mpat.modpat_skip * (cbe_obj->get(/time_bin_def))[det_index] ;;;;;;;;;;;;;;;;;;;;;
    ;
    ; Compute ang_to_map - This gives the map index for an index into the roll_angle
    ;roll_at_map = [dphi*((findgen(total_phi/modpat_skip+1)*modpat_skip)  )]
    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)
    ;stop
    (*phz_ptr).map_index[0] = floor(map_index)
    (*phz_ptr).map_index[1] = (*phz_ptr).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_ptr).map_index[0])
    wcoef[1,*]     = 1.0 - wcoef[0,*]

    (*phz_ptr).map_index = (*phz_ptr).map_index * modpat_skip < (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_ptr)[select].dphase = ((*cbe_ptr)[select].phase_map_ctr mod twopi) - $
    ;		(phase_mctr[(*phz_ptr)[select].map_index]*(*phz_ptr)[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_ptr)[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_ptr)[select].map_index]              ;Phase at center of Sine and Cosine Maps.
    ;	coef   = (*phz_ptr)[select].scoef)                              ;1 or -1 depending on fix(roll_angle/!pi)
    ;	dphase = cbe_phase_map_ctr - phase * coef
    ;       (*phz_ptr)[select].dphase = dphase
    ;
    ;

    pi_change = where(  ((*phz_ptr)[1:*].map_index[0]-(*phz_ptr).map_index[0]) lt 0, num_pi_change)

    if num_pi_change eq 0 then pi_change = [-1,nbin-1] else pi_change = [-1,pi_change, nbin-1]
    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 mod twopi) - $
        (phase_mctr[(*phz_ptr)[select].map_index[icoef]] * scoef[icoef,select])

        endfor
    sincos = ['cos','sin']
    for i=0,1 do for j=0,1 do $
    case j of
        0: (*phz_ptr).phz_coef[i,j] = reform(call_function( sincos[j],dphase[i,*])*wcoef[i,*] )
        1: (*phz_ptr).phz_coef[i,j] = reform(call_function( sincos[j],dphase[i,*])*wcoef[i,*]*scoef[i,*])
        endcase

    ; cos_factor[0,*] = float(cos((*phase_ptr).dphase[0])*(*cbe).modamp*vrate)* (*phase_ptr).wcoef[0]
    ; sin_factor[0,*] = float(sin((*phase_ptr).dphase[0])*(*cbe).modamp*vrate)* (*phase_ptr).wcoef[0]
    ; cos_factor[1,*] = float(cos((*phase_ptr).dphase[1])*(*cbe).modamp*vrate)* (*phase_ptr).wcoef[1]
    ; sin_factor[1,*] = float(sin((*phase_ptr).dphase[1])*(*cbe).modamp*vrate)* (*phase_ptr).wcoef[1]


    ;Now identify all the uniq values of map_index

    hist_map_index = histogram( (*phz_ptr).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_ptr).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_ptr).isum = temp


    mpat.phase_ptr = phz_ptr

    map_ptr[det_index, harmonic] = ptr_new(mpat)

    endif

return, mpat
end