; NAME:
;   readnoise.pro (version 1.6)
;
; PURPOSE:
;   This program calculates read noise in a specific region of two
;   images obtained with an IR array. It assumes that the images have
;   already been collapsed to 2 dimensions, i.e. individual CDS images should
;   already be subtracted before using this code.
;
; CALLING SEQUENCE (all procedures):
;   rnval = readnoise(im1, im2, region)
;
; INPUTS (all procedures)
;   im1, im2    - input image arrays (collapsed from CDS or Fowler to 2D plane)
;   region      - region of 2D plane over which readnoise should be calculated.
;                 This should be a 4-vector in the form [x1, x2, y1, y2] to
;                 specify an array subsection A[x1:x2, y1:y2]
;   colinc		- column increment for region
;   rowinc		- row increment for region
;
; KEYWORD PARAMETERS:
;   none
;
; EXAMPLES:
;   Calculate read noise
;   rnval = readnoise(im1, im2, region)
;
; REFERENCE:
;
; MODIFICATION HISTORY:
;   Written by:
;       Don Figer, IDTL
;   Modified by:
;       Michael Telewicz, IDTL, April 3, 2002
;           Added header to script
;       Don Figer, IDTL, April 4, 2002
;           Added realign parameter to choose data shifting in makecds.
;           Note that makecds will probably be moved out of readnoise.pro
;       Michael Telewicz, IDTL, April 12, 2002
;           Added batchreadnoise procedure
;           Added keywords to makecds and readnoise to enhance
;           functionality of batchreadnoise
;       Ernie Morse, IDTL, July, 2002
;           Added pixel rejection functionality
;           Modified to write output to FITS file
;           Modified batchreadnoise so that input file can contain
;           both CDS and Fowler images
;       Ernie Morse, IDTL, Nov. 15, 2002
;           removed sigma clipping and bad pixel mask input parameters
;           changed calculation to place pixels into histogram, fit to Gaussian
;           and get std. dev. from that fit instead of taking a std. dev. after
;           sigma clipping and bad pixel rejection
;       Ernie Morse, IDTL, April 14, 2003 (v1.1)
;           Removed the bin representing 0 from the gaussfit to protect against
;           spurious zeros caused by railed pixels.
;       Ernie Morse, IDTL, April 16, 2003 (v1.2)
;           Took out the zero-bin rejection--it has been replaced by using NaN
;           to mark railed pixels.
;       Ernie Morse, IDTL, February 6, 2004 (v1.3)
;           Edited to ensure compliance with IDTL coding standards.
;       Ernie Morse, IDTL, May 18, 2004 (v1.4)
;           Added NaN flag for histogram function
;       Ernie Morse, IDTL, June 30, 2004 (v1.5)
;           Changed _min and _max to long integer values to prevent problems
;           with integer overflow.
;       Ernie Morse, IDTL, August 9th, 2004 (v1.6)
;           Added rowinc and colinc input variables for region
;       Don Figer, IDTL, October 9th, 2004 (v1.6)
;           Added traps to capture cases when good data are not available. In
; 			these cases, a zero will be returned.
;
function readnoise, im1, im2, region, colinc, rowinc

    ; difference image is target minus bias
    diffarray = im2 - im1

    ; save region of interest into its own array
    roi = diffarray[region[0]:region[1]:colinc, region[2]:region[3]:rowinc]

    ; find standard deviation and median of difference array
    ; stddev doesn't like NaN, so don't consider points with this value
    ; note that we only perform stddev when good data are available
    goodroi=where(finite(roi),ngoodroi)
    if (ngoodroi ne 0) then begin
    	sigma = stddev(roi[goodroi])
    endif else begin
    	sigma=0
    endelse

    med = median(diffarray[region[0]:region[1]:colinc,region[2]:region[3]:rowinc], /even)

    ; use sigma to set histogram limits.  Sigma*5 should be inclusive enough.
    ; Set binsize to 0.2 for now.
;    _size = 1.0
;    _min = long(med - (sigma * 5)) + _size / 2.0
;    _max = long(med + (sigma * 5)) - _size / 2.0
    _min = long(med - (sigma * 5))
    _max = long(med + (sigma * 5))
    _nbins = 128

	if (_max le _min) then _max=_min+1

    math_err_stat = check_math(mask = 208)
    math_errcount = 0
    catch, errval
    if (errval ne 0) then begin
        _nbins = _nbins / 2
        message, /reset
        math_errcount = math_errcount + 1
    endif

    ; make histogram and array of corresponding bin values.
    ; bins are set up to be centered around integer values
    ;rnhist = histogram(roi, min = _min, max = _max, binsize = _size, /nan)
    ;rnvals = findgen(n_elements(rnhist)) * _size + _min + (_size / 2.0)
    rnhist = histogram(roi, min = _min, max = _max, nbins = _nbins, location = rnvals, /nan)

    ; adjust bins so that values are at center of bin
    rnvals = rnvals + ((_max - _min) / (_nbins -1)) / 2.0

    ; filter out bins with no entries.  This is necessary in the event that
    ; the input arrays have not been reference pixel corrected, as all the
    ; pixels will have integer values, and all bins in between integer values
    ; will therefore be empty.
    goodbins = where(rnhist ne 0,ngood)

    ; fit histogram to a Gaussian, saving the coefficients in array coeff
    ; note that we only perform a fit when good data are available
    if (ngood ne 0) then rnfit = gaussfit(rnvals[goodbins], rnhist[goodbins], coeff)
    math_err_stat = check_math(mask = 208)
    if (math_err_stat ne 0 and math_errcount lt 5) then begin
        message, !error_state.msg
    endif

    ; readnoise is the width of the distribution, divided by sqrt(2)
    ; (because we are using the difference of two images)
    ; note that we handle the case in which there was no valid fit by setting
    ; readnoiseval to zero.
    readnoiseval=0.
    if (ngood ne 0) then readnoiseval = coeff[2] / sqrt(2)

    return, readnoiseval

end