pro combine_flats_roi56,roi56list,flatz,savnam

;  Routine to read in the map-average flat images and average them over all operations
;  This is done for ROI56 data only

;  INPUTS:
;   roi56list = array of path names to the IDL-save files of lists of good dark
;		images for ROI=56 dark data, one path for each operation considered
;	flatz = array of pathnames to save files of good flat path names.  Used to get the number
;		of good files averaged in each sequence.
;	savnam = string specifying the date of the data, i.e. '_july2011_roi56'


;  !!!MUST BE RUN IN SOLARSOFT IDL

;  first read in all the data and store it for comparison among maps
nmap = n_elements(roi56list)

;  first get the numbers of good files used in the averages
nfile = fltarr(nmap)
for kk = 0,nmap-1 do begin
	restore,flatz(kk)
	nfile(kk) = float(n_elements(derf))
	nftot = total(nfile)
endfor

;  read in the variables
for kk = 0,nmap-1 do begin
	restore,roi56list(kk)
	if kk eq 0 then begin
		sz = size(mflat)
		nxx = sz(1) & nyy = sz(2)
		gavprf1 = fltarr(nxx,nmap)
		gavprf2 = fltarr(nxx,nmap)
		gclrsp1 = fltarr(nxx,nyy,nmap)
		gclrsp2 = fltarr(nxx,nyy,nmap)
		gdelsh = fltarr(nmap)
		gfltav = fltarr(nxx,nyy,2,nmap)
		gmflat = fltarr(nxx,nyy,2,nmap)
		gshfit1 = fltarr(nyy,nmap)
		gshfit2 = fltarr(nyy,nmap)
		gshft1 = fltarr(nyy,nmap)
		gshft2 = fltarr(nyy,nmap)
		gshftsl = fltarr(nxx,2,nmap)
		gspecrfit = fltarr(nyy,nmap)
		gspecshft = fltarr(nyy,nmap)
		gsvrr = fltarr(nyy,2,nmap)
	endif
	gavprf1(*,kk) = avprf1
	gavprf2(*,kk) = avprf2
	gclrsp1(*,*,kk) = clrsp1
	gclrsp2(*,*,kk) = clrsp2
	gdelsh(kk) = delsh
	gfltav(*,*,*,kk) = fltav
	gmflat(*,*,*,kk) = mflat
	gshfit1(*,kk) = shfit1
	gshfit2(*,kk) = shfit2
	gshft1(*,kk) = shft1
	gshft2(*,kk) = shft2
	gshftsl(*,*,kk) = shftsl
	gspecrfit(*,kk) = specrfit
	gspecshft(*,kk) = specshft
	gsvrr(*,*,kk) = svrr

endfor

;  use original quantites for summing
	avprf1(*) = 0.
	avprf2(*) = 0.
	clrsp1(*,*) = 0.
	clrsp2(*,*) = 0.
	delsh = 0.
	fltav(*,*,*) = 0.
	mflat(*,*,*) = 0.
	shfit1(*) = 0.
	shfit2(*) = 0.
	shft1(*) = 0.
	shft2(*) = 0.
	shftsl(*,*) = 0.
	specrfit(*) = 0.
	specshft(*) = 0.
	svrr(*,*) = 0.

;  average the relevant quantities weighted by number of samples
for kk = 0,nmap-1 do begin
	avprf1 = avprf1 + nfile(kk)*gavprf1(*,kk)
	avprf2 = avprf2 + nfile(kk)*gavprf2(*,kk)
	clrsp1 = clrsp1 + nfile(kk)*gclrsp1(*,*,kk)
	clrsp2 = clrsp2 + nfile(kk)*gclrsp2(*,*,kk)
	delsh = delsh + nfile(kk)*gdelsh(kk)
	fltav = fltav + nfile(kk)*gfltav(*,*,*,kk)
	mflat = mflat + nfile(kk)*gmflat(*,*,*,kk)
	shfit1 = shfit1 + nfile(kk)*gshfit1(*,kk)
	shfit2 = shfit2 + nfile(kk)*gshfit2(*,kk)
	shft1 = shft1 +nfile(kk)*gshft1(*,kk)
	shft2 = shft2 + nfile(kk)*gshft2(*,kk)
	shftsl = shftsl + nfile(kk)*gshftsl(*,*,kk)
	specrfit = specrfit + nfile(kk)*gspecrfit(*,kk)
	specshft = specshft + nfile(kk)*gspecshft(*,kk)
	svrr = svrr +nfile(kk)*gsvrr(*,*,kk)
endfor

;  renormalize
	avprf1 = avprf1/nftot
	avprf2 = avprf2/nftot 
	clrsp1 = clrsp1/nftot
	clrsp2 = clrsp2/nftot
	delsh = delsh /nftot
	fltav = fltav/nftot
	mflat = mflat/nftot
	shfit1 = shfit1/nftot
	shfit2 = shfit2/nftot
	shft1 = shft1/nftot
	shft2 = shft2/nftot
	shftsl = shftsl/nftot
	specrfit = specrfit/nftot
	specshft = specshft/nftot
	svrr = svrr/nftot

;  get standard deviations in terms of fractional values of parameters
print,nftot
save,filename='clrsp1_2010avg.save',clrsp1


mdn = median(mflat(*,*,1))
tvwin,((mflat(*,*,1)-mdn)>(-.02*mdn))<.02*mdn

savefile = strcompress('gaintable' + savnam + '.save',/remove_all)

save,filename=savefile,mflat,shftsl, $
	delsh,specshft,specrfit,avprf1,shft1,shfit1,clrsp1, $
	avprf2,shft2,shfit2,clrsp2,fltav,svrr

;  save the needed data in ssw geny format
savegenx,file='gaintable'+savnam+'.geny',mflat,shftsl,specrfit,specshft,svrr

return
end