PRO hsi_vis_display, vis, SUBCOLL=subcoll
;
; PURPOSE: Displays RHESSI visibilities, as calculated by of hsi_vis_fit
;
; INPUT: vis = visibility stucture constaining calibrated visibilities and ancillary data for all subcollimators.
;
; If stack is provided, stacked count rates are displayed; otherwise the DC term is plotted.
;
;   rate = A*cos(phase) + B*sin(phase) + C
;
; 22-Mar-05 	Initial version. (ghurford@ssl.berkeley.edu)
; 23-Mar-05 gh	Add GRID_ORIENTATION keyword
; 14-Apr-05 gh	Plot info on analysis page too.
; 19-Apr-05 gh	Modify to accept revised vis input structure
;				Eliminate keywords, CALIB_OBJ and GRID_ORIENTATION which are no longer needed
; 21-Apr-05 EJS Fixed OS-dependent statements to make runnable on unix
; 28-Apr-05 gh	Adjusted polar plots so that position angle=0 (solar north) corresponds to +ve Y axis.
;
; Some preliminaries
!P.MULTI	= [0,3,2]
nsc = N_ELEMENTS(subcoll)
t1 = ANYTIM(vis.time_range[0], /ECS)
t2 = ANYTIM(vis.time_range[1], /ECS)
allang 		= FINDGEN(361)*!DTOR		; Set of angles for smooth plot
amptab 		= FLTARR(128,9)-1			; Arrays to hold amplitudes and posnangles for replotting. -1==fill
thetatab	= FLTARR(128,9)-1
avamptab	= FLTARR(9)-1
;
; Begin loop over subcollimators
FOR n=0, nsc-1 DO BEGIN
	isc					= subcoll[n] - 1
	pa					= REFORM(vis.posn_angle[isc,*])
	x					= REFORM(vis.x[isc,*])
	y 					= REFORM(vis.y[isc,*])
	xyerr				= REFORM(vis.xyerr[isc, *])
	dc					= REFORM(vis.dc[isc,*])
	amp 				= (x^2 + y^2)^0.5
	phz 				= ATAN(y,x) * !RADEG
	amax				= MAX(amp)*1.1					; add 10% allowance for plotting
	nrb 				= N_ELEMENTS(vis.ok[isc,*])
	avamp				= REPLICATE(MEAN(amp), 361)
	ampp 				= [amp, amp[0]]
	theta				= (90.- [pa, pa[0]])*!DTOR		; Convert position angle to a theta for x-y plotting
	amptab[0:nrb,n]		= ampp							; save for amplitude summary plot
	thetatab[0:nrb,n] 	= theta
	avamptab[n]			= avamp[0]
	xx					= [x, x[0]]
	yy					= [y, y[0]]
;
; Generate plots
	PLOT,  pa,      amp,      PSYM=6, XSTYLE=1,  XRANGE=[0,360], XTITLE='POSN ANGLE (deg)', TITLE='AMPLITUDE'
	OPLOT, pa,    xyerr,            LINESTYLE=0
	PLOT,  pa,      phz,      PSYM=6,  XSTYLE=1,  XRANGE=[0,360], XTITLE='POSN ANGLE (deg)', TITLE='PHASE', 					$
                                       YSTYLE=1,  YRANGE=[-180,+180], YTITLE='PHASE (deg)'
	OPLOT, [0,360],    [0,0],      LINESTYLE=1
	PLOT,  pa, dc,      	  PSYM=-6, XSTYLE=1,  XRANGE=[0,360], XTITLE='POSN ANGLE (deg)', TITLE='DC COMPONENT'
	OPLOT, [0,360],    [0,0],      LINESTYLE=1
	PLOT,   xx,           yy,      PSYM=-6,  XSTYLE=1, XRANGE=[-amax,+amax], TITLE='SINE,COSINE COMPONENTS', 	$
 	                                         YSTYLE=1, YRANGE=[-amax,+amax], /ISOTROPIC
	OPLOT,   [0,0], [-amax,+amax], LINESTYLE=1
	OPLOT, [-amax,amax], [0,0],    LINESTYLE=1
	PLOT,    ampp,      theta,      PSYM=-6,   XSTYLE=1, XRANGE=[-amax,+amax], TITLE='AMPLITUDE, POSN_ANGLE', $
									         YSTYLE=1, YRANGE=[-amax,+amax], /POLAR, /ISOTROPIC
	OPLOT, avamp,     allang,      LINESTYLE=1, /POLAR
;
; Add labels to this plot page and the next
	PLOT, [0,1], [0,1], /NODATA, XSTYLE=4, YSTYLE=4
	XYOUTS, 0, .8, STRMID(t1,  0, 10)
	XYOUTS, 0, .7, STRMID(t1, 11, 8)+'-'+STRMID(t2, 11, 8)
	XYOUTS, 0, .5, STRING(isc+1, FORMAT="('Subcoll', I2)")
	XYOUTS, 0, .3, STRCOMPRESS(STRING(vis.energy_band, FORMAT="(I5,'-', I5)"),/REMOVE_ALL) + ' keV'
	XYOUTS, 0, .1, STRING(nrb, FORMAT="(I3, ' roll bins')")
ENDFOR
;
; Now make a single page of polar amplitude plots
!P.MULTI = [0,3,3]
amax = MAX(amptab)
avmax = MAX(avamptab)
FOR n=0, nsc-1 DO BEGIN
	ok = WHERE(amptab[*,n] GE 0, nok)
	IF nok GT 1 THEN BEGIN
		PLOT, [-amax, amax], [-amax,amax], /NODATA, XSTYLE=1, YSTYLE=1, /ISOTROPIC
;		PLOT,  amptab[ok,n], thetatab[ok,n], PSYM=-6,   XSTYLE=1, XRANGE=[-amax,+amax], TITLE='AMPLITUDE, POSN_ANGLE', $
;									         YSTYLE=1, YRANGE=[-amax,+amax], /POLAR, /ISOTROPIC
		OPLOT,  amptab[ok,n], thetatab[ok,n], PSYM=-6, /POLAR
		OPLOT, REPLICATE(avmax,361),     allang,      LINESTYLE=1, /POLAR
		XYOUTS, 0.8*amax, 0.75*amax, STRING(subcoll[n], FORMAT='(I1)')
	ENDIF
ENDFOR
RETURN
END