FUNCTION HSI_VISIBILITY_FIT, stack, 	$
	EMPTY_FLAG		= empty_flag,	 	$
    DOUBLE			= double, 			$
	NOERR			= noerr,			$
    gtran_div_2nd	= gtran_div_2nd
;+
; NAME: hsi_visibility_fit
;
; PURPOSE: Function returns an array of visibilities (and associated error estimates) that have
;   been fitted to a NRB x NFB array of stacked count rates, in the form  generated by
;   hsi_phz_stacker. (NRB is the number of roll bins and NFB is the number of phase bins.)
;	Visibilities are calculated for just a single subcollimator and single energy and time intervals.
;
; For each roll bin, it assumes that the count rate dependence on phase is:
;    rate = A*cos(phase(*))+B*sin(phase(*))+C		(eqn 1)
;
;
; METHOD:
;     For each roll bin, calculate the set of count rates using the supplied livetime and counts.
;	  Then, excluding points with no live time, identify cos(phase) and sin(phase) as a pair
;     independent variables and use the IDL REGRESS routine to do the fitting and error estimation.
;
; 	For constant errors, this should be equivalent to E.Schmahl's hsi_vis_fit routine which, in turn has
;		been shown to be mathematically equivalent to G. Hurford's program hsi_remove_sinusoid.pro (Apr 7,'04).
;
;	NB - For future reference, this module could be generalized to fit multiple harmonics simultaneously
;			by increasing the number of terms in eqn 1 and by increasing the dimensionality of the output structure.
;
; INPUTS:
;
;       stack = array of structures produced by the program hsi_phz_stacker.pro
;   	         using the /reform flag (i.e, stack.count is 2-dimensional)
;       	     The reform flag should be the default in hsi_phz_stacker.
;
; KEYWORDS:
;         /DOUBLE   	 Makes calculations in double precision
;         EMPTY_FLAG 	 if set, is the nonzero value to which stack.count has been set if there was no livetime.
;         /GTRAN_DIV_2ND forces division by gridtran to be after the fit. Default is to divide before fitting.
;		  /NOERR		 ignores statistical error in input data.
;
; OUTPUTS:
;   vis = a structure with tags:
;		x: 				FLTARR(nrb)		; x component of visibility for each roll angle (=A)
;		y: 				FLTARR(nrb)		; y component (=B)
;		dc: 			FLTARR(nrb)		; phase-independent component (=C)
;		sigmax:			FLTARR(nrb)		; standard deviations in x,y
;		sigmay:			FLTARR(nrb)
;		errflag:		INTARR(nrb)		; Indicates errors in fitting process.  (0 = good)
;		chisq:			FLTARR(nrb)
;
;		errflag = 0 if successful
;				= 1 if REGRESS indicated a singular array
;				= 2 if REGRESS indicated a small pivot element was used
;				= 3 if there are not enough good phase bins (<3)
;		Other tags should be disregarded unless errflag = 0
;
; HISTORY
; 22-Apr-05 	Initial version, adapted from Ed Schmahl's hsi_vis_fit (ghurford@ssl.berkeley.edu)
; 28-Apr-05 gh	Correct normalization error in DC term in fit.
;-
default, empty_flag,	0 		; Should be set to the same as used by hsi_phz_stacker
default, gtran_div_2nd,	0
;
; Set up output structure for a single subcollimator.
IF SIZE(stack.count, /N_DIMENSIONS) NE 2 THEN MESSAGE,'Stack was not created using reform keyword.'		; This is a fatal problem.
nrb			= N_ELEMENTS(stack[0,*].count) 		; number of roll bins
nfb			= N_ELEMENTS(stack[*,0].count)		; number of phase bins
modeffcorr  = SIN(!PI/nfb) / (!PI/nfb)  		; Correction factor for finite number of phase bins.(=0.99 for nfb=12)
visib		= { x: 			FLTARR(nrb),	$
				y: 			FLTARR(nrb),	$
				dc: 		FLTARR(nrb),	$
				sigmax:		FLTARR(nrb),	$
				sigmay:		FLTARR(nrb),	$
				errflag: 	INTARR(nrb)+3, 	$		; Shows fitting errors (0 = good fit). Preset to missing.
				chisq:		FLTARR(nrb)     }
;
; Begin loop to independently process each roll bin
FOR j=0,nrb-1 DO BEGIN
	ok = WHERE((stack[*,j].livetime) GT 0 AND stack[*,j].count NE empty_flag, nok)
  	IF nok LT 3 THEN CONTINUE				; If there are too few good phase bins, just leave errflag=3 and bail.
;
; Calculate count rates and errors, optionally normalizing by gridtran at this stage
; Input statistical errors are estimated on the basis of the total counts (summed over phase bins) and
;   redistributing them on the basis of the relative live time in each phase bin.
; Reason for doing this is to handle the case where modamp>gridtran which would give abnormally small
;   statistical errors to phase bins with minimal counts.  This way, the regression fit closely matches
;   the technique used in the calculation of grid modulation parameters.
	countsum			= TOTAL(stack[ok,j].count)
	livetimesum			= TOTAL(stack[ok,j].livetime)
	rate 				= stack[ok,j].count / stack[ok,j].livetime
	sigmarate 			= SQRT (countsum / stack[ok,j].livetime / livetimesum)	; =SQRT(count*livetime/sumlt)/livetime
	IF NOT KEYWORD_SET(gtran_div_2nd) THEN BEGIN			; this is the default
		rate 		= rate      / stack[ok,j].gridtran
		sigmarate 	= sigmarate / stack[ok,j].gridtran
	ENDIF
  	cosine 	= COS(stack[ok,j].phase_map_ctr)
	sine	= SIN(stack[ok,j].phase_map_ctr)
	xyin	= FLTARR(2, nok)
	xyin[0,*] = cosine
	xyin[1,*] = sine			; Form a 2 x nok input array
	IF KEYWORD_SET(noerr) THEN	$
			fitparm = REGRESS(xyin, rate, CONST=cparm,									$		; Errors not supplied
				SIGMA=fiterr, CHISQ=chisq, STATUS=status, DOUBLE=KEYWORD_SET(double))	$
	ELSE 	fitparm = REGRESS(xyin, rate, CONST=cparm, MEASURE_ERRORS=sigmarate, 		$		; Use measurement errors
				SIGMA=fiterr, CHISQ=chisq, STATUS=status, DOUBLE=KEYWORD_SET(double))
;
; Apply some normalization factors, including gridtran, if it was not done already
	normfactor = modeffcorr* AVG(stack[ok,j].modamp)		; modulation efficiency
	avgridtran = AVG(stack[ok,j].gridtran)
	IF KEYWORD_SET(gtran_div_2nd) THEN normfactor = normfactor * avgridtran		; this is not the default
;
; load results into output structure for a single subcollimator.
	IF KEYWORD_SET(gtran_div_2nd) THEN visib.dc[j] = cparm / avgridtran ELSE visib.dc[j] = cparm
	visib.x[j] 			= fitparm[0] / normfactor
	visib.y[j]			= fitparm[1] / normfactor
	visib.sigmax[j]		= fiterr[0]  / normfactor
	visib.sigmay[j]		= fiterr[1]  / normfactor
	visib.chisq[j]		= chisq
	visib.errflag[j]	= status
ENDFOR
RETURN, visib
END