pro findgainf,flat1,flat2,bias1,bias2, $
	gain,readnoise,meandiff,varflat,varbias, $
	fcoadd=fcoadd,bcoadd=bcoadd, $
	xx=xx,yy=yy, $
	median=median,silent=silent

; finds gain and readnoise for a camera given 
; 2 flat fields (or sky frames) and 2 bias frames
; (no need to bias subtract from the flats)
;
; able to handle flats and biases with diff # of coadds
;   (though each pair of flat must have same # of coadds,
;   likewise with the bias frames)
;
; INPUTS
;	flat1, flat2:   flatfields
;	bias1, bias2:   bias frames
;
; OUTPUTS
; 	[gain] = e-/ADU
; 	[readnoise] = e- (per coadd)
; 	meandiff:  calculated in the process (per coadd)
; 	varflat:   variance between diff of 2 flat fields (per coadd)
; 	varbias:   variance between diff of 2 bias frames (per coadd)
;
; OPTIONAL INPUTS
;	fcoadd:  # of coadds for flats
;	bcoadd:	 # of coadds for biases
;	xx: 	 2 element array for x-dimensions
; 	yy: 	 "		   " y-dimensions
; 	/median: use median instead of mean (recommended)
;
; ALGORITHM
;   given flats with Cf coadds, biases with Cb coadds,
;
;   	varflat  = sigma(F1-F2)^2.0/Cf
;   	varbias  = sigma(B1-B2)^2.0/Cb
;
;	meandiff = mean(F1 + F2)/Cf - mean(B1 + B2)/Cb
;	gain = meandiff / (varflat - varbias)
;	readnoise = gain * (varbias/2.)^0.5	
;
;
; 2/20/95 MCL
;
; Future Fixes:
;   - when find mean, probably should use iterstat
;


if n_params() lt 4 then begin
	print,'pro findgainf,flat1,flat2,bias1,bias2,' 
	print,'             gain,readnoise,meandiff,varflat,varbias,'
	print,'             [fcoadd=],[bcoadd=],[xx=],[yy=],[median],[silent]'
	retall
endif


; defaults: regions (whole chip) and coadds (1)
sz = size(flat1)
if keyword_set(xx) eq 0 then xx=[0,sz(1)-1]
if keyword_set(yy) eq 0 then yy=[0,sz(2)-1]
if keyword_set(fcoadd) eq 0 then fcoadd = 1.0
if keyword_set(bcoadd) eq 0 then bcoadd = 1.0


; get diff of flat frames, find variance per coadd
; (use diff b/c removes structure from images)
flatdif = flat1(xx(0):xx(1),yy(0):yy(1)) - flat2(xx(0):xx(1),yy(0):yy(1))
iterstat,flatdif,out,/silent
varflat = out(3)^2.0/fcoadd


; get diff of bias frames, find variance per coadd 
biasdif = bias1(xx(0):xx(1),yy(0):yy(1)) - bias2(xx(0):xx(1),yy(0):yy(1))
iterstat,biasdif,out,/silent
varbias = out(3)^2.0/bcoadd


; calculate the difference of the variances in ADU
vardiff = varflat - varbias


; calculate the mean of differences
if keyword_set(median) eq 0 then begin
   if keyword_set(silent) eq 0 then print,'* using mean values *'
   meandiff = (mean(flat1+flat2)/fcoadd - mean(bias1+bias2)/bcoadd) 
endif else begin
   if keyword_set(silent) eq 0 then print,'* using median values *'
   meandiff = (median(flat1+flat2)/fcoadd - median(bias1+bias2)/bcoadd) 
endelse


; compute gain and readnoise
gain = meandiff/vardiff
readnoise = gain * (varbias/2.)^0.5


; print out results
if keyword_set(silent) eq 0 then begin
print,'chip region (',strc(xx(0)),':',strc(xx(1)),',',strc(yy(0)),':',strc(yy(1)),')'
print,'   gain = ',strc(gain),' e-/ADU'
print,'   readnoise = ',strc(readnoise),' e-'
endif

end