pro welldepth_reduce12, testdir, infile = infile, outdir=outdir, $
        region = region, rowinc = rowinc

	; set directory path delimiter based on
	; operating system in use (non-unix is assumed to be Windows).
	if (!version.os_family eq 'unix') then begin
		pathdelim = '/'
	endif else begin
		pathdelim = '\'
	endelse

	; save the current plot type
	plottype = !d.name

	; check to see if input directory name ends in path delimiter--if not,
	; add it to the end
	if (strmid(testdir, strlen(testdir)-1, 1) ne pathdelim) then begin
	    testdir = testdir + pathdelim
	endif

	; set outdir to testdir + 'Results' if not set by user
	if (keyword_set(outdir) eq 0) then begin
		outdir = testdir + 'Results' + pathdelim
	endif

	; check to see if outdirectory name ends in path delimiter--if not,
	; add it to the end
	if (strmid(outdir, strlen(outdir)-1, 1) ne pathdelim) then begin
		outdir = outdir + pathdelim
	endif

    if (not keyword_set(infile)) then begin
        infile = testdir + 'ramps.txt'
    endif

    if (not keyword_set(rowinc)) then begin
        rowinc = 1
    endif

    ; read the list of ramps
    readcol, infile, ramps, vref, offset, format = '(A, F, Z)'

    ; get information from the header of the 1st image
    ; images are assumed to be 3d image cubes
    fits_read, testdir + ramps[0], 0, h1, /header_only

	xsize = fix(sxpar(h1, 'naxis1'))
	ysize = fix(sxpar(h1, 'naxis2'))
	zsize = fix(sxpar(h1, 'naxis3'))
	scaName = sxpar(h1, 'detname')
	dettemp = strtrim(string(sxpar(h1, 'dettemp'), format = '(F5.1)'), 2)
	wheel1 = sxpar(h1, 'wheel1')
	wheel2 = sxpar(h1, 'wheel2')

    ; find the number of pixels in one frame
    framepix = long(xsize) * long(ysize)

    ; get the exposure times for each read of the ramp
    exptime = (sxpar(h1, 'exp_time') / (zsize - 1)) * findgen(zsize)

	; extract SCA type from SCA name
	scaType = strmid(scaName, 0, strpos(scaName, '-'))

    ; get voltages from the header
    if (strmid(scaType, 0, 1) eq 'S') then begin
        biases = ['VDETCOM', 'VDDUC', 'VDDCL', 'VROWOFF', 'VSSUC', $
                'ISLEWREF', 'VDDOUT', 'VSUB', 'VSSOUT1', 'VP', 'VNROW', $
                'VGGCL', 'VNCOL']
    endif
    if (strmid(scaType, 0, 1) eq 'H') then begin
        biases = ['DSUB', 'VRESET', 'VBIASPWR', 'VBIASGTE', 'DRAIN', $
                'VPULLA7', 'VPULLA15', 'OUTWINA', 'VREF', 'VPULLA23', $
                'VPULLA31', 'VDDA', 'VDD']
    endif

	; if the user didn't set the region, define default regions based on the
	; SCA type
	if (not keyword_set(region)) then begin
		case scaType OF
			'H1RG' : begin
				region = [4, 1019, 4, 1019]
			end
		  	'H2RG' : begin
		  		region = [4, 2043, 4, 2043]
		  	end
			'SB304': begin
				region = [4, 2051, 0, 2047]
			end
			ELSE   : begin
				print, 'Detector type ' + scaType + $
				    ' not supported.  Exiting procedure.'
				return
			end
		endcase
	endif

    ; adjust VREF and offsets to be relative to values for 1st ramp
    vref = vref - vref[0]
    offset = offset - offset[0]

   ; adjust VREF into ADU, assuming 0.1V corresponds to 25K shift in ADU
   vref = vref * 2.5e4 / 0.1

   ; adjust offset into ADU, assuming 1.0V = 7400 ADU (by experiment)
   offset = offset / 4095.0 * 10.0 * 7.4e3

    ; find where images are off both rails. Be conservative
    ramp_medians = fltarr(n_elements(ramps), zsize)

    for rampcount = 0, n_elements(ramps) - 1 do begin
        for readcount = 0, zsize - 1 do begin
            fits_read, testdir + ramps[rampcount], r, $
                    first = readcount * framepix, $
                    last = (readcount + 1) * framepix - 1
            r = reform(r, xsize, ysize, /overwrite)
            ramp_medians[rampcount, readcount] = median(r[region[0]:region[1], $
                    region[2]:region[3]:rowinc])
        endfor
    endfor

    ; find offset between each of the ramps
    ramp_offsets = fltarr(n_elements(ramps) - 1)

    for i = 1, n_elements(ramps) - 1 do begin

        ; find nonrailed reads of 1st ramp
        nonrailed_1 = where(ramp_medians[i - 1, *] gt -31000.0 $
                and ramp_medians[i - 1, *] lt 31000.0)

        ; find nonrailed reads of 2nd ramp
        nonrailed_2 = where(ramp_medians[i, *] gt -31000.0 $
                and ramp_medians[i, *] lt 31000.0)

	    ; use a linear fit to the last 3 good points of the first ramp and
	    ; the first 3 points of the second ramp.
	    p1 = poly_fit(exptime[nonrailed_1[n_elements(nonrailed_1) - 3:*]], $
	            ramp_medians[i - 1, $
	            nonrailed_1[n_elements(nonrailed_1) - 3:*]], 1, /double)
	    p2 = poly_fit(exptime[nonrailed_2[0:2]], $
	            ramp_medians[i, $
	            nonrailed_2[0:2]], 1, /double)

	    ; match to the midpoint between the ends of the two ramps
	    midtime = mean([exptime[nonrailed_1[n_elements(nonrailed_1) - 1]], $
	            exptime[nonrailed_2[0]]])

	    ramp_offsets[i - 1] = poly(midtime, p1) - poly(midtime, p2)

    endfor

    ; find out where saturation occurs in the last ramp
    ; define saturation as the first read whose median signal is within 200 ADU
    ; of the previous read.
    ; start by calculating difference of median of each read in last ramp with
    ; the previous read.
    nonrailed = where(ramp_medians[n_elements(ramps) - 1, *] gt -31000.0)
    med_diff = ramp_medians[n_elements(ramps) - 1, nonrailed[0] + 1:*] $
            - ramp_medians[n_elements(ramps) - 1, nonrailed[0]:zsize - 2]

    ; presatread is the read # of the last read prior to saturation
    presatread = (where(med_diff lt 200.0))[0] + nonrailed[0]

    ; read in the last non-saturated read and the 1st saturated read
	fits_read, testdir + ramps[n_elements(ramps) - 1], topend, $
	        first = presatread * framepix, $
	        last = (presatread + 2) * framepix - 1
    topend = double(topend)
    topend = reform(topend, xsize, ysize, 2, /overwrite)

    ; find saturated pixels in region of interest (roi)
    diff = abs(topend[*,*,1] - topend[*,*,0])
    roi = diff[region[0]:region[1], region[2]:region[3]:rowinc]
    satpix = where(roi lt 200)

    ; get the baseline read of the 1st ramp, which is the first read
    ; that is not railed at the low end of the A/D range
    baseread = (where(ramp_medians[0,*] gt -31000.0))[0]
    fits_read, testdir + ramps[0], bottomend, $
            first = baseread * framepix, $
            last = (baseread + 1) * framepix - 1

    bottomend = double(bottomend)
    bottomend = reform(bottomend, xsize, ysize, /overwrite)

    ; calculate the offset in ADU between the 1st and last ramp and subtract
    ; it from the first ramp.
;    total_offset = vref[n_elements(ramps) - 1] + offset[n_elements(ramps) - 1]
    total_offset = total(ramp_offsets)
    bottomend = bottomend - total_offset

    ; find accumulated signal between 1st and last ramps
    sig_accum = topend[*,*,1] - bottomend
    roi = sig_accum[region[0]:region[1], region[2]:region[3]:rowinc]

    ; make a histogram of the accumulated signal in the saturated pixels in the
    ; region
    binsize = 100
    satpixhist = histogram(roi, binsize = binsize, locations = signal, $
            omax = signal_max, omin = signal_min)

    ; convert from raw pixel counts to percentage of pixels in region
    satpixhist = satpixhist / float(n_elements(roi)) * 100.0

    ; adjust signal bin values to show the center value of each bin
    signal = signal + binsize / 2.0

    ; find cumulative percentages in each bin
    cumulative = fltarr(n_elements(signal))
    cumulative[0] = satpixhist[0]
    for bincount = 1, n_elements(signal) - 1 do begin
        cumulative[bincount] = cumulative[bincount - 1] + satpixhist[bincount]
    endfor

    ; find the exposure times and corresponding values of non-railed medians
    nonrailed = where(ramp_medians gt -31000.0 and ramp_medians lt 31000.0)

    ; convert to 2-D indices
    nonrailed_2d = array_indices(ramp_medians, nonrailed)

    ; get the exposure times for use in fitting
    fit_times = exptime[nonrailed_2d[1,*]]

    ; apply the offset to each of the ramp medians
;    rev_offset = reverse(vref + offset)
    for rampcount = 0, n_elements(ramps) - 2 do begin
;        ramp_medians[rampcount, *] = ramp_medians[rampcount, *] $
;                 - rev_offset[rampcount]
        ramp_medians[rampcount, *] = ramp_medians[rampcount, *] $
                 - total(ramp_offsets[rampcount:*])
    endfor

    fit_signal = ramp_medians[nonrailed]
    fit_signal = fit_signal[1:*] - fit_signal[0]

    frac_rate = [fit_signal / fit_times[1:*]]

    fitvals = where(fit_signal lt max(fit_signal) * 0.8)

    p = poly_fit(fit_signal[fitvals], frac_rate[fitvals], 2, /double)
    p1 = p / p[0]

    ugood = where(fit_signal gt 0.0)

	fname = outdir + scaName + '_linsat.jpg'

    ; transform bias values into a string for plot output
    bias_out = strarr(n_elements(biases))

    for i = 0, n_elements(biases) - 1 do begin
        bias_out[i] = strtrim(string(sxpar(h1, biases[i]), format = '(F7.3)'), 2)
    endfor

    black = 0
    white = 255

	; Open the z-buffer device
	set_plot, 'z'

	; set resolution and font
	device, set_resolution=[8000,6000]
	device, set_font='Courier'

	!P.MULTI = [0, 1, 2] 	; Set up for two plots per page
	fsub = fname  ; set text to print at bottom of plot

	; set plot margin & tick marks
	!Y.MINOR=0
	!y.margin=[6,4]
	!x.margin=[12,12]

	; set some size values
	!p.charsize = 8.0
	!p.charthick = 15.0
	!p.thick = 10.0
	!x.thick = 20.0
	!y.thick = 20.0

    plottitle = scaName + ' Well Depth Histogram at ' + dettemp + ' K (' $
            + wheel1 + ' + ' + wheel2 + ')'

    ; define x and y coordinates of histogram peak
    ymax = max(satpixhist, ymax_sub)
    x_ymax = signal[ymax_sub]

    ; set x-axis range so that peak is in center of plot
    delta_x = (x_ymax - signal_min) > (signal_max - x_ymax)
    xmax = x_ymax + delta_x
    xmin = x_ymax - delta_x
    xwidth = xmax - xmin

    plot, signal, satpixhist, $
            psym = 10, $
            xtitle = 'Well Depth (ADU)', $
            ytitle = '% Number of Pixels', $
            title = plottitle, $
            xrange = [xmin, xmax], $
            yrange = [0, ymax * 1.1], $
            ystyle = 9, $
            xstyle = 1, $
            background = white, $
            color = black

    ; draw vertical line through histogram peak
    plots, [x_ymax, x_ymax], [0, ymax], linestyle = 2, color = black

    ; add annotation of peak value
    xyouts, x_ymax + xwidth * 0.012, ymax * 0.1, $
            strcompress(string(x_ymax), /rem), $
            orient = 90, color = black

    ; add axis for cumulative percentage
    axis, yaxis = 1, yrange = [0, 100], /save, ytitle = 'Cumulative %', $
            ytickv = indgen(11) * 10., $
            yticks = 10, $
            yminor = 5, $
            color = black

    ; plot the cumulative percentage
    oplot, signal, cumulative, lines = 3, color = black

    ; add biases to plot
    for i = 0, n_elements(biases) - 1 do begin
        xyouts, x_ymax + xwidth * 0.3, 85. - (i*5), biases[i]+' = ' $
                + bias_out[i] + ' V', color = black
    endfor

    ; plot linearity
    plot, fit_signal[ugood], frac_rate[ugood]/p[0], yra = [0.6, 1.15], $
        xra = [0, max(fit_signal)*1.1], xstyle = 1, $
        ystyle = 1, title = scaName + ' Linearity', xtitle = 'Total DN in Well', $
        ytitle = 'Fractional Countrate', xmar = [8, 6], $
        background = white, color = black

    circsym, 5, 0
    oplot, fit_signal[ugood], frac_rate[ugood]/p[0], psym = 8, color = black
;    circsym, 0.6, 1
;    oplot, fit_signal[ugood], frac_rate[ugood]/p[0], psym = 8, color = black
    circsym, 5, 1
    oplot, fit_signal[fitvals], frac_rate[fitvals]/p[0], psym = 8, color = black

    oplot, findgen(max(fit_signal) * 1.1), $
            poly(findgen(max(fit_signal) * 1.1), p1), linestyle = 3, $
                    color = black

    plots, [0, max(fit_signal) * 1.1], [1.0, 1.0], linestyle = 2, $
            color = black

    a1 = strcompress(string(p1[1], format = '$(e20.10)'), /rem)
    a1 = strmid(a1, 0, 5)+strmid(a1, strpos(a1, 'e'))
    a2 = strcompress(string(p1(2), format = '$(e20.10)'), /rem)
    a2 = strmid(a2, 0, 5)+strmid(a2, strpos(a2, 'e'))

    xyouts, max(fit_signal) * 0.1, 0.8, 'fit = 1 + (' + a1 + ' * DN) + (' + $
            a2 + ' * DN^2)', color = black
    xyouts, max(fit_signal) * 0.2, 0.77, ' =  included in fit', /data, $
            color = black
    xyouts, max(fit_signal) * 0.2, 0.74, ' =  not included in fit', /data, $
            color = black
    circsym, 5, 1
    oplot, [max(fit_signal) * 0.17], [0.775], psym = 8, color = black
    circsym, 5, 0
    oplot, [max(fit_signal) * 0.17], [0.745], psym = 8, color = black
;    circsym, 0.6, 1
;    oplot, [max(fit_signal) * 0.17], [0.745], psym = 8, color = black


    ; for some single pixels
    ;plot, r2(1000, 1100, *)-float(r2(1000, 1100, 0)), $
    ;    ((r2(1000, 1100, *)-float(r2(1000, 1100, 0)))/(findgen(10)*30.))/ $
    ;    280., yra = [0.6, 1.1]
    ;for i = 1000, 1100 do oplot, r2(i, 1100, *)-float(r2(i, 1100, 0)), $
    ;    ((r2(i, 1100, *)-float(r2(i, 1100, 0)))/(findgen(10)*30.))/280.

    ; write out jpg image
	jpgimg = tvrd()
	write_jpeg, fname, $
	        congrid(jpgimg, 1600, 1200, /center, /interp), quality=100

	; finish up by closing device
	device, /close


    !p.multi = 0


	!p.font = (-1)

end