;+
;  HSI_POLARDEMO -- a brief script to demonstrate polar mapping
;    Simulates a double source with 2000 photons/s/SC for 6
;    subcollimators, makes 6 polar maps, adds them and displays the sum.
;    The model map is in rectangular coordinates, but dmap is in
;    polar coordinates (1st index=azimuth, 2nd=radius)
;    With the present settings, the 1024x512 map spans 360 degrees of azimuth
;    and 512 arcsec in radius. 
;           Bottom of the map is r=394 asec 
;           Top of the map is r=906 asec
;           Left edge is azimuth=-180 (from the positive x axis)
;           Right edge is azimuth=+180 (from the positive x axis)
;    At the mean radius of 650 arcsec,
;    the azimuth pixels are 4 arcsec (sufficient for det_index GE 3).
;    The radial coverage is changed by a keyword, but the azimuthal
;    coverage is probably cast in concrete.

;    EJS -- Nov 24, 1999
;    EJS -- Jan 13, 2000 --Some corrections, cleaning up


       xy_offset=[650,0]
;      2-asec pixels in the model map: (It's true!  Why, I don't know.
; POINT SOURCES:
       model=fltarr(256,256) & model(128+0,128+0)=100 
       & model(128-40,128-0)=100. ;&  model(128,128+8)=5.
; POINT AND EXTENDED SOURCE:
;      gau=2.^(-shift(dist(5),2,2)^2)  ; gaussian with fwhm=2x2
;      gau8x8=rebin(gau,20,20)  ; fwhm=8x8
;      model(127,127)=total(gau8x8) & model(127-10:127+9,127+20:127+39)=gau8x8

       sim_xy_offset=xy_offset+10               ; src position
       sim_photons=2000
       ;det_mask=intarr(27) & det_mask(det_index)=1
       bin_size=2L^12 ; microseconds  (this is ok for det_index=3-8)
;       bin_size=2L^11 ; microseconds  (this is ok for det_index=2-8)
       num_bins=(2L^22)/bin_size   
       time_range=[0,4]  ; If this is changed, polar_mapper won't work
       o = obj_new( 'hsi_calib_eventlist')
       timebindef = (o->get( /time_bin_def )) * 0 + bin_size
       ce = o->getcalibeventlist( SIM_TIME_range=time_range, $
             TIME_RANGE=time_range, SIM_PHOT=sim_photons, $
             TIME_BIN_DEF=timebindef,SIM_MODEL=model, $
             SIM_XYOFFSET=sim_xy_offset, XYOFFSET=xy_offset,pixel_size=1)  

for det_index=3,8 do begin

       ce_deref = *ce[det_index,0]
       dmap=hsi_polar_mapper(det_index, ce_deref,xy_offset) ;,/plot)

       w=locate_val(dmap,max(dmap)) & R_peak=w(0) & AZ_peak=w(1)
       message,'dmap size='+string((size(dmap))(1))+$
         string((size(dmap))(2)),/inform
       message,'AZ(peak),R(peak)='+string(AZ_peak)+string(R_peak),/inform
       message,'1st index is azimuth (-!pi to !pi), 2nd is  radius',/inform

       print,'Finished map ',det_index
       if (kmap eq 0) then dmap_sum=dmap else dmap_sum=dmap+dmap_sum
endfor ; det_index loop

window,0,xsize=n_elements(dmap(*,0)),ysize=n_elements(dmap(0,*))


tvscl,dmap_sum
dmap_xy=hsi_pmap2xy(dmap_sum,960-(size(dmap_sum))[2],960, box_size=256, pixel_size=2.0) 
tvscl, dmap_xy
stop    ; keep the variables in MAIN

end