c
c name:		CRX_ped_darkE2a.for
c purpose:	determines/removes pedestal & dark current, finds black/white measles & cosmic rays
c inputs:	*.nic data files listed in "namelistCRX.txt"
c outputs:	before and after *.grd files + binary files for each frame, the latter --> Aaron's program
c    "		Note here that the pedestal, dark, and CR information incorporated into the "trailer"
c    "		The output file xxx.dat presents graphs of the various quantities
c calls:	subroutines cosmicray, measles (+ plane), and INDEXX 
c mods:		July 2003, Andrew Buffington (UCSD; (858)-534-6630; abuffington@ucsd.edu)
c-
	program main	!program ped_dark
	implicit none
	integer*1 trail1(382),tdum1(12),trail2(21),tdum2(28),trail3(69)
	integer*2 iframe(318,64),head(3)
	integer*4 h,hh,m,mm,s,ss,i,j,k,kmax,kmin,iarg,i1,i2,ineigh,pedrecover,darkrecover,trash,ipcnt(318,64)
	integer*4 idarkcnt,ibaddark,ipedcnt,ibadped,ihist(250),icrcnt,meascnt,ikill(318,64)
	integer*4 icrneig,n,inegmeas,ipixsum,ipixdif
	integer*4 indx(20352)
	integer*4 mode,id,jd
	real*4 dark,darklast,darksum,darkinit,darkinit2,ped,pedlast,pedinit,pedinit2,pedsum,darkslope,peddiff
	real*4 frame(318,64),framein(318,64),frameout(318,64),badsum,pattern(318,64),fleast(318,64),arg,arg1,patcr
	real*4 minresponse(1524),mingrit(1524),t,t0,pixsum,pixdif,cut
	character*37 infile
	character*40 out3file,out2file,out1file  
	character*15 nname
	character*10 hhead(5)
	character*8 ddate(3),dddate
	
	data infile   /'G:/SEAN/Nic/c2frm_2003_day_hhmmss.nic'/		
	data out1file /'G:/SEAN/Crx/c2frm_2003_day_hhmmsscr1.grd'/	!cosmic-ray not removed file...
	data out2file /'G:/SEAN/Crx/c2frm_2003_day_hhmmsscrx.grd'/	!cosmic-ray removed file...
	data out3file /'G:/SEAN/Bin/c2frm_2003_day_hhmmsscrx.bin'/	!cosmic-ray removed binary file...
	
	data ihist /250*0/
	data i1,i2 /228,301/
	data ddate /'05/29/03','05/30/03','05/31/03'/
	data badsum /0./
	data pattern /20352*0./
	data ipcnt /20352*0/
	data n /20352/
	data fleast /20352*100000./
	data trash /0/
	data patcr /6.94/	!6.94 adu's = 500 electrons
	data minresponse /1524*100000./
	data mingrit /1524*100000./

c
	mode = 4	
	id = 1272
	jd = 256
	if(mode.eq.4)then
		id=id/4
		jd=jd/4
	endif
c
c	Bring in minimum response table
c
	print *, 'reading minimum response table in...'
	open (17,file='min2_4.dat',readonly)	
	open (18,file='mingrit2_4.dat',readonly)
	open (19,file='aurora.txt')
		
	do i=1,1524
	read(17,17)i,minresponse(i)
	read(18,17)i,mingrit(i)
17	format(i6,f10.2)
	enddo
	close(17)	
	close(18)	
	
	open (10,file='namelist2CRX.txt')
	open (12,file='CRX_E2(5-29-03).dat')

	pedrecover=0
	pedinit=1013.
	pedinit2=1017.
	darkinit=4.5
	darkinit2=8.5
	darkrecover=0
	icrcnt=0
c********************************************************************
	kmin=1	
	kmax=26609	!18262
	
	do k=1,kmax
	if(k.eq.11)then
		pedlast=pedinit
		darklast=darkinit
	endif

	read (10,10)nname,h,hh,m,mm,s,ss
10	format(39x,a15,6i1)
	t=3600.*float(10*h+hh)+60.*float(10*m+mm)+float(10*s+ss)
	if(k.gt.18262)t=t+86400.
	if(k.eq.11)t0=t	
	if(k.gt.10.and.k.lt.kmin)go to 98
	if(iabs(k-8).le.1)go to 98		!three open-shutter frames, kill...
c***	
	write(infile(13:27),'(a15)')nname
	write(infile(28:33),'(6i1)')h,hh,m,mm,s,ss
c	print *,infile
	open(11,file=infile,form='binary',readonly)
	read(11)head,((iframe(i,j),i=1,318),j=1,64),trail1,tdum1,trail2,tdum2,trail3
	close(11)
	ped=0.
	dark=0.
	idarkcnt=0
	darksum=0.
	ipedcnt=0
	pedsum=0.
	if(k.gt.10)icrcnt=0
	icrneig=0
	meascnt=0
	ineigh=0
	inegmeas=0
	ibadped=0
	ibaddark=0
	do j=1,64
	do i=1,318
	ikill(i,j)=0
	enddo
	enddo
c**********************************
	do j=1,64
	darkslope=2.*float(j-32)/32.
	do i=4,5
	if(iframe(i,j).ne.0)then
		iarg=int4(iframe(i,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)/0.75
		if(i.eq.4)then
		    if(abs(arg-pedlast).lt.2.0.or.k.eq.1)then
			pedsum=pedsum+arg
			ipedcnt=ipedcnt+1
		    endif	
		endif
		if(i.eq.5)then
		    if(abs(arg-pedlast-darklast-darkslope).lt.13.88.or.k.eq.1)then
			darksum=darksum+arg
			idarkcnt=idarkcnt+1
		    endif
		endif
	endif
	enddo
	enddo
	do j=1,63
	darkslope=2.*float(j-32)/32.
	do i=317,318
	if(iframe(i,j).ne.0)then
		iarg=int4(iframe(i,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)/0.75
		if(i.eq.318)then
		    if(abs(arg-pedlast).lt.2.0.or.k.eq.1)then
			pedsum=pedsum+arg
			ipedcnt=ipedcnt+1
		    endif	
		endif
		if(i.eq.317)then
		    if(abs(arg-pedlast-darklast-darkslope).lt.13.88.or.k.eq.1)then
			darksum=darksum+arg
			idarkcnt=idarkcnt+1
		    endif
		endif
	endif
	enddo
	enddo
c
	if(ipedcnt.lt.80)ibadped=1
	if(idarkcnt.lt.20)ibaddark=1

	if(ipedcnt.gt.0) ped=pedsum/float(ipedcnt)
	if(idarkcnt.gt.0) dark=darksum/float(idarkcnt)-ped
	
	peddiff=ped-pedlast
c	print *,k,ped,dark

	if((ibadped.eq.0.and.ibaddark.eq.0).or.k.le.10)then
		pedlast=ped
		darklast=dark
		pedrecover=0
		darkrecover=0
	else	
c		print*,'Bad ped or dark ',k,ped,dark,ipedcnt,idarkcnt,infile
		trash=trash+1
		if(ibadped.eq.1)then
			print*,'BAD PEDESTAL',k,'ped= ',ped,',   ipedcnt= ',ipedcnt,infile
			pedrecover=pedrecover+1			
			if(pedrecover.gt.12)then		 !don't lose more than 12 frames consecutively
				pedlast=pedinit
				if(mod(k,2).eq.0)pedlast=pedinit2	!vary recovery pedlast value
				pedrecover=0
				print*, 'Attempting Pedestal Recovery '
			endif
		endif
		if(ibaddark.eq.1)then
			print*,'BAD DARK CURRENT',k,'dark= ',dark,'   idarkcnt= ',idarkcnt,infile
			darkrecover=darkrecover+1
			if(darkrecover.gt.12)then		 !don't lose more than 12 frames consecutively
				darklast=darkinit
				if(mod(k,2).eq.0)darklast=darkinit2	!vary recovery darklast value
				darkrecover=0
				print*, 'Attempting Dark Current Recovery '
			endif
		endif
		go to 98
	endif
c************************************************************
	if(k.gt.10)go to 96
		do j=1,64
		do i=5,317
		if(iframe(i,j).ne.0.and.(i.ne.317.or.j.ne.64))then
			iarg=int4(iframe(i,j))
			if(iarg.lt.0)iarg=iarg+65536
			arg=float(iarg)/0.75
			if(i.eq.6.or.i.eq.316)arg=arg*1.3333
			arg=arg-dark-ped
			if(k.eq.1)then
				pattern(i,j)=pattern(i,j)+arg
				ipcnt(i,j)=ipcnt(i,j)+1
			else
			    arg1=arg-pattern(i,j)/float(ipcnt(i,j))
				if(k.eq.2)then
					if(arg1.lt.-patcr)then
					icrcnt=icrcnt+1
c					print *,'1st one bad',k,i,j,pattern(i,j),arg
						pattern(i,j)=arg
						go to 50
					endif
					if(arg1.gt.+patcr)then
					icrcnt=icrcnt+1
c					print *,'2nd one bad',k,i,j,arg,pattern(i,j)
					go to 50
					endif
				else
				if(arg1.gt.+patcr)then
				icrcnt=icrcnt+1
c				print *,'>=3rd one bad',k,i,j,arg,pattern(i,j)/float(ipcnt(i,j))
				go to 50
				endif
				endif
c				
				pattern(i,j)=pattern(i,j)+arg
				ipcnt(i,j)=ipcnt(i,j)+1
50				continue			
			endif
		endif
		if(k.eq.10.and.ipcnt(i,j).gt.0)pattern(i,j)=pattern(i,j)/float(ipcnt(i,j))
		enddo
		enddo
		if(k.eq.10)print *,'Done forming initial pattern, ',icrcnt,' cosmic rays removed'
		go to 98
96	continue	
**************************************************************
c	do j=1,jd
c	do i=1,id
c		frame(i,j)=float(iframe(i,j))/0.75
c		if(frame(i,j).ne.0.and.frame(i,j).le.fleast(i,j))fleast(i,j)=frame(i,j)
c	enddo
c	enddo
c*******
	pixsum=0.
	ipixsum=0
	pixdif=0.
	ipixdif=0	
	do j=1,64
	do i=1,318
		iarg=int4(float(iframe(i,j)))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)/0.75
		if(i.eq.6.or.i.eq.316)arg=arg*1.3333
		arg=arg-ped
		if(j.gt.4.and.j.lt.318)arg=arg-dark-pattern(i,j)
		frame(i,j)=arg!	-fleast(i,j)	!remove bland background
		if(iarg.eq.0)frame(i,j)=0.
		framein(i,j)=frame(i,j)		!for unremoved measles/C.R. file
	enddo
	enddo
c*******	
	call INDEXX(n,frame,indx)
	call measles(n,frame,indx,meascnt,ineigh,ihist,ikill)
	do j=1,64
		do i=1,318
c			if(fleast(i,j).eq.100000)then
c				fleast(i,j)=0
c				frame(i,j)=-10
c			endif
			if(frame(i,j).ne.0.)then
c			frame(i,j)=frame(i,j)+fleast(i,j)	!replace bland bckgrnd 
			endif
			if(ikill(i,j).eq.1)frame(i,j)=0.	!KILL measles pixels
			frameout(i,j)=frame(i,j)
			if(frame(i,j).ne.0..and.frame(i,j).lt.5000..and.i.gt.1.and.j.gt.1)then
				ipixsum=ipixsum+1
				pixsum=pixsum+frame(i,j)
				if(frame(i-1,j).ne.0..and.frame(i-1,j).lt.5000..and.i.gt.6.and.i.lt.317)then
					ipixdif=ipixdif+1
					pixdif=pixdif+abs(frame(i-1,j)-frame(i,j))
				endif
				if(frame(i,j-1).ne.0..and.frame(i,j-1).lt.5000..and.i.ge.6.and.i.lt.317)then
					ipixdif=ipixdif+1
					pixdif=pixdif+abs(frame(i,j-1)-frame(i,j))
				endif
			endif
		enddo	
	enddo
c
	pixsum=pixsum/float(ipixsum)
	pixdif=pixdif/float(ipixdif)
	
	iarg=mod(nint((t-t0)/4.),1524)+1
	pixsum=pixsum-minresponse(iarg)					!for use when minresponse table read in
	pixdif=pixdif-mingrit(iarg)					!for use when mingrit table read in
c	if(minresponse(iarg).gt.pixsum) minresponse(iarg)=pixsum	!for use in forming minresponse table
c	if(mingrit(iarg).gt.pixdif) mingrit(iarg)=pixdif		!for use in forming mingrit table
	pixsum=pixsum-0.025*float(meascnt)				!for use when minresponse table read in
	cut=pixdif+float(meascnt)/20.				!cut at line between pixdif=10 and meascnt=200
	if(cut.gt.10.)then
c		print 6,infile,meascnt,ineigh,pixdif
6		format(' file = ',a46,': Bad measles + grit, ',2i4,f10.1)
		trash=trash+1
		go to 98
	endif
c	do j=1,64
c		do i=6,316
c			if(frame(i,j).gt.0.)frame(i,j)=frame(i,j)-0.025*float(meascnt)
c		enddo
c	enddo

	
	dddate=ddate(1)
	if(k.gt.18262)dddate=ddate(2)
	if(pixsum.gt.3.0)then
		trash=trash+1
		write(19,19)dddate,h,hh,m,mm,s,ss,meascnt,pixsum,pixdif,k
		print 19,dddate,h,hh,m,mm,s,ss,meascnt,pixsum,pixdif,k
19		format('"',a8,1x,2i1,':',2i1,':',2i1,'"',2x,i6,2f10.2,i6)
		go to 98
	endif		
	write(12,12)dddate,h,hh,m,mm,s,ss,ped-pedinit,dark,ipedcnt,idarkcnt,meascnt,ineigh,pixsum,pixdif
	if(mod(k,100).eq.0)then
	print 12,dddate,h,hh,m,mm,s,ss,ped-pedinit,dark,ipedcnt,idarkcnt,meascnt,ineigh,pixsum,pixdif,k
12	format('"',a8,1x,2i1,':',2i1,':',2i1,'"',2x,2f10.3,4i6,2f10.2,i6)
	endif
c	write(out1file(24:33),'(a10)')infile(24:33)
c	write(out2file(24:33),'(a10)')infile(24:33)
c	write(out3file(24:33),'(a10)')infile(24:33)
c	open (14,file=out1file)
c	open (15,file=out2file)
c	open (16,file=out3file,form='binary')
c	write(14,140)
c	write(15,140)
140	format('DSAA'/'318 64'/'0 318'/'0 64'/'-10 1000')
c	write(14,15)((framein(i,j),i=1,318),j=1,64)
c	write(15,15)((frameout(i,j),i=1,318),j=1,64)
14	format(20i6)
	head(3)=+4
c	write(16)head,((frame(i,j),i=1,318),j=1,64),trail1,ped,dark,pixdif,trail2,ipedcnt,idarkcnt,inegmeas,
c     &	meascnt,ineigh,icrcnt,icrneig,trail3
c     	close(14)
c	close(15)
c	close(16)	
c	print*,k,'ped = ',ped, 'dark = ',dark,'measles = ',meascnt	

98	continue
	enddo
c	open (13,file='pattern2.grd')
c	write(13,13)
13	format('DSAA'/'318 64'/'0 318'/'0 64'/'-10 100')	
c	open (13,file='least2.grd')
c	write(13,13)
c13	format('DSAA'/'318 64'/'0 318'/'0 64'/'0 70')
c	do j=1,64
c	write(13,15)(pattern(i,j),i=1,318)
15	format(20f9.1)
c	write(13,15)(fleast(k,j),k=1,318)
c	enddo
c	close(13)
c	
	close(10)
c	close(12)
c	print *,'Start writing output files'
c	open (17,file='min2_4.dat')
c	do i=1,1524
c	write(17,17)i,minresponse(i)
c17	format(i6,f10.2)
c	enddo
c	close(17)	
c	open (17,file='mingrit2_4.dat')
c	do i=1,1524
c	write(17,17)i,mingrit(i)
c	enddo
c	close(17)
	close(19)	
	
	print *,'reached end,',trash,' discarded out of ',kmax-kmin+1
	stop 
	end
c********************************************************************************************
	subroutine measles(n,frame,indx,meascnt,ineigh,ihist,ikill)
	implicit none
c
c	This subroutine is a version of a cosmic-ray spot remover which uses the ADUs departure of
c	a pixel from its neighbors-average z0 to trigger its replacement/elimination. The response must 
c	be above a fixed value "dcut" and be at least double the 8-pixel average response z0. Also, the
c	departure from the neighbors average is compared to the difference of neighboring pixels on opposing 
c	sides of the center pixel.
c 		
	implicit none
	integer*4 i,j,ii,jj,meascnt,neigh,ineigh,nfix(3,3),ihist(250),iarg,n,m
	integer*4 indx(20352),ikill(318,64)
	real*4 frame(318,64),z(3,3),z0,z0sv
	real*4 dcut,sumzsq,devcut
	
	data dcut /13.88/				!13.88 is a 1000-electron cut
	data devcut /400./				!20 standard deviation cut
	
	meascnt=0

	do m=1,n
		iarg=1+n-m
		j=(indx(iarg)-1)/318
		i=indx(iarg)-318*j
		j=j+1
		if(i.le.6.or.i.ge.316.or.j.lt.2.or.j.gt.63)go to 98	!bail if near edges
				
		if(frame(i,j).ne.0.and.frame(i,j).ne.-10..and.frame(i,j).lt.dcut)go to 99	!quit if below threshold


		do jj=1,3
		do ii=1,3
			if(frame(i-2+ii,j-2+jj).eq.0.)go to 98		!bail if 3x3 has any empty pixel
			z(ii,jj)=frame(i-2+ii,j-2+jj)
		enddo
		enddo
		call plane(z,z0)
		if(z(2,2).lt.dcut)go to 98		!bail if not above threshold
		do jj=1,3
		do ii=1,3
			if(z(ii,jj).gt.z(2,2))go to 98
		enddo
		enddo
		sumzsq=0.
		do jj=1,3
		do ii=1,3
			if(ii.ne.2.and.jj.ne.2)sumzsq=sumzsq+z(ii,jj)*z(ii,jj)
		enddo
		enddo
		sumzsq=sumzsq/5.	!note altho 8 contribute, 3 params are fit...
		if(z(2,2)*z(2,2)/sumzsq.lt.devcut)go to 98
		z0sv=z0
		meascnt=meascnt+1
c		print *,i,j,'measle',z(2,2),z(2,2)/dev,dev
		neigh=0
		do jj=1,3
		do ii=1,3
		    nfix(ii,jj)=0
		    if((ii.ne.2.or.jj.ne.2).and.z(ii,jj).gt.dcut)then
		        z0=z0-z(ii,jj)/8.	!take away the neighbor's contribution
		        nfix(ii,jj)=1 
		        neigh=neigh+1
		    endif
		enddo
		enddo
		nfix(2,2)=1
		if(neigh.le.3)then
		    do jj=1,3
		    do ii=1,3
		        if(nfix(ii,jj).eq.1)then
		        	frame(i-2+ii,j-2+jj)=z0 	!Put in neighbors average...
		        	ikill(i-2+ii,j-2+jj)=1
		        endif
		    enddo
		    enddo
		    if(neigh.gt.0)ineigh=ineigh+1
		else				!too many neighbors, just fix original pixel
		    frame(i,j)=z0sv
		    ikill(i,j)=1
		endif
98	continue		
	enddo
99	continue
c	print *,'measles m = ',m
	return
	end
c*****************************************************************************
	subroutine plane(z,z0)
	implicit none
c
c	This subroutine solves for a least-squares planar fit through the 8 z-values
c	surrounding [2,2], then reduces all 9 z's relative to the plane,
c	and returns with the reduced values.
c
	real*4 z(3,3),sumz,sumx,sumy,a,b,z0
	integer*4 i,j
c
	sumz=z(1,2)+z(3,2)
	sumx=0.
	sumy=0.
	do i=1,3
		sumz=sumz+z(i,1)+z(i,3)
		sumx=sumx+z(1,i)-z(3,i)
		sumy=sumy+z(i,1)-z(i,3)
	enddo
	a=sumx/6.
	b=sumy/6.
	z0=sumz/8.
	do i=1,3
	    do j=1,3
		z(i,j)=z(i,j)-z0+a*float(i-2)+b*float(j-2)
	    enddo
	enddo
	return
	end
c*********************************************************************************
      SUBROUTINE INDEXX(N,ARRIN,INDX)	!from Numerical Receipes...
      implicit none
      integer*4 INDX(20352),N,I,J,L,IR,INDXT	!Absoft Compiler dislikes INDX(N), etc.
      real*4 ARRIN(20352),Q
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END