PRO hsi_vis_consistency_check, vis, PS=ps
;
; Generates displays to confirm that that amplitudes with equal and opposite position angles are equal and
;   that their phases are complementary.
;
; Input is the visibility structure created by hsi_visibility_analysis
;
; Begin loop over subcollimators
;
; 20-Apr-05 	Initial version (ghurford@ssl.berkeley.edu)
; 21-Apr-05 EJS Fixed OS-dependent statements to make compatible with unix
; 26-Apr-05 gh	Do not restore plotting device.
; 29-Apr-05 gh	Minor plot format change.
;
histphzbin	= 4								; histogram phase bin size (degrees)
maxphz		= 30
nrbmax = 128								; Upper limit to number of roll bins
phzsumtab	= FLTARR(128,9) 				; Will hold phase sums (degrees)
nokps		= INTARR(9)						; Number of ok phase sums for each detector
IF KEYWORD_SET(ps) THEN set_plot, 'ps'
!P.MULTI = [0,2,3]
FOR isc=0,8 DO BEGIN
	dummy 		= WHERE(vis.ok[isc,*] NE 0, nrb)		; ok = -1 if no good, or 1 if ok.  0 if it was never used.
	IF nrb EQ 0 THEN CONTINUE					; no data for this subcollimator
	nrb2 		= nrb/2								; opposite half rotations are separated by nrb2
	ok 			= REFORM(vis.ok[isc,*])
	okp 		= WHERE(ok * SHIFT(ok, nrb2) GT 0, nok)	; indices of partnered visibilities
	IF nok EQ 0 THEN CONTINUE
	x 			= REFORM(vis.x[isc, okp])
	y 			= REFORM(vis.y[isc, okp])
	pa 			= REFORM(vis.posn_angle[isc, okp])
	nok2		= nok/2
	amp 		= SQRT(x^2 + y^2)
	ampp		= SHIFT(amp, nok2)
	amprat		= (amp-ampp) / (amp+ampp)
	rmsamp 		= STDDEV(amprat)
	PLOT, pa, (amprat<0.49)>(-0.49),TITLE=STRING(isc+1, rmsamp, FORMAT="('SC', I2, ' RMS=', F5.3)"), 	$
				XRANGE=[0,360],      XSTYLE=1, XTITLE='POSN_ANGLE(deg)', 							$
				YRANGE=[-0.5, +0.5], YSTYLE=1, YTITLE='(A0-A180) / (A0+A180)',						$
		  		CHARSIZE=1.0, PSYM=6, SYMSIZE=0.5
	OPLOT, [0,360], [0,0]
	phz 		= ATAN(y,x) *!RADEG
	phzz		= SHIFT(phz, nok2)
	phzsum		= ((phz+phzz+540) MOD 360) - 180.		; forces result into +-180 range
	nokps[isc] = nok									; save results for histogram
	phzsumtab[0:nok-1, isc] = phzsum							; save results for histogram
	phzsumav	= MEAN(phzsum)
	phzsumsem	= STDDEV(phzsum) / SQRT(nok2)
	PLOT, pa, (phzsum<59)>(-59), CHARSIZE=1.0, PSYM=6, SYMSIZE=0.5,				$
		TITLE=STRING(isc+1, phzsumav, phzsumsem, FORMAT="('SC', I2, ' <PHZSUM>=', F5.1, ' +-', F4.1, ' deg')"), 	$
				XRANGE=[0,360],    XSTYLE=1, XTITLE='POSN_ANGLE(deg)', 		$
				YRANGE=[-60, +60], YSTYLE=1, YTITLE='PHZ0+PHZ180 (deg)'
	OPLOT, [0,360], [0,0]
ENDFOR
!P.MULTI = [0,3,3]
nbin= (2*maxphz/histphzbin)
binedges = (FINDGEN(nbin)+0.5)*histphzbin -maxphz
FOR isc=0,8 DO BEGIN
	IF nokps[isc] EQ 0 THEN CONTINUE
	PLOT, binedges, HISTOGRAM(phzsumtab[0:(nokps[isc]-1), isc], BINSIZE=histphzbin, MIN=-maxphz, MAX = maxphz)/2, $
			psym=10, TITLE = STRING(isc+1, FORMAT='("SUBCOLL ", I2)'), XTITLE='PHASE SUM (degrees)', YTITLE='NUMBER'
ENDFOR
stop
RETURN
END