PRO hsi_visibility_analysis, oce, subcoll, DPHZ=dphz, PS=ps, NRB=nrb, HARM=harm, NFB=nfb, ATTEN=atten, RATE_PLOT=rate_plot
;
; Converts a calibrated eventlist object into a set of calibrated visibitilities.
;
; 13-Dec-04 	Adapted from original vis_test. (ghurford@ssl.berkeley.edu)
;				Hardwired for subcollimator 6.
; 11-Mar-05 gh	Add SUBCOLL and HARM keywords & remove hardwiring
;				Add provision to handle livetime=0 bins.
; 14-Mar-05 gh	Streamline arguments per RS suggestions
; 15-Mar-05 gh	Minor reordering of plotted output.
;				Improve handling of empty roll/phase bins
; 22-Mar-05 gh	Use EJS routine, hsi_vis_fit to convert stack to visibilities using LS fitting.
; 23-Mar-05 gh	Add handling for multiple subcollimators
;				Add a kluge for SAVEing visibility output.
;				Renormalize visibilities to allow for time bin size.
;				Include grid_orientation in call to hsi_vis_display.
; 14-Apr-05 gh	Add correction for front segment detector efficiency.
; 15-Apr-05 gh	Use auto_time_bin=1 with cbe_digital_quality=0.95 instead of hardwired time bins
;				Determine attenuator state automatically
; 19-Apr-05 gh	Reorganize output visibility structure to include descriptive data
; 20-Apr-05 gh	Debug handling of cases where some subcollimators are missing
; 21-Apr-05 EJS Fixed OS-dependent statements to make compatible with unix
;               Fixed error in call to message.
; 22-Apr-05 gh	Replace call to hsi_vis_fit with call to hsi_visibility_fit to accomodate error-weighted fits
; 25-Apr-05 gh	Add RATE_PLOT switch to provide optional plots of stacked data.
; 26-Apr-05 gh	Use reset_plot to restore plotting device.
; 28-Apr-05 gh	Minor cosmetic changes
;
detector_label	= ['1f', '2f', '3f', '4f', '5f', '6f', '7f', '8f', '9f', 	$
						'1r', '2r', '3r', '4r', '5r', '6r', '7r', '8r', '9r']
twopi 			= 2 * !PI
root2 			= 2.^0.5
bignumber 		= 1.E8				; arbitrary big number to serve as a flag
nsc 			= N_ELEMENTS(subcoll)
;
 OPTION = 2		; =0,1,2 for original hsi_vis_fit, hsi_vis_fit without errors, or hsi_vis_fit with errors
IF KEYWORD_SET(ps) 			THEN set_plot, 'ps'			; Set up a single .ps file for all subcoll
IF KEYWORD_SET(rate_plot) 	THEN !P.MULTI	= [0,3,3]
IF N_ELEMENTS(nrb) EQ 0 	THEN nrb 		= 32			; default = 32 roll bins
IF N_ELEMENTS(nfb) EQ 0 	THEN nfb 		= 12
;
det_index 						= subcoll-1
oce-> set, use_auto_time_bin	= 1
oce-> set, cbe_digital_quality 	= 0.95
;
; Define structure to hold the visibility output for multiple subcollimators.
IF N_ELEMENTS(nrb) EQ 1 THEN nrbmax = nrb ELSE nrbmax = MAX(nrb)
vis = { x: 				FLTARR(9, nrbmax), 		$
		y: 				FLTARR(9, nrbmax), 		$
		dc: 			FLTARR(9, nrbmax), 		$
		xyerr:			FLTARR(9, nrbmax), 		$		; = quadratic average of errors in x and y
		chisq:			FLTARR(9, nrbmax),		$
		posn_angle: 	FLTARR(9, nrbmax), 		$		; in degrees
		ok: 			INTARR(9, nrbmax)-1,	$		; preset to -1 ==> visibility missing
		time_range: 	DBLARR(2),		   		$
		energy_band: 	FLTARR(2) }
;
; Calculate the calibrated event list for all subcollimators.
cbe_ptr 	= oce->getdata()															; = [9,3] array of pointers
time_bin 	= oce->get(/time_bin_min) * oce->get(/time_bin_def)		; Get array of time bin durations (bmicrosec)
obs_time	= oce->get(/obs_time_interval)
energy_band = oce->get(/energy_band)
atten1 		= hsi_atten_state(ANYTIM(obs_time[0], /ECS))			; returns attenuator state (0,1,2,3)
atten2 		= hsi_atten_state(ANYTIM(obs_time[1], /ECS))
IF atten1 NE atten2 OR atten1 LT 0 OR atten1 GT 3 THEN MESSAGE, STRCOMPRESS(atten1)+STRCOMPRESS(atten2)
;
; Preload non visibility parameters, whether or not there is data.
vis.time_range	= obs_time
vis.energy_band = energy_band
grid_parameters = hsi_grid_parameters()			; --> structure with grid parameters
;
; Begin loop over requested subcollimators.
; N.B.  Assumes front segments and equal number of roll and phase bins for all subcollimators.
FOR n=0, nsc-1 DO BEGIN
	isc 					= subcoll[n]-1					; this subcollimator index
	hessi_drm_4image, drm, energy_band, detector_label[isc], atten1		; get drm value
	timebinseconds			= time_bin[isc] / 2.^20						; duration of time bin (seconds)
	timearea 				= timebinseconds * drm[0] * 40.				; = timebin(seconds) x effective area(cm^2)
	grid_orientation 		= grid_parameters[isc].orient * !RADEG					; grid orientation (degrees)
; Convert roll angle index to grid position angle.
; See figure 4 in imaging concepts paper, Sol Phys 210,71.
; Roll angle is the direction of the spacecraft Y axis, measured clockwise from solar North
; Grid orientation is the direction of slits with respect to s/c Y axis, positive counterclockwise looking at front of grids.
; Position angle of grids (parallel to lines of maximum transmission) is defined as increasing clockwise from solar north.
; Note that the conventional polar angle in xy plots is +ve counterclockwise from the +ve X axis = 90 degrees - posn angle
	roll_angle				= FINDGEN(nrb)*360./nrb
	pa 						= (roll_angle + grid_orientation) MOD 360.
	ipa 					= SORT(pa)		; sequence of indices to have position angles monatonically increase
	vis.posn_angle(isc,*)	= pa[ipa]
	;
	; Get calibrated event list for this subcollimator
	cbe 					= (*cbe_ptr[isc])								; cbe for selected subcoll
	cbe.phase_map_ctr 		= cbe.phase_map_ctr + dphz						; Correct phase if necessary
	;
	; Stack the time bins, calculate visibilities and load output structure
	IF nfb EQ 0 THEN stack = hsi_phz_stacker( cbe, N_ROLL_BINS=nrb,                    EMPTY_FLAG=-1, /REFORM) $
				ELSE stack = hsi_phz_stacker( cbe, N_ROLL_BINS=nrb,  N_PHASE_BINS=nfb, EMPTY_FLAG=-1, /REFORM)
	IF KEYWORD_SET(rate_plot) THEN BEGIN
		rate 	= F_DIV(stack.count>0, stack.livetime>0)
		maxbin	= N_ELEMENTS(stack.count)
		PLOT, rate/timebinseconds, PSYM=10,  XRANGE=[0, maxbin], XSTYLE=1,	$
					XTITLE='Phase/Roll Bin', YTITLE='Counts / second'
		maxrate = MAX(rate/timebinseconds)
		XYOUTS, maxbin*0.9, maxrate*0.9, STRING(subcoll[n], FORMAT="(I1)")
	ENDIF
	;
	; branch to use old hsi_vis_fit
	IF option EQ 0 THEN BEGIN
		visib					= hsi_vis_fit(stack)
		vis.x [isc, 0:nrb-1] 	= visib[0, ipa] / timearea
		vis.y [isc, 0:nrb-1]	= visib[1, ipa] / timearea
		vis.dc[isc, 0:nrb-1]	= visib[2, ipa] / timearea
		vis.ok[isc, 0:nrb-1]	= 1
	ENDIF ELSE BEGIN
		visib	= hsi_visibility_fit(stack, NOERR=2-option)
		vis.x [isc, 0:nrb-1]	= visib.x [ipa] / timearea
		vis.y [isc, 0:nrb-1]	= visib.y [ipa] / timearea
		vis.dc[isc, 0:nrb-1]	= visib.dc[ipa] / timearea
		vis.ok[isc, 0:nrb-1]	= 1-visib.errflag[ipa]					; temporary kluge to match conventions
		vis.xyerr[isc,0:nrb-1]	= SQRT(0.5*(visib.sigmax[ipa]^2 + visib.sigmay[ipa]^2)) / timearea
		vis.chisq[isc, 0:nrb-1]	= visib.chisq[ipa]
	ENDELSE
ENDFOR
;
; Now display, use and save the resulting visibilities
hsi_vis_display, vis, SUBCOLL=subcoll
hsi_vis_avg_size, vis
hsi_vis_consistency_check, vis
SAVE, vis, FILENAME='vis.dat'				; Temporary expedient
;
reset_plot									; Restores plot device to screen
RETURN
END