pro egain, indir, outdir, filelist, outps=outps, skipfirst=skipfirst

	; NAME:
	;      pro egain.pro
	;
	; PURPOSE:
	;	This program reduces data for measuring electroning gain
    ;
	; CALLING SEQUENCE:
	;      egain, filelist
	;
	; INPUTS:
    ;      filelist  - The name of a text file containing the images to process
    ;      skipfirst - Skip first n points in fitting
    ;      outps     - Write PostScript
	;
	; KEYWORD PARAMETERS:
	;      outps - Write postscript output
	; EXAMPLE
	;      Run egain to reduce results from an electronic gain experiment. Make PostScript output
	;      egain, filelist, /outps
	;
	; REFERENCE:
	;
	;
	; MODIFICATION HISTORY:
	;       Written by:  B.J. Rauscher, IDTL, October 9, 2002
	;       Modified: D. Figer, IDTL, March 31, 2003
	;				  E. Morse, IDTL, April 3, 2003
	;					Changed x-axis title to VRESET
    ;
    ;				  E. Morse, IDTL, April 7, 2003 (v1.3)
    ;					explicity specified the third dimension value of 0 in the array region definition
    ;					changed the y-axis label to expunge reference to CDS
	;-

	; Set directory string according to handle UNIX and Windows
	if (!version.os_family eq 'unix') then begin
		slashstring='/'
	endif else begin
        slashstring='\'
    endelse

	; set line and text characteristics for publication quality
	!p.charsize = 1.4
    !x.charsize = 1.4
    !y.charsize = 1.4
    !p.charthick = 4.0
    !x.thick = 8.0
    !y.thick = 8.0
    !p.thick = 8.0
    !x.minor = 10
    !y.minor = 5

    ; Declare variables
    FIRST = 1 ; A flag

    ; Check for command line errors
    if (n_params() eq 0) then stop, 'Usage: egain [, indir, outdir, filelist, skipfirst=skipfirst, /outps]'

    ; Set defaults
    if (not keyword_set(indir)) then indir='.'+slashstring	; Default is current directory
    if (not keyword_set(outdir)) then outdir = indir		; Default is indir
    if (not keyword_set(filelist)) then filelist='file.lst'
    if (not keyword_set(outps)) then outps = 0  			; Default does not write PostScript
    if (not keyword_set(skipfirst)) then skipfirst = 0 		; Default fits all points

	; add a slash to the directory names if the last character is not one.
	if (strlen(indir) ne strpos(indir,slashstring,/reverse_search)+1) then indir=indir+slashstring
	if (strlen(outdir) ne strpos(outdir,slashstring,/reverse_search)+1) then outdir=outdir+slashstring

	;save current device name so it can be reset later
    dname = !D.NAME

    ; Read in the file list

    ; Open the file list if it exists
    if (file_test(indir+filelist) eq 0) then begin
    	print,indir+filelist+' does not not exist.'
    	return
    endif
    openr, fileunit, indir+filelist, /get_lun

    ; Calibrate each input file. These are assumed to be CDS
    ; images
    _astring = ''          ; An empty string for file names
    inputfiles = ''        ; An empty string array for file names
    while (not eof(fileunit)) do begin
	    readf, fileunit, _astring
	    inputfiles = [inputfiles, _astring]
    endwhile

    ; Free up fileunit since it is no longer needed
    free_lun, fileunit

    ; Extract data for fitting from each file in the list.

    ; Create data arrays
    vreset = fltarr(n_elements(inputfiles))
    signal = fltarr(n_elements(inputfiles))

    ; Operate on each file in turn. Carefule with the array index
    ; in this loop. Due to the way the file names were read in,
    ; the first element is in position 1.
    for i = 1, n_elements(inputfiles) - 1 do begin

        ; Read in the FITS file
	    fits_read, indir+inputfiles[i], data, header

        ; Get the name, sca_type, date, time, etc. - only need to do this once.
        if (FIRST) then begin
	        FIRST = 0
	        detname = ''
	        getdetname, header, detname
	        sca_type = '' ; Create an empty string to hold the sca_type
	        getscatype, header, sca_type ; Get the sca type
	        thedate = ''
	        getdate, header, thedate
	        thetime = ''
	        gettime, header, thetime
	        print, 'Reducing electronic gain data...'
	        print, 'DETNAME = '+detname
	        print, 'SCA_TYPE = '+sca_type
	        print, 'DATE = '+thedate
	        print, 'TIME = '+thetime
        endif

		; Determine array region to be used for calculating statistics
		if (sca_type eq 'H1RG') then region=data[4:1019,4:1019,0]
		if (sca_type eq 'H2RG') then region=data[4:2043,4:2043,0]
		if (sca_type eq 'SB304') then region=data[4:2051,0:2047,0]

        ; Get the VRESET voltage. This depends on the sca type
        ; Note that for Raytheon devices we do not actually
        ; change VRESET, we change VddUC. For the sake of
        ; simplicity, we will use vreset here.
        if (sca_type eq 'H1RG' or sca_type eq 'H2RG') then begin
	        vreset[i] = sxpar(header, 'VRESET')
        endif else begin
	        if (sca_type eq 'SB304') then begin
	            vreset[i] = sxpar(header, 'VDDUC')
	        endif else begin
	            print,'SCA type not recognized'
	            return
	        endelse
        endelse

        ; Get the mean signal.
        ; The djs_iterstat uses sigma clipping to reject
		; statistical outliers. Here we use sigrej = 5 for 5
		; sigma rejection.
		djs_iterstat, region, sigrej = 5, mean=themean
		signal[i] = themean
	endfor

    ; The data files currently start at index =1. Adjust this so
    ; that they are zero-offset as per IDL spec.
    vreset = vreset[1:n_elements(vreset)-1]
    signal = signal[1:n_elements(signal)-1]

    ; Find the best fitting line. Here we also allow the user to skip points.
    _vreset = fltarr(n_elements(vreset) - skipfirst)
    _signal = fltarr(n_elements(signal) - skipfirst)
    for j = skipfirst, n_elements(vreset) - 1 do begin
        _vreset[j-skipfirst] = vreset[j]
        _signal[j-skipfirst] = signal[j]
    endfor
    theLine = poly_fit(_vreset, _signal, 1, SIGMA = sigma)
    ; Plot the data

    ; Define a special plotting symbol to plot big, fat points
	; define vectors for plotting symbol (circle)
	A = FINDGEN(30) * (!PI*2/29.)		; make a vector of 16 points, A(i) = 2pi/16.
	; set radius for circle.
	radius = 1.5		; define bubble radius
	USERSYM, radius*COS(A), radius*SIN(A),thick=2.,/fill		; define plotting symbol as filled circle

	; select postscript output device if outps=1
    if (outps) then begin
        set_plot, 'ps'
        device, filename=outdir+'egain.ps', /LANDSCAPE,bits_per_pixel=8,/encapsul
    endif else begin
	    ; set display device according to whether using a PC or workstation
	    if (!version.os_family eq 'unix') then begin
			set_plot,'x'
	    endif else begin
	        set_plot,'win'
	    endelse
    endelse

    ; Plot the result
    _xrange = vreset[n_elements(vreset) - 1] - vreset[0]
    _yrange = signal[n_elements(signal) - 1] - signal[0]
    xlimits = [vreset[0] - 0.2 * _xrange, vreset[n_elements(vreset) - 1] + 0.2 * _xrange]
	TITLE='Electronic Gain of SCA '+detname
;	XTITLE='VddUC Voltage (Volts)'
	XTITLE='VRESET Voltage (Volts)'
	YTITLE='Signal in First Read (ADU)'

    ; y-limits depend on sca type (actually, I don't think
    ; so. Uncomment the following line if you have problems.
;        if (sca_type = 'H1RG') then $
    ylimits = [signal[n_elements(signal) - 1] + 0.2 * _yrange, signal[0] - 0.2 * _yrange]

	; stop if SCA type is not recognized
    if (sca_type ne 'H1RG' and sca_type ne 'H2RG' and sca_type ne 'SB304') then begin
		stop, 'SCA type not recognized'
	endif

	; reverse limits in the case that SCA is SB304
    if (sca_type eq 'SB304') then begin
		ylimits=reverse(ylimits)
    endif

	; make plot
	plot, vreset, signal, TITLE=title, XTITLE=xtitle, YTITLE=ytitle, psym=8, XRANGE=[xlimits[0], xlimits[1]], XSTYLE=1, YRANGE=[ylimits[0], ylimits[1]], YSTYLE=1, BACKGROUND=255, COLOR=0,position=[0.2,0.2,0.9,0.9]

    ; Overlay a line fit
    _x0 = vreset[0]
    _y0 = theLine[0] + theLine[1] * _x0
    _x1 = vreset[n_elements(vreset) - 1]
    _y1 = theLine[0] + theLine[1] * _x1
    plots, _x0, _y0, color=0
    plots, _x1, _y1, color=0, /CONTINUE

    ; Calculate gain and error
    elecGain = 1.0e+6 / theLine[1] ; Calculate gain in uV/ADU
    sigma = 1.0e+6 * sigma[1] / theLine[1] / theLine[1]

    ; Add labels
    xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.5+!X.CRANGE(0), (!Y.CRANGE(1)-!Y.CRANGE(0))*0.92+!Y.CRANGE(0), 'Electronic gain = '+strtrim(string(elecGain,format='(F8.3)'),2)+' !4l!XV/ADU', color=0
    xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.5+!X.CRANGE(0), (!Y.CRANGE(1)-!Y.CRANGE(0))*0.88+!Y.CRANGE(0), 'Sigma = '+strtrim(string(sigma,format='(F6.3)'),2)+' !4l!XV/ADU', color=0
    xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.5+!X.CRANGE(0), (!Y.CRANGE(1)-!Y.CRANGE(0))*0.84+!Y.CRANGE(0), 'Skipped first = '+strtrim(string(skipfirst),2)+' points', color=0
    xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.5+!X.CRANGE(0), (!Y.CRANGE(1)-!Y.CRANGE(0))*0.80+!Y.CRANGE(0), 'DATE = '+thedate, color=0
    xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.5+!X.CRANGE(0), (!Y.CRANGE(1)-!Y.CRANGE(0))*0.76+!Y.CRANGE(0), 'TIME = '+thetime, color=0

    ; Close postscript device if using it
	if (outps) then device, /close

	; Write plot into a graphic image file
	set_plot,'z'
	device,set_resolution=[8000,5000]
	DEVICE, SET_FONT='Helvetica'
	radius = 6.		; define bubble radius
	USERSYM, radius*COS(A), radius*SIN(A),thick=2.,/fill		; define plotting symbol as filled circle
	plot, vreset, signal, TITLE=title, XTITLE=xtitle, YTITLE=ytitle, psym=8, XRANGE=[xlimits[0], xlimits[1]], XSTYLE=1, YRANGE=[ylimits[0], ylimits[1]], YSTYLE=1, BACKGROUND=255, $
		COLOR=0,position=[0.2,0.2,0.9,0.9],charsize=10.,charthick=16.,xthick=20.,ythick=20.
    ; Add labels
    xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.5+!X.CRANGE(0), (!Y.CRANGE(1)-!Y.CRANGE(0))*0.92+!Y.CRANGE(0), 'Electronic gain = '+strtrim(string(elecGain,format='(F8.3)'),2)+' !4l!XV/ADU', color=0,charsize=10.,charthick=16.
    xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.5+!X.CRANGE(0), (!Y.CRANGE(1)-!Y.CRANGE(0))*0.88+!Y.CRANGE(0), 'Sigma = '+strtrim(string(sigma,format='(F6.3)'),2)+' !4l!XV/ADU', color=0,charsize=10 ,charthick=16
    xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.5+!X.CRANGE(0), (!Y.CRANGE(1)-!Y.CRANGE(0))*0.84+!Y.CRANGE(0), 'Skipped first = '+strtrim(string(skipfirst),2)+' points', color=0,charsize=10,charthick=16
    xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.5+!X.CRANGE(0), (!Y.CRANGE(1)-!Y.CRANGE(0))*0.80+!Y.CRANGE(0), 'DATE = '+thedate, color=0,charsize=10,charthick=16
    xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.5+!X.CRANGE(0), (!Y.CRANGE(1)-!Y.CRANGE(0))*0.76+!Y.CRANGE(0), 'TIME = '+thetime, color=0,charsize=10,charthick=16
    ; Overlay a line fit
    oplot,[vreset(0),vreset(n_elements(vreset) - 1)],[theLine(0) + theLine(1) * _x0,theLine(0) + theLine(1) * _x1],thick=15.,color=0
	WRITE_JPEG,outdir+'egain.jpg',congrid(tvrd(0),1600,1000,/interp),quality=100
	device,/close

	; restore plotting device that was set before program was run
	set_plot,dname

end