	subroutine cosmicra(mode,ic,ifm,id,jd,iframemask,frame,avg,adu2e,icrcnt,ihist,ihit,neighbors)
	implicit none
c
c	This version thought to be sort of good for all modes, all flight cameras, but see below
c
	integer*4 mode,ic,ifm,id,jd,iframemask(id,jd)
	integer*4 i,j,imin,imax,ii,jj,ibig,icrcnt,ihist(250),jarg,ihit(id,jd),neighbors,iflag
c	real*4 avgbin(id,jd)
	real*4 frame(id,jd),avg(id,jd),avgbin(636,128),adu2e
	real*4 crmin,crtest,z(3,3),z0,sumzsq,dev,devmax,arg,amodesq,abin,summax
c
	data crmin,devmax /1000.,5.0/	!1000 e- threshold, 5 sigma	!2x2 binning only
c
	icrcnt=0
	neighbors=0
	amodesq=float(mode)*float(mode)
	crtest=crmin/(amodesq*adu2e)	!convert to desired response threshold for cosmic rays
	abin=amodesq*adu2e/100.		!multiplying response by this yields 100 electron bins
	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
c
	do j=1,jd			!subtract background in preparation for CR removal
	do i=imin-1,imax+1
c		if(iframemask(i,j).eq.1)then
		if(frame(i,j).gt.0.)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.crtest)then
c			print *,'Suspected...',i,j,avgbin(i,j)
		    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
			if(z(ii,jj).gt.z(2,2))go to 99		!Quicker bail out on this pixel
		    enddo
		    enddo
		    if(ibig.eq.0)then
c		    	print *,i,j,'Passed FOV cut & biggest one cut...'
			sumzsq=0.		!calculate standard deviation
			summax=0.
			do ii=1,3
			do jj=1,3
			    if(ii.ne.2.or.jj.ne.2)then
			         arg=z(ii,jj)*z(ii,jj)
			         sumzsq=sumzsq+arg
				 if(arg.gt.summax)summax=arg
			    endif
			enddo
			enddo
			dev=sqrt(sumzsq/5.)	!Note altho 8 bins contribute, 3 free params, 5 deg of freedom...
			iflag=0
			if(summax.gt.sumzsq/2.)then
c				print 1000,i,j,z
1000				format(' i,j,z = ',2i5,9f6.0)
				sumzsq=sumzsq-summax
				dev=sqrt(sumzsq/4.)
				iflag=1
			endif
c			print *, 'Dev = ',dev, ', z(2,2) = ',z(2,2)
c			if(dev.lt.2.)print *,i,j,dev
			if(dev.lt.2.)dev=2.	!Limit unusually small values of std dev...
			arg=z(2,2)/dev
			if(arg.ge.devmax)then	!OK, so it's a 2 x 2 bin cosmic ray...
c				print *,'Got one: ',i,j,avgbin(i,j),'-->',z0
			    jarg=1+nint(abin*z(2,2))
			    if(jarg.gt.250)jarg=250
			    ihist(jarg)=ihist(jarg)+1
			    ihit(i,j)=ihit(i,j)+1
			    icrcnt=icrcnt+1
			    if(iflag.eq.1)neighbors=neighbors+1
			    avgbin(i,j)=z0	!"fix it"
			endif
		    endif
		endif
99	    continue
	    enddo
	enddo
	do j=1,jd				!Replace background after CR removal
	do i=1,id
	frame(i,j)=avgbin(i,j)+avg(i,j)
	enddo
	enddo
	return
	end

