; NAME:
;			cosmic.pro (v1.4)
;
; PURPOSE:
;			Identifies cosmic ray hits in an input frame that is the difference of two Fowler-16 processed
;			frames (for IDTL data) or a processed CDS frame (REAG).
; CALLING SEQUENCE:
;			cosmic, pic, counts, cross, diag, mcount, mask, lowcut, highcut
;
; INPUTS:
;			pic		input difference frame
;			counts	empty on input, holds counts of pixels in pic that are identified as cosmic ray hits on output
;			cross	empty on input, on output it holds the values of the nearest-neighbor non-diagonal pixels for each
;					pixel in pic identified as cosmic ray hit
;			diag	empty on input, on output it holds the values of the pixels that are located diagonally to each
;				pixel hit by a cosmic ray
;			mcount	empty on input, on output it contains the index values of all the pixels in pic identified as
;				cosmic ray hits
;			mask	on input, it holds an initial cosmic ray mask--on output, it is the final cosmic ray mask
;				(1 = hit, 0 = non-hit)
;			lowcut	lower limit (in ADU) of cosmic ray intensity to consider in input difference frame
;			highcut	upper limit (in ADU) of cosmic ray intensity to consider in input difference frame
;
; KEYWORD PARAMETERS:
;			none
;
;
; EXAMPLE:
;			cosmic, pic, counts, cross, diag, mcount, mask, lowcut, highcut
;
; REFERENCE:
;
;
; MODIFICATION HISTORY:
;       Written by:  Don Hall's UH group
;
;		Ernie Morse, IDTL, May 23, 2003 (v1.1)
;		Added mask input value, eliminated upper count limit (crit2) for selecting cosmic rays,
;		added rejection of pixels around borders of image, added condition that final output mask is
;		only the pixels that are found in common in the input mask and that satisfy the criteria of this procedure
;
;		Ernie Morse, IDTL, June 3, 2003 (v1.2)
;		Added header and additional comments
;
;		Ernie Morse, IDTL, June 5, 2003 (v1.3)
;		Added lowcut and highcut criteria as input parameters--they were formerly hardwired into procedure
;
;		Ernie Morse, IDTL, August 12, 2003 (v1.4)
;		Modified to deal gracefully with a frame that has no cosmic ray hits
;

pro cosmic,pic,counts,cross,diag,mcount, mask, lowcut, highcut

	; get dimensions of input image
	dimpic = size(pic)
	xsize = dimpic[1]
	ysize = dimpic[2]

	;
	; pic is a dark subracted or differenced image
	;

	crossrat=6   ; maximum acceptable ratio of cross pixels
	;
	; first find high pixels (> crit < crit2) and store in vector counts
	;

	mcount = where(pic ge lowcut and pic le highcut)

	if (mcount[0] ne -1) then begin

	    ;add 1 to the mask in the locations of the hot pixels
	    mask[mcount] = mask[mcount] + 1

	    ; only consider the pixels that were both in the input mask and had a value within the critical range
	    ; in the difference frame
	    mcount = where(mask eq 2)


	    if (mcount[0] ne -1) then begin
	        counts  = pic(mcount)
	        print,'initial number of high points ',n_elements(mcount)

	        ; now find neighbor pixel values relative to center and store in array cross

	        cross = fltarr(n_elements(mcount),4)
	        cross(0,0) = pic(mcount+1)/counts
	        cross(0,1) = pic(mcount-1)/counts
	        cross(0,2) = pic(mcount+xsize)/counts
	        cross(0,3) = pic(mcount-xsize)/counts

	        ; now find maximum and minimum values of cross pixels

	        crmin = cross(*,0)<cross(*,1)<cross(*,2)<cross(*,3)
	        crmax = cross(*,0)>cross(*,1)>cross(*,2)>cross(*,3)

	        ; restrict to cases where cross pixels vary muliplicatively by < crossrat

	        m=where(crmax le crossrat*crmin)
	        if (m[0] ne -1) then begin
	            counts=counts(m)
	            cross=cross(m,*)
	            mcount=mcount(m)

	            ; now form diagonal pixel array

	            diag=cross-cross
	            diag(0,0) = pic(mcount+xsize+1)/counts
	            diag(0,1) = pic(mcount-xsize+1)/counts
	            diag(0,2) = pic(mcount+xsize-1)/counts
	            diag(0,3) = pic(mcount-xsize-1)/counts

	            ; set the mask to 1 in the location of the accepted events
	            mask[*] = 0
	            mask[mcount] = 1
	            print,'final number of points',n_elements(counts)
	        endif else begin
	        	; make a NaN array if none of the hit pixels satisy the uniformity criterion
	            cross = make_array(1,4, /float, value=!values.f_nan)
	            diag = make_array(1,4, /float, value=!values.f_nan)
	        endelse
	    endif else begin
	    	; make a NaN array if no pixels were in both the input mask and the critical range
	        cross = make_array(1,4, /float, value=!values.f_nan)
	        diag = make_array(1,4, /float, value=!values.f_nan)
	    endelse
	endif else begin
		; make a NaN array if no pixels were identified as hits in the input mask
	    cross = make_array(1,4, /float, value=!values.f_nan)
	    diag = make_array(1,4, /float, value=!values.f_nan)
	endelse
end