;+
; Project     :	SOHO - CDS
;
; Name        :	FIND_COSMIC_RAYS
;
; Purpose     :	Find cosmic rays and determine their strengths
;
; Category    :	Class2, Instrument, Engineering, Statistics
;
; Explanation :	This procedure finds all the cosmic rays above a certain
;		threshold value which are in a CDS NIS raster.  Each cosmic ray
;		is counted, and the total number of counts it added to the
;		image over all the pixels it affected is calculated.
;
; Syntax      :	FIND_COSMIC_RAYS, INPUT, STRENGTH
;
; Examples    :	FIND_COSMIC_RAYS, 's7105r00', STRENGTH
;
;		QLDS = READCDSFITS('s7105r00')
;		FIND_COSMIC_RAYS, QLDS, STRENGTH
;
; Inputs      :	INPUT	 = Either the name of CDS FITS file, or a quicklook
;			   data structure from READCDSFITS.
;
; Opt. Inputs :	None.
;
; Outputs     :	STRENGTH = A list of strengths for all cosmic rays found.  If
;			   none are found, then this will contain the single
;			   value of 0.  Statistics can then be performed on
;			   these strength values.
;
; Opt. Outputs:	None.
;
; Keywords    :	THRESHOLD = An optional threshold value to use in searching for
;			    cosmic rays.  The default value is 10.
;
;		ERRMSG = If defined and passed, then any error messages will be
;			 returned to the user in this parameter rather than
;			 depending on the MESSAGE routine in IDL.  If no errors
;			 are encountered, then a null string is returned.  In
;			 order to use this feature, ERRMSG must be defined
;			 first, e.g.
;
;				ERRMSG = ''
;				FIND_COSMIC_RAYS, ERRMSG=ERRMSG, ...
;				IF ERRMSG NE '' THEN ...
;
; Calls       :	READCDSFITS, CDS_CLEAN, GT_WINDATA
;
;		Also calls the internal routines FCR_AN_IMAGE and
;		FCR_REM_NEIGHBORS.
;
; Common      :	None.
;
; Restrictions:	Cannot be used with data taken with either the SUMLINE or
;		SUMWIN compression.  Also, cannot be used with GIS data.
;
; Side effects:	None.
;
; Prev. Hist. :	None.
;
; History     :	Version 1, 06-Mar-1997, William Thompson, GSFC
;
; Contact     :	WTHOMPSON
;-
;
	PRO FCR_REM_NEIGHBORS, MASK, X, Y, IMAGE, STRENGTH, IMAX, JMAX
;
;  Reiteratively remove adjacent pixels from the mask that are part of the same
;  cosmic ray.
;
	FOR XX = (X-1) > 0, (X+1) < IMAX DO BEGIN
	    IF MASK(XX,Y) EQ 1 THEN BEGIN
		STRENGTH = STRENGTH + IMAGE(XX,Y)
		MASK(XX,Y) = 0
		FCR_REM_NEIGHBORS, MASK, XX, Y, IMAGE, STRENGTH, IMAX, JMAX
	    ENDIF
	ENDFOR
;
	FOR YY = (Y-1) > 0, (Y+1) < JMAX DO BEGIN
	    IF MASK(X,YY) EQ 1 THEN BEGIN
		STRENGTH = STRENGTH + IMAGE(X,YY)
		MASK(X,YY) = 0
		FCR_REM_NEIGHBORS, MASK, X, YY, IMAGE, STRENGTH, IMAX, JMAX
	    ENDIF
	ENDFOR
;
	RETURN
	END
;
;==============================================================================
;
	PRO FCR_AN_IMAGE, IMAGE, THRESHOLD, STRENGTH
;
	ON_ERROR, 2
;
;  Generate a mask of all points in the image which are greater than the
;  threshold.
;
	MASK = IMAGE GE THRESHOLD
	NUMBER = 0
	W = WHERE(MASK NE 0, COUNT)
	STRENGTH = 0
;
;  Get the size of the image.
;
	SZ = SIZE(IMAGE)
	IMAX = SZ(1) - 1
	JMAX = SZ(2) - 1
;
;  Step through each cosmic ray, and keep track of the total strengths.
;
	WHILE COUNT NE 0 DO BEGIN
		NUMBER = NUMBER + 1
		X = W(0) MOD SZ(1)
		Y = W(0) / SZ(1)
		MASK(X,Y) = 0
		STR = IMAGE(X,Y)
		FCR_REM_NEIGHBORS, MASK, X, Y, IMAGE, STR, IMAX, JMAX
		STRENGTH = [STRENGTH, STR]
		W = WHERE(MASK NE 0, COUNT)
	ENDWHILE
;
	IF NUMBER GT 0 THEN STRENGTH = STRENGTH(1:*)
;
	RETURN
	END
;
;==============================================================================
;
	PRO FIND_COSMIC_RAYS, INPUT, STRENGTH, THRESHOLD=K_THRESHOLD,	$
		ERRMSG=ERRMSG
;
	ON_ERROR, 2
;
;  Determine the threshold value.
;
	IF N_ELEMENTS(K_THRESHOLD) EQ 1 THEN THRESHOLD = K_THRESHOLD ELSE $
		THRESHOLD = 10
;
;  Initialize STRENGTH to 0.
;
	STRENGTH = 0
;
;  Check the input parameters.
;
	If N_PARAMS() NE 2 THEN BEGIN
	    MESSAGE = 'Syntax:  FIND_COSMIC_RAYS, INPUT, STRENGTH'
	    GOTO, HANDLE_ERROR
	ENDIF
;
;  If a data structure was passed, then use it.  Otherwise, read the FITS file
;  whose name was passed.
;
	IF DATATYPE(INPUT,1) EQ 'Structure' THEN BEGIN
	    A = INPUT
	END ELSE BEGIN
	    IF DATATYPE(INPUT,1) NE 'String' THEN BEGIN
		MESSAGE = 'INPUT must be either a structure or ' +	$
			'a character string'
		GOTO, HANDLE_ERROR
	    ENDIF
;
	    A = READCDSFITS(INPUT)
	    IF DATATYPE(A,1) NE 'Structure' THEN BEGIN
		MESSAGE ='Unable to read FITS file ' + INPUT
		GOTO, HANDLE_ERROR
	    ENDIF
	ENDELSE
;
;  Make sure the data is NIS data.
;
	IF A.HEADER.DETECTOR NE 'NIS' THEN BEGIN
	    MESSAGE = 'Can only be used with NIS data'
	    GOTO, HANDLE_ERROR
	ENDIF
;
;  Check the dimensions of the data.  If summed in either the wavelength or
;  slit direction, then do not attempt to find any cosmic rays.
;
	SZ = SIZE(GT_WINDATA(A,0))
	IF (SZ(1) EQ 1) OR (SZ(3) EQ 1) THEN BEGIN
	    MESSAGE = 'Unable to extract cosmic rays from summed data'
	    GOTO, HANDLE_ERROR
	ENDIF
;
;  Clean the images, so that a difference image can be formed, leaving only
;  cosmic rays.
;
	B = A
	CDS_CLEAN, B, MINSIG=THRESHOLD
;
;  Step through each window, and generate a difference array CC.
;
	N_WINDOWS = N_ELEMENTS(A.DETDESC)
	FOR I_WINDOW = 0,N_WINDOWS-1 DO BEGIN
	    AA = GT_WINDATA(A, I_WINDOW)
	    BB = GT_WINDATA(B, I_WINDOW)
	    CC = AA - BB
	    SZ = SIZE(CC)
;
;  Make sure that the data is in the format (Wavelength, Solar_Y, Exposure).
;
	    IF SZ(0) EQ 4 THEN BEGIN
		CC = REFORM(CC)
		N_EXP = SZ(4)
	    END ELSE BEGIN
		CC = REFORM(REARRANGE(CC,[1,3,2]))
		N_EXP = SZ(2)
	    ENDELSE
;
;  Step through the exposures and extract the CR strengths.
;
	    FOR I_EXP = 0,N_EXP-1 DO BEGIN
		FCR_AN_IMAGE, CC(*,*,I_EXP), THRESHOLD, STR
		IF STR(0) NE 0 THEN BEGIN
		    IF (I_WINDOW EQ 0) AND (I_EXP EQ 0) THEN STRENGTH = STR $
			    ELSE STRENGTH = [STRENGTH, STR]
		ENDIF
	    ENDFOR
	ENDFOR
	GOTO, FINISH
;
;  Error handling point.
;
HANDLE_ERROR:
	IF N_ELEMENTS(ERRMSG) EQ 0 THEN MESSAGE, MESSAGE ELSE	$
		ERRMSG = 'FIND_COSMIC_RAYS: ' + MESSAGE
;
FINISH:
	RETURN
	END