pro bias, istart, iend, bias, darks, $
          fdir=fdir,datadir = datadir, fnam = fnam, ext = ext, $
          avgclip = avgclip, rn = rn, $
          noweight = noweight, $
          outdir = outdir, nosave = nosave, noscan = noscan, $
          keck = keck, lick = lick, gemini = gemini, ctio = ctio, ircal = ircal, kermit=kermit,$
		  aobir=aobir,outfile=outfile,auto=auto,instrument=instrument

;+
; PURPOSE:
; quick script for making a bias frame and displaying useful info
; 
; USAGE:
; default is median-scaled median averaging to combine frames
;   (which is not appropropriate for correlated double sampling
;   images, which will have the bias level subtracted: better to use
;   med-averaging *w/o* any weighting (/noweight): see 12/12/97 notes.
;   For CTIO data, procedure automatically sets /noweight)
;
; default filename format is Keck
;
; see MKBIAS.IDL for a wrapper to do multiple bias sets
;
; HISTORY: Written by M. Liu (UCB) 01/29/97 
; 02/24/97 (MCL): converted from BIAS.IDL script to BIAS.PRO
;                 saves 1st image header to final image
;                 uses SXADDHIST add history lines to header
; 05/04/97 (MCL): added /lick and /gemini
; 12/11/97 (MCL): added /avgclip option (must give the readnoise)
; 01/24/99 (MCL): checks if all images have same # of coadds
;                 (useful to ensure input list is correct)
; 06/07/99 (MCL): checks if output directory exists
; 07/07/99 (MCL): added /ircal
; 2003-07-21 MDP: Added /auto
;-


; default file paths
if keyword_set(fdir) then datadir=fdir ; back compatibility
if not(keyword_set(datadir)) then datadir = 'data/'
if not(keyword_set(outdir)) then outdir = 'calib/'
	checkdir,datadir,fallback="./"
	checkdir,outdir

; default names
if not(keyword_set(fnam)) then fnam = 's'
if keyword_set(lick) then fnam = 'd'
if ~(keyword_set(instrument)) then begin
	if keyword_set(ircal) then instrument="IRCAL"
	if keyword_set(kermit) then instrument="Kermit"
ENDIF

; user info
if n_params() lt 2 then begin
    print, 'pro bias,istart,iend,bias,'
    print, '         [datadir="', datadir, '"], [fnam="', fnam, '"], '
    print, '         [avgclip],'
    print, '         [outdir="', outdir, '"], [nosave],'
    print, '         [lick=0],[gemini=0],[keck=1],[ctio],[ircal],[aobir]'
	stop
endif


; for double-differenced datasets, set /noweight 
if keyword_set(ctio) or keyword_set(ircal) then noweight = 1


; starting and ending file numbers


; get images
;loadims, istart[0], iend[0], darks, head0, imlist, ncoadds, $
 ; fd = datadir, fnam = fnam, ext = ext, $
  ;lick = lick, gemini = gemini, ctio = ctio, ircal = ircal,aobir=aobir

fitsloader,transpose([[istart],[iend]]),darks,headers,datadir=datadir,instrument=instrument,/divcoadds,imlist=imlist
  head0 =headers[*,0]

  ncoadds = sxpararr(headers,"NCOADDS")
  exptimes = sxpararr(headers,"EXPTIME")
  
if (n_elements(ncoadds) gt 1)  then $
	if(stddev(ncoadds) ne 0) then begin
    message, 'bias frames have diff number of coadds!!', /info
    print, imlist+'  '+string(ncoadds)+' coadds'
    if getyn('do you wish to abort?',0,noquery=auto) then begin
        message, 'returning to MAIN', /info
		return
    endif
endif

if n_elements(uniqvals(exptimes)) gt 1 then begin
	message, 'bias frames have different exposure times!!', /info
    print, imlist+'  '+string(ncoadds)+' coadds, '
    if getyn('do you wish to abort?',0,noquery=auto) then begin
        message, 'returning to MAIN', /info
		return
    endif
endif


; history string
spawn, 'date', date
spawn, 'pwd', pwd
spawn, 'echo $USER@$HOST', id
sxaddhist, '<<BIAS.PRO: '+date+', '+id+'>>', head0
sxaddhist, '  directory = "'+datadir+'"', head0
;sxaddhist, '  file prefix = "'+fnam+'"', head0
;sxaddhist, '  file suffix = "'+ext+'"', head0
sxaddhist, '  loaded images '+strc(istart[0])+' to '+strc(iend[0]), head0

mdarks = meds(darks)
if keyword_set(kermit) then begin
	; occasionally Kermit produces images that are completely
	; wrong values due to some as-yet-undetermined intermittent bug
	; in the readout. Looks like the CDS final read is being dropped??
	; automatically detect if this has happened and discard those frames.
	
	; we can do this because the bad frames appear to have median values
	; less than -120 or so.
	;
	; cutoff seems to instead by -70 or so in 2004.06.12 data??
	
	cutoff = -70
	plot, ps = 10, /ynoz, mdarks,title="Median Value for Dark Frames per coadd",ytitle="ADU",xtitle="Image Number"
	oplot,!x.crange,[cutoff,cutoff],lines=1
	wshow
	
	wbad = where(mdarks le cutoff,badcount,comp=wgood)
	if badcount gt 0 then begin
		for i=0L,badcount-1 do begin
			print,"Frame "+imlist[wbad[i]]+" appears to be bad. Discarding it."
			sxaddhist,'   skipping bad frame  '+imlist[wbad[i]],head0
		endfor 
		darks0 = darks
		darks = darks[*,*,wgood]
		mdarks= meds(darks)
		stop
	endif	

endif



; average the images
if keyword_set(avgclip) then begin
    if not(keyword_set(rn)) then begin
        rn = 0
        message, '/avgclip option set', /info
        read, 'please enter readnoise per image in electrons = ', rn
    endif
    avgclip, darks, bias, rn, 0, map, /darks, sigrej = sj
    sxaddhist, $
      '  combined with AVGCLIP: /darks, sigrej='+strc(sj), head0
endif else begin
   ; avgmed, bias, darks, /evenmed, noweight = noweight
   ;  AVGMED is horribly slow for large data sets. Use median instead!
   ;  (and besides I don't really understand the weighting in avgmed...)
	bias = median(darks,dim=3,/even)
	sxaddhist,"combined with IDL median(...,/even)",head0
	
;    if keyword_set(noweight) then $
;     sxaddhist, $
;    '  combined with AVGMED - no wts, evenmed averaging', head0 $
; else $
;   sxaddhist, $
;  '  combined with AVGMED - median wts, evenmed averaging', head0
endelse


; print info
print
message, 'bias images = '+strc(istart[0])+' - '+strc(iend[0]), /info
message , 'median bias level = ' + strc(median(bias)), /info
message, 'max - min bias level = ' + strc(max(mdarks)-min(mdarks)), /info
message, 'bias level std dev = ' + strc(stddev(mdarks)), /info
message, '(sorted) bias medians: ' + commalist(mdarks(sort(mdarks))), /info
plot, ps = 10, /ynoz, mdarks,title="Median Value for Dark Frames per coadd",ytitle="ADU",xtitle="Image Number"
sxaddhist, '  max - min bias median = ' + strc(max(mdarks)-min(mdarks)), head0
sxaddhist, '  bias level std dev = ' + strc(stddev(mdarks)), head0
wait,0.5


; save the image
if(!version.os eq 'Win32') then sep_char = '.' else sep_char = ':'   
; if Gemini data, note which channel in the output file name
if keyword_set(gemini) then begin
    fext = '_' + string(sxpar(head0, 'ITIME'), form = '(G0.0)') + 's'
    fext = strlowcase(strc(sxpar(head0, 'CHANNEL'))) + fext 
	outfile = outdir + 'bias' + strc(istart[0]) + fext+".fits"
endif else if keyword_set(lick) then begin
    fext = '_' + string(sxpar(head0, 'EXPOSURE'), form = '(G0.0)') + 's'
	outfile = outdir + 'bias' + strc(istart[0]) + fext+".fits"
endif else if keyword_set(ircal) then begin
; note that loadims puts NREADS into NREADS_O... for NREADS > 1
; jpl: replace with MODE fitscard...


    if (sxpar(head0,'NREADOLD') ne 0) then begin
        fext = '_' + strtrim(string(sxpar(head0, 'COADDTIM')),2) + sep_char $
                   + strtrim(string(sxpar(head0, 'NREADOLD')),2) + sep_char $
                   + strtrim(string(sxpar(head0, 'NRESETS')),2) 
    endif else begin 
        fext = '_' + strtrim(string(sxpar(head0, 'COADDTIM')),2) + sep_char $
                   + strtrim(string(sxpar(head0, 'NREADS')),2) + sep_char $
                   + strtrim(string(sxpar(head0, 'NRESETS')),2) 
    endelse
	outfile = outdir + 'bias' + strc(istart[0]) + fext+".fits"
endif else if keyword_set(kermit) then begin
		modestr  = sxpar(head0,"EXP_MODE",count=count) 
		if (count lt 1) then begin
			modestr = sxpar(headers[*,1],"EXP_MODE",count=count)
		endif 
		if (count lt 1) then begin
			modestr = "COULDNT_FIND_MODESTR"
		endif
		; use - not _ as a delimeter since the modestr can contain _.
	    fext = '-' + strtrim(modestr)
		; store reduced Julian date in the bias file name.
		date = sxpar(head0,"DATE-OBS")
		datev = date_conv(date,'V')
		juldate,datev,rjd
		sxaddpar,head0,"RJD",rjd,"Reduced Julian date of dark frame start"
		noweight=1
	outfile = outdir + 'bias.'+strtrim(modestr)+".rjd"+strc(rjd)+".fits"
endif else if keyword_set(ctio) then begin
    fext = '_' + string(sxpar(head0, 'INT_S'), form = '(G0.0)') + 's'
	outfile = outdir + 'bias' + strc(istart[0]) + fext+".fits"
endif else if instrument eq "NIRC2" then begin
	fext = '_' + strc(string(sxpar(head0,'EXPTIME'),format="(F6.2)"))+sep_char+strc(getnirc2nreads(head0))
	outfile = outdir + 'bias' + string(istart[0],format="(I4.4)") + fext+".fits"
endif else begin

    fext = '_' + string(sxpar(head0, 'EXPTIME'), form = '(G0.0)') + 's'
	outfile = outdir + 'bias' + string(istart[0],format="(I4.4)") + fext+".fits"
endelse


print,"output file name will be "+outfile

; scan if desired
if not(keyword_set(noscan)) then begin
    print
    if getyn('scan difference of indiv image with final bias frame?',0,noquery=auto) then begin
        skysub, darks, bias, diff, /silent
		for k=0,(size(diff))[3]-1 do begin
			print,"Dark "+strc(k)+" has stddev diff "+strc(stddev(diff[*,*,k]))
		endfor 
        scan, diff, wait = 0.5
  ;      med = median(bias)
  ;      sd = stdev(mdarks)
	;	scan_put, diff, min=-3*sd, max=3*sd
    endif
endif

;; check if output file exists -> doesn't work
;if not(exist(outdir)) then begin
;    print, 'output directory "', outdir, '" does not exist!'
;    read, 'enter new output directory: ', outdir
;endif
display2,bias,title="Bias file for "+fext
wshow

if keyword_set(aobir) then begin
	; this is needed to change the FITS header since the
	; input stuff from CFHT has BZERO=32768 but our output here doesn't.
	print,"Updating BZERO keyword in output file"
	sxaddpar,head0,"BZERO",0,"offset"
endif


if ~(keyword_set(auto)) then begin
	atv,[[[bias]],[[darks]]],/bl,names=["bias",imlist]
	stop
endif

message, 'saving image: '+outfile, /info
if not(keyword_set(nosave)) then $
  if filecheck(outfile,noquery=auto,/overwrite) then writefits, outfile, bias, head0
;  writefits, "N"+outfile, bias


end