;+
; PROJECT:
;	HESSI
; NAME:
;	HSI_ANNSEC_BPROJ_WEIGHT_MAP
;
; PURPOSE:
;	This procedure computes the weighting maps for back projection.
;
; CATEGORY:
;	HESSI, UTIL, IMAGE
;
; CALLING SEQUENCE:
;	hsi_annsec_bproj_weight_map, cbe_obj, weight_map, map_ptr=map_ptr
;
; CALLS:
;	HSI_ANNSEC_MAP
;
; INPUTS:
;       Cbe_obj - object containing calibrated eventlist.
;
; OPTIONAL INPUTS:
;	none
;
; OUTPUTS:
;       none explicit, only through commons;
;
; OPTIONAL OUTPUTS:
;	none
;
; KEYWORDS:
;	none
; COMMON BLOCKS:
;	none
;
; SIDE EFFECTS:
;	none
;
; RESTRICTIONS:
;	none
;
; PROCEDURE:
;	Implements the weighting factors given in
;	http://hesperia.gsfc.nasa.gov/~schmahl/bproj_correction/corrbproj_gh.html
;
; MODIFICATION HISTORY:
;	Version 1, richard.schwartz@gsfc.nasa.gov, renamed from hsi_annsec_bproj_weight.
;	26-may-2000.
;	
;
;-
pro hsi_annsec_bproj_weight_map, cbe_obj, this_mpat,  $
	  force=force, $
	this_harmonic=this_harmonic, this_det_index=this_det_index



det_index = fcheck(this_det_index,0)
harmonic  = fcheck(this_harmonic,0)
twopi = 2.0 * !pi
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]


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



rmap_dim  = this_mpat.rmap_dim
phase_ptr = this_mpat.phase_ptr
total_phi = this_mpat.total_phi

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

out2  = out
vrate = f_div((*cbe).gridtran,float((*cbe).livetime)) 

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

;cos_factor = float(cos(dphase)*(*cbe).modamp*vrate)
;sin_factor = float(sin(dphase)*(*cbe).modamp*vrate)*scoef
c2_factor  = cos_factor^2
s2_factor  = sin_factor^2
cs_factor  = cos_factor*sin_factor
c_factor   = cos_factor * vrate
s_factor   = sin_factor * vrate


max_sum = max( (*phase_ptr).isum )
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]
	c2_factor [add_to_these]  =  c2_factor[add_to_these] +  c2_factor[select]
	s2_factor [add_to_these]  =  s2_factor[add_to_these] +  s2_factor[select]
	cs_factor [add_to_these]  =  cs_factor[add_to_these] +  cs_factor[select]
	c_factor  [add_to_these]  =   c_factor[add_to_these] +   c_factor[select]
	s_factor  [add_to_these]  =   s_factor[add_to_these] +   s_factor[select]
	endif
	endfor

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

for ii=0L, nmap-1 do begin
        i = select0[ii]
	index1 = (*phase_ptr)[i].map_index * rmap_dim[0]
	index2 = index1 + npixel
	cmap = (*this_mpat.cmap_ptr)[index1:index2]
        smap = (*this_mpat.smap_ptr)[index1:index2]
	out    = out + cmap * cos_factor[i] - $
				smap * sin_factor[i]
	
	out2 = out2 + 2*(  cmap * c_factor[i] - $
				smap * s_factor[i] )
	
	endfor
;stop 
temp = (*this_mpat.cmap_ptr)^2
for ii = 0L, nmap-1 do begin
        i = select0[ii]
	index1 = (*phase_ptr)[i].map_index * rmap_dim[0]
	index2 = index1 + npixel
	out2 = out2 +  c2_factor[i]*temp[index1:index2] 
	
	endfor

temp = (*this_mpat.smap_ptr)^2
for ii = 0L, nmap-1 do begin
        i = select0[ii]
	index1 = (*phase_ptr)[i].map_index * rmap_dim[0]
	index2 = index1 + npixel
	out2 = out2 +  s2_factor[i]*temp[index1:index2] 
	
	endfor

temp = (*this_mpat.cmap_ptr)*(*this_mpat.smap_ptr)

for ii = 0L, nmap-1 do begin
        i = select0[ii]
	index1 = (*phase_ptr)[i].map_index * rmap_dim[0]
	index2 = index1 + npixel
	out2 = out2 - 2* cs_factor[i]*temp[index1:index2]
	
	endfor



out = out + total(vrate)
out2 = out2 + total(vrate^2)

this_mpat.weight_map_ptr = ptr_new(reform( [[out],[out2]]/(npat),rmap_dim[0], rmap_dim[1],2))

end