pro avg_darkdata,input_dir,roistart,darkavg,skipexam=skipexam

;  read, process eclipse season dark data from a single directory input_dir
;
;  INPUTS:
;	input_dir	= input directory for raw dark data
;	roistart	= spectral ROI start pixel
;	skipexam	= if set, then skip manual examination of images
;  OUTPUT:
;	darkavg		= average dark image

;  get the string representation of the start pixel of the ROI
rstrt = strcompress(string(roistart),/remove_all)
;  get the time ID for the filename
ppos = strpos(input_dir,'/20')
plgth = strlen(input_dir)-ppos - 2
timid = strmid(input_dir,ppos+1,plgth)
;  generate the output filename for the selected files for averaging
outfilnam = 'goodarks_' + rstrt + '_' + timid + '.save'

;if outfilnam eq 'goodarks_56_20150528_120037.save' then skipexam = 1

;  bypass manual examination
if keyword_set(skipexam) then begin
	restore,outfilnam
	nfile = n_elements(derf)
	print,'bypassing the exam of each file, read filenames from save file'
	goto,bypass39
endif

startover:


;  First cycle through the files, displaying each one and rejecting bad
;  files
;  find the file names in the current directory
derf = file_search(input_dir + '*.fits')
nfile = n_elements(derf)
print,'directory: ',input_dir,'; number of files = ',nfile
fileindx = intarr(nfile)	;  accept index for each file (1=accept)
goodfile = derf	; sequential list of good file path names
goodfile(*) = ' '

;  loop through files examining them
countd = 0
for kk = 0,nfile-1 do begin
	sp4fits = derf(kk)
;  First read in the data and the FITS header
	print,' Begin processing file: ',kk,' of ',nfile-1,'  ',sp4fits
        hdr=headfits(sp4fits)
        read_sot,sp4fits,index,dat
	dat = float(dat)

;  get some header values
	bitshft = sxpar(hdr,'SPBSHFT')
        nxs2 = sxpar(hdr,'NAXIS2')
        nxs1 = sxpar(hdr,'NAXIS1')

; DIGITAL WRAP, BIT SHIFTING CORRECTIONS
; correct for digital wrap-around of Stokes I, if wrap-around is present
        temp = dat(*,*,*,0)
        whr = where(temp lt 0., countwrap)
        if countwrap ne 0 then temp(whr) = temp(whr) + 65536.
        dat(0,0,0,0) = temp
;  account for bit shifting of the data by multiplying the appropriate
;  Stokes images by 2
        switch bitshft of
                3:      dat(*,*,*,1:2) = 2.*dat(*,*,*,1:2)
                2:      dat(*,*,*,3) = 2.*dat(*,*,*,3)
                1:      dat(*,*,*,0) = 2.*dat(*,*,*,0)
        endswitch

;  initialize dark sums, display area
	if countd eq 0 then begin
		disp = fltarr(8*nxs1,nxs2)
	endif

;  fill display
	for iside = 0,1 do begin

;  get average Stokes I for each side
		for istks = 0,3 do begin
			istrt = iside*nxs1*4 + istks*nxs1
			temp = dat(*,*,iside,istks)
			;if istks eq 0 then temp = temp-min(temp)-100.
			if istks eq 0 then begin
				avgval = mean(temp,/double)
				rms = stddev(temp)
				;temp = temp-avgval(iseq,iside,0)-3.*rms(iseq,iside,0)
				temp = temp-avgval-3.*rms
			endif
			disp(istrt:istrt+nxs1-1,*) = temp
		endfor
	endfor
	strnfile = strcompress(string(nfile-1))
	strkk = strcompress(string(kk))
	titl = 'File ' + strkk + ' of ' + strnfile
	if kk eq 0 then begin
		tvwin,(disp<100.)>(-100.),title=titl
	endif else begin
		tvwin,(disp<100.)>(-100.),title=titl,/currw
	endelse


	accept = ' '
	jump114:
	read,' accept this image for average dark? (y/n), or go back (b) ',accept
;  case to accept the file
	if accept eq 'y' then begin
		fileindx(kk) = 1
;  case to back up sequence
	endif else if accept eq 'b' then begin
                read,' enter number of files to go back: ',bkup
                kk = kk - fix(bkup) - 1
;  be sure not to back up beyond beginning of file
				if kk lt 0 then kk = 0
;  case of rejecting current one
	endif else if accept eq 'n' then begin
		fileindx(kk) = 0
	endif else begin
		goto, jump114
	endelse

;  end loop over selected files
endfor

;  readjust the list of filenames for good files to be averaged
countd = 0
for kk = 0,nfile-1 do begin
	flnm = derf(kk)
	if fileindx(kk) eq 1 then begin
		goodfile(countd) = derf(kk)
		countd = countd + 1
	endif
endfor
goodfile = goodfile(0:countd-1)
derf = goodfile

print,'number of orig files= ',nfile,', no of selected files= ',n_elements(derf)
;save,filename=outfilnam,derf

goto,kilroy1
;  skip manual examination
bypass39:countd = n_elements(derf)
skipexam=0



kilroy1:
;  Now, second pass through the selected files to get average dark images
;  loop through files
nfile = countd
for kk = 0,nfile-1 do begin
	sp4fits = derf(kk)
;  First read in the data and the FITS header
	print,' Begin processing file: ',sp4fits
        hdr=headfits(sp4fits)
        read_sot,sp4fits,index,dat
	dat = float(dat)

;  get some header values
	bitshft = sxpar(hdr,'SPBSHFT')
        nxs2 = sxpar(hdr,'NAXIS2')
        nxs1 = sxpar(hdr,'NAXIS1')

; DIGITAL WRAP, BIT SHIFTING CORRECTIONS
; correct for digital wrap-around of Stokes I, if wrap-around is present
        temp = dat(*,*,*,0)
        whr = where(temp lt 0., countwrap)
        if countwrap ne 0 then temp(whr) = temp(whr) + 65536.
        dat(0,0,0,0) = temp
;  account for bit shifting of the data by multiplying the appropriate
;  Stokes images by 2
        switch bitshft of
                3:      dat(*,*,*,1:2) = 2.*dat(*,*,*,1:2)
                2:      dat(*,*,*,3) = 2.*dat(*,*,*,3)
                1:      dat(*,*,*,0) = 2.*dat(*,*,*,0)
        endswitch

;  initialize dark sums, display area
	if kk eq 0 then begin
		darkavg = fltarr(nxs1,nxs2,2)
		darkavg(*,*,*) = 0.
		meandark = fltarr(nfile,2)
	endif

;  fill display
	for iside = 0,1 do begin
		darkavg(*,*,iside) = darkavg(*,*,iside) + dat(*,*,iside,0)
;  get average Stokes I for each side
		meandark(kk,iside) = mean(dat(*,*,iside,0),/double)
	endfor

;  end loop over files
endfor


;  renormalize dark average
darkavg = darkavg/float(nfile)



;  plot average dark level for CCDSIDE0
meandark0 = median(meandark(*,0))
stdevark0 = stddev(meandark(*,0)-meandark0)
plot,(meandark(*,0)-meandark0),title='CCDSIDE0, less average, DN'
accept = ''
accept2 = ''
print, 'Should be numbers of a few.  type any character to continue: '
jump224:
read, 'To eliminate more data points, type e. Type g to continue: ',accept
if accept eq 'e' then begin
	elim1 = 0 & elim2 = 0
	kilroy2: read,'enter range of points to eliminate: ',elim1,elim2
	print,'eliminate files ',elim1,elim2,' from present sequence?'
	read,'enter either y or n: ',accept2
	if accept2 eq 'n' then goto,kilroy2

;  readjust the list of filenames for good files to be averaged
	countd = 0
	nfile = n_elements(derf)
	for kk = 0,nfile-1 do begin
		flnm = derf(kk)
		if kk lt elim1 or kk gt elim2 then begin
			goodfile(countd) = derf(kk)
			countd = countd + 1
		endif
	endfor
	goodfile = goodfile(0:countd-1)
	derf = goodfile
;	save,filename=outfilnam,derf
	goto, kilroy1
endif
if accept ne 'g' then goto, jump224


		

;  plot average dark level for CCDSIDE0
meandark1 = median(meandark(*,1))
stdevark1 = stddev(meandark(*,1)-meandark1)
plot,(meandark(*,1)-meandark1),title='CCDSIDE1, less average, DN'
print, 'Should be numbers of a few.  type any character to continue: '
read, ' type b to redo this directory, type any other character to continue: ',accept
if accept eq 'b' then goto,startover

;  save good file names, average dark data
save,filename=outfilnam,derf,darkavg

return
end