C+
C NAME:
C	cosmicra
C PURPOSE:
C	To remove cosmic-ray hits from SMEI data.
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call cosmicra(mode,ic,ifm,id,jd,iframemask,frame,avg,icrcnt)
C INPUTS:
C	mode		   integer	1 - "Engineering" frame records, no onboard binning
C					2 - "Hi-res" frame records, binned 2 x 2 pixels
C					4 - "Normal Data" frame records, binned 4 x 4 pixels
C	ic		   integer	Camera 0,1,2, or 3
C	ifm		   integer	frame mask mode - 0 mask not imposed, 1 frame mask imposed
C	id		   integer	data frame max x-value (along rows)
C	jd		   integer	data frame max y-value (down columns)
C	iframemask(id,jd)  integer	data mask - 0 data are invalid 1 data are valid
C	frame(id,jd)   	   real		data frame that needs to have cosmic rays removed
C	avg(id,jd)	   real		Expected bland background for this frame *** see below
C OUTPUTS:
C	frame(id,jd)   	   real		data frame that now has cosmic rays removed
C	icrcnt		   integer	number of cosmic rays found & fixed in this frame
C CALLS:
C	plane
C PROCEDURE:
C	Uses the procedure described in the UCSD memo dated 23 January 2002 to remove cosmic rays
C	This version applies only the less stringent "threshold + std dev" method developed for 4 x 4 bins
C	However, it here applies this method for all modes
C	Each frame subtracts a bland background map (avg(id,jd)) 
C	***This background assumption needs to be revisited, as do the threshold/std dev values.
C	These latter may be set differently for each mode
C	Rapid variation of zodiacal light over the course of each orbit may require modification
C	of the background map for each frame as it is brought in -- yuck!
C MODIFICATION HISTORY:
C	Jan, 2003 A. Buffington (UCSD)
C-
	subroutine cosmicra(mode,ic,ifm,id,jd,iframemask,frame,avg,icrcnt)
c
	integer*4 i,j,imin,imax,ii,jj,ibig,icrcnt
	real*4 frame(id,jd),avg(id,jd),avgbin(id,jd),z(3,3),z0,sumzsq,dev,devmax,arg
c
	data threshold,devmax /1111.,20./	
	icrcnt=0
	if(mode.eq.1)then
		imin=24
		imax=1262
	endif
	if(mode.eq.2)then
		imin=12
		imax=631
	endif
	if(mode.eq.4)then
		imin=7
		imax=315
	endif
	do j=1,jd			!subtract background
	do i=imin-1,imax+1
		if(iframemask(i,j).eq.1)then
		    avgbin(i,j)=frame(i,j)-avg(i,j)
		else
		    avgbin(i,j)=0.
		endif
	enddo
	enddo
	do j=2,jd-1			!test for CR's: > threshold, >neighbors, > devmax...
	    do i=imin,imax
		if(avgbin(i,j).gt.threshold)then
		    do ii=1,3
		    do jj=1,3
		    	arg=avgbin(ii-2+i,jj-2+j)
		    	if(arg.eq.0.)go to 99			!Bail out if any of the 9 out of FOV
			z(ii,jj)=arg
		    enddo
		    enddo
		    call plane(z,z0)
		    ibig=0
		    do ii=1,3
		    do jj=1,3
			if(z(ii,jj).gt.z(2,2))ibig=ibig+1
c			if(z(ii,jj).gt.z(2,2))go to 99		!Quicker bail out on this pixel
		    enddo
		    enddo
		    if(ibig.eq.0)then
			sumzsq=0.		!calculate standard deviation
			do ii=1,3
			do jj=1,3
			    if(ii.ne.2.or.jj.ne.2)sumzsq=sumzsq+z(ii,jj)*z(ii,jj)
			enddo
			enddo
			dev=sqrt(sumzsq/5.)	!Note altho 8 bins contribute, 3 free params, 5 deg of freedom...
			if(dev.lt.20.)dev=20.	!Limit unusually small values of std dev, this rarely happens...
			arg=z(2,2)/dev
			if(arg.ge.devmax)then	!OK, so it's a 4 x 4 bin cosmic ray...
			    icrcnt=icrcnt+1
			    avgbin(i,j)=z0	!"fix it"
			endif
		    endif
		endif
99	    continue
	    enddo
	enddo
	return
	end