;+
; PROJECT:
;	HESSI
; NAME:
;	HSI_ANNSEC_BPROJ
;
; PURPOSE:
;	Returns the back-projection of the hessi calibrated eventlist.
;
; CATEGORY:
;	HESSI, UTIL
;
; CALLING SEQUENCE:
;	map=hsi_annsec_bproj(cbe_obj, map_ptr=fifo, /sum, det_index_mask=det_index_mask)
;	
;	Call it this way when the modulation pattern maps are carried outside of the object.
;	This mode is used for testing.
;
;	map=hsi_annsec_bproj(cbe_obj, map_ptr=map_ptr,/sum)
;
;
; CALLS:
;	none
;
; INPUTS:
;       CBE_OBJ - Object containing the calibrated eventlist or a pointer array to same.  If computing this for 
;		a set of detectors without the THIS_DET_INDEX and THIS_HARMONIC keywords, then either
;		the cbe containing object or the cbe pointer + DET_INDEX_MASK + rmap_dim must be defined.
;
; OPTIONAL INPUTS:
;	
;
; OUTPUTS:
;       Returns the map as a pointer (if SUM eq 0) or as a 2-d image in polar coordinates.
;
; OPTIONAL OUTPUTS:
;	none
;
; KEYWORDS:
;	MAP_PTR - keyword only to be used in testing.  Normally modulation pattern maps are carried
;	in the hsi_modul_pattern object, but for testing purposes (see hsi_annsec_test.pro) they
;	are passed around through this keyword.
;	SUM - If set, sum the map over the detectors.
;	THIS_DET_INDEX- detector index for single collimator/harmonic pair.
;	THIS_HARMONIC- harmonic index for single collimator/harmonic pair.
;	DET_INDEX_MASK - Specify array of collimators and harmonics. See hsi_calib_eventlist.
;	DATA_PTR - Alternate count rate input. 
;		E.G. Used to create point spread functions and intermediate MEM images.
;		Counts per time_bin corrected for live_time. May be a vector if not for
;		an array of detectors.
;	COUNTS_SUMMED - sum of livetime corrected counts. Returned as a 9 x 3 array
;	unles this_det_index is used, when it comes back as a single number.
;	
; COMMON BLOCKS:
;	none
;
; SIDE EFFECTS:
;	none
;
; RESTRICTIONS:
;	none
;
; PROCEDURE:
;	none
;
; MODIFICATION HISTORY:
;	Version 1, richard.schwartz@gsfc.nasa.gov
;	30-mar-2000
;	Version 2, richard schwartz, det_index_mask and a2d_index_mask are anded. 2-may-2000.
;	Version 3, ras, uses new direct indexing (see hsi_annsec_map) and phase_ptr.isum field.
;	Version 3.1, ras, 25-may-2000
;	- Added protection against no data in calib_eventlist for this_det_index & this_harmonic.
;
;-
function hsi_annsec_bproj,  cbe_obj,  $
map_ptr=map_ptr, det_index_mask=det_index_mask, $
sum=sum,  time_test=time_test, weight=weight, $
this_harmonic=this_harmonic, this_det_index=this_det_index, data_ptr=data_ptr, $
counts_summed=counts_summed

if not exist(this_det_index) then begin

    if not exist(det_index_mask) then det_index_mask = cbe_obj->get(/det_index_mask)
    det_index_mask = det_index_mask and cbe_obj->get(/a2d_index_mask)    ;;;;  To be removed after Release 5.

    if not exist(weight) then weight = cbe_obj->get(/weight)
    
    if not keyword_set(sum) then out = ptrarr(9,3) else out = 0
    checkvar, counts_summed, fltarr(9,3)

    a2d_index = Where( cbe_obj->get(/a2d_index_mask) NE 1, count )
    IF count NE 0 THEN BEGIN 
        det_index_mask[a2d_index, *] = 0
    ENDIF
    
    det_list = where( det_index_mask[*,0], n_det)
    
    for idet_index=0,n_det-1 do begin
        det_index = det_list[idet_index]
        harmonic_list = where( det_index_mask[det_index, *], n_harmonic)
        for jharmonic=0, n_harmonic-1 do begin
            harmonic = harmonic_list[jharmonic]
                                ;
            
            bproj = hsi_annsec_bproj( cbe_obj, $
            time_test = time_test, $
            map_ptr=map_ptr, this_det_index=det_index, this_harmonic=harmonic,$
            data_ptr=data_ptr, counts_summed= counts, weight=weight)
            counts_summed[det_index,harmonic]=counts
            if keyword_set(sum) and n_elements(out) eq 1 then begin
		rmap_dim = (*map_ptr[det_index, harmonic]).rmap_dim
		out = fltarr(rmap_dim[0],rmap_dim[1])
		endif
            if keyword_set(sum) then out = out + bproj else $
            out[det_index,harmonic]=ptr_new(bproj)
            endfor
        endfor
    return, out
    endif

t0=systime(1)
det_index = this_det_index
harmonic  = this_harmonic
twopi = 2.0 * !pi

if not exist(weight) and size(/tname,cbe_obj) eq 'OBJREF'  then weight = cbe_obj->get(/weight)

if size(/tname,cbe_obj) eq 'POINTER' then cbe = cbe_obj else $
cbe = cbe_obj->getdata(class='hsi_calib_eventlist')
if n_elements(cbe) gt 1 then cbe = cbe[det_index,harmonic]
cbe_valid = ptr_valid( cbe )

i1=0
i2=0
if cbe_valid then i2=n_elements(*cbe)-1
npat = i2-i1+1

out = 0.0

if not cbe_valid then	return, out

if size(/tname,map_ptr) ne 'POINTER' then $
  map_ptr = cbe_obj->getdata(class='hsi_modul_pattern')

mpat=*map_ptr[det_index,harmonic]

;mpat=hsi_annsec_map(cbe_obj,map_ptr,det_index, harmonic=harmonic, $
;period=period) 


phase_ptr = mpat.phase_ptr
rmap_dim  = mpat.rmap_dim
out = fltarr(long(rmap_dim[0])*rmap_dim[1])


livetime = (*cbe).livetime ;* (cbe[1].time-cbe[0].time) / 2.0^20

dtype = size(/tname,data_ptr)
case 1 of
    dtype eq 'POINTER': data = *data_ptr[det_index,harmonic]
    dtype eq 'UNDEFINED': data  = f_div((*cbe).count,float(livetime))
    else: data = data_ptr
    endcase

vrate      = data * (*cbe).gridtran
counts_summed = total( data )

cos_factor = float(cos((*phase_ptr).dphase)*(*cbe).modamp*vrate)
sin_factor = float(sin((*phase_ptr).dphase)*(*cbe).modamp*vrate)


max_sum = max( (*phase_ptr).isum )

;
; Sum the phases for intervals which share modulation pattern maps.
; I.E. we're adding phasors here.
;
select0 = where( (*phase_ptr).isum eq 0, nmap)

if max_sum ge 1 then for this_sum=1L,max_sum do begin
    select = where((*phase_ptr).isum eq this_sum, nselect)  
    if nselect ge 1 then begin
	map_index = (*phase_ptr)[select].map_index    
	if idl_release(lower=5.3,/inclusive) then find_ix='value_locate' else find_ix='find_ix' 
	sindex0 = call_function( find_ix, (*phase_ptr)[select0].map_index,map_index)
	add_to_these = select0[sindex0]              
    	cos_factor[add_to_these]=cos_factor[add_to_these]+cos_factor[select]
    	sin_factor[add_to_these]=sin_factor[add_to_these]+sin_factor[select]*(*phase_ptr)[select].scoef
	endif
    endfor


npixel = long(rmap_dim[0]*rmap_dim[1])-1

t1=systime(1)
;
; Apply the summed phases to the modulation pattern maps and sum them up to obtain
; the back projection.
;
for ii=0L, nmap-1 do begin
    i = select0[ii]
    index1 = (*phase_ptr)[i].map_index * rmap_dim[0]
    index2 = index1 + npixel
    out  = out + (*mpat.cmap_ptr)[index1:index2] * cos_factor[i] - $
    (*mpat.smap_ptr)[index1:index2] * sin_factor[i]
    endfor
t2=systime(1)

;help,det_index, nmap

tot_vrate = total(vrate)
out = out + tot_vrate
if keyword_set(time_test) then print, t2-t1,(t2-t1)/(t2-t0)

out = reform(/over, out, rmap_dim[0],rmap_dim[1])
if keyword_set(weight) then begin
	wmap_ptr =  mpat.weight_map_ptr
	if not ptr_valid(wmap_ptr) then begin
	
	if size(/tname, cbe_obj) ne 'OBJREF' then return, 0

	hsi_annsec_bproj_weight_map, cbe_obj, mpat,  $
	this_harmonic = harmonic, this_det_index= det_index
	
	(*map_ptr[det_index,harmonic]).weight_map_ptr = mpat.weight_map_ptr
	endif

	out = hsi_annsec_bproj_weight( out, counts_summed, mpat.weight_map_ptr)
	endif

if fcheck(error_weight,0) then message,/info,'Weight map not computed. Cbe_obj must be object.'

return, out
end