c
c name:		CRX_ped_dark_cam2.for
c purpose:	determines 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 
c    "		Note here that the pedestal, dark, and CR information incorporated into the "trailer"
c    "		The output file xxx.dat presents information for graphs of the various quantities
c calls:	subroutines  measles (+ plane), and INDEXX 
c mods:		November 2003, Andrew Buffington (UCSD; (858)-534-6630; abuffington@ucsd.edu)
c-
	program CRX_ped_dark_cam2
	implicit none
	integer*1 trail1(75),trail2(275),tdum1(12),trail3(21),tdum2(28),trail4(69)
	integer*2 iframe(318,64),head(3)
	integer*4 h,m,s,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,icrx,ibin,igraph,iffix
	integer*4 indx(20352),k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,hh,mm,ss
	integer*4 mode,id,jd,offset,ii,jj,icorrcnt,jarg,iii,jjj,iqbad,iqcnt,iorbit,iarglast
	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),arg,arg1,patcr
	real*4 minresponse(1524),mingrit(1524),t,t0,pixsum,pixdif,acorr(1280,301),acorrsum,flatcorr(318,64)
	real*4 torbit,tbegin,targ,pednaive
c	real*4 cut
	real*8 q(4),qq(4)
c	character out3file*48,out2file*48,out1file*48
	character infile*64
	character doy(15)*3,ddoy*3
	character nname*15,outname*25
	character hhead(5)*10
	character ddate(15)*8,dddate*8,pauldate*10
	data infile   /'../../../../../mnt/big/Orbits/cam2_DOY/c2frm_2003_DOY_hhmmss.nic'/
	data outname  /'c2frm_2003_DOY_hhmmss.nic'/
c	data out1file /'G:/SEAN/buffer2/Crx/c2frm_2003_DOY_hhmmsscr1.grd'/	!cosmic-ray not removed file...
c	data out2file /'G:/SEAN/buffer2/Crx/c2frm_2003_DOY_hhmmsscrx.grd'/	!cosmic-ray removed file...
c	data out3file /'G:/SEAN/buffer2/Bin/c2frm_2003_DOY_hhmmsscrx.bin'/	!cosmic-ray removed binary file...
	
	data ihist /250*0/
	data i1,i2 /228,301/
	data ddate /'05/24/03','05/25/03','05/26/03','05/27/03','05/28/03','05/29/03','05/30/03','05/31/03','06/01/03',
     &	'06/02/03','06/03/03','06/04/03','06/05/03','06/06/03','06/07/03'/
	data doy /'144','145','146','147','148','149','150','151','152','153','154','155','156','157','158'/
	data badsum /0./
	data pattern /20352*0./
	data ipcnt /20352*0/
	data n /20352/
	data trash,iffix,iqcnt /3*0/
	data patcr /6.94/	!6.94 adu's = 500 electrons
	data minresponse /1524*100000./
	data mingrit /1524*100000./
	data torbit /6096.297/	!value determined from quaternion analysis (see orbitcompare.for & badquarts.for)
c	data tbegin /3900./	!starts orbit #1 while cam 2 looks at -90 deg dec. (see orbitcompare.for output)
	data tbegin /0/		!starts orbit #1 at DOY144 0h 24m 51s, presumably the same as above
	data iarglast /1524/
c
	icrx=0			!icrx=1 means make before and after files
	ibin=0			!ibin=1 means make binary CRX files
	igraph=1		!igraph=1 means make GRAPHER output file
	mode = 4	
	id = 1272
	jd = 256
	id=id/mode
	jd=jd/mode
c
c	Bring in Large Scale Flatfield Table, fill empty entries with neighbors average
c	Offsets for cameras 1, 2, and 3 are respectively 16, 10, and 12 (see Table I, memo dated 20 April 2001)
c	Here we add one to these each because "offset" here is the first column
c
	offset=11	!Camera 2
	print *, 'reading flatfield table...'
	open (10,file='DATA_SMEI/newfltgrd2.GRD',readonly)
	read (10,1)hhead
1	format(a10)
	read (10,3)((acorr(i,j),i=1,1280),j=1,301)
3	format(20f10.5)
	do i=5,id-5
	do j=5,jd-5
		iarg=2+mode*(i-1)
		jarg=offset+mode*(j-1)
		if(acorr(iarg,jarg).eq.0.)then
			acorrsum=0.
			icorrcnt=0
			do ii=-1,1
			do jj=-1,1
				if(ii.ne.0.or.jj.ne.0)then
					if(acorr(iarg+mode*ii,jarg+mode*jj).ne.0.)then
						acorrsum=acorrsum+acorr(iarg+mode*ii,jarg+mode*jj)
						icorrcnt=icorrcnt+1
					endif
				endif
			enddo
			enddo
			if(icorrcnt.ge.6)then
				acorrsum=acorrsum/float(icorrcnt)
				iffix=iffix+1
				do iii=0,mode-1
				do jjj=0,mode-1
				acorr(iarg+iii,jarg+jjj)=acorrsum
				enddo
				enddo
c				print *,iarg,jarg+49,', replaced with ',acorrsum
			endif
		endif
	enddo
	enddo
	close(10)
	print *,iffix,' large scale flatfield entries filled in with neighbors average'	
	do i=1,id
	do j=1,jd
		iarg=2+mode*(i-1)			!Want every 4th column starting with (2)
		jarg=offset+mode*(j-1)			!Want every 4th row starting with (11)
		flatcorr(i,j)=acorr(iarg,jarg)
	enddo
	enddo		
c
c	Bring in minimum response table
c
	print *, 'reading minimum response table in...'
	open (17,file='DATA_SMEI/min2_4.dat',readonly)	
	open (18,file='DATA_SMEI/mingrit2_4.dat',readonly)
c		
	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)	
c	
	open (10,file='DATA_SMEI/namelist2.txt',readonly)
	open (30,file='DATA_SMEI/mayevent_cam2_fixed.dat',readonly)
	if(igraph.eq.1)open (12,file='cam2_144to157_CRXd.dat')
	if(igraph.eq.1)open (21,file='cam2_144to157_CRXd.txt')
c
	pedrecover=0
	pedinit=1016.
	pedinit2=1020.
	darkinit=4.5
	darkinit2=8.5
	darkrecover=0
	icrcnt=0
c********************************************************************
	kmax =  271323	!73921
	k1 = 	10097	!start day 145
	k2 = 	31052	!146
	k3 = 	52022	!147
	k4 = 	72972	!148
	k5 = 	93912	!149
	k6 = 	114862	!150
	k7 = 	135746	!151
	k8 = 	156671	!152
	k9 = 	177600	!153
	k10 = 	198524	!154
	k11 = 	219447	!155
	k12 = 	240374	!156
	k13 = 	261294	!157
	k14 = 	271324	!17:29:47 of day 157 (end of data frames)
c
c	Here make the above list compliant with the first 10 frames being pattern frames
c
c	kmax=kmax+10
c	k1=k1+10
c	k2=k2+10
c	k3=k3+10
c	k4=k4+10
c	
	kmin=1	
	do k=1,kmax
	if(k.eq.11)then
		pedlast=pedinit
		darklast=darkinit
	endif
	if(k.lt.k4)read (10,10)nname,h,m,s	!dir on Andy's machine makes a different list than ls on Aaron's machine
	if(k.ge.k4)read (10,100)nname,h,m,s
10	format(39x,a15,3i2.2)
100	format(a15,3I2.2)
c***
c	if(k.eq.10)then
c	nname='c2frm_2003_158_'
c	h=2
c	m=17
c	s=7
c	endif
c***
c32	if(k.gt.10)then
c		read (30,30)pauldate,hh,mm,ss,qq
c30		format(a10,i2.2,1x,i2.2,1x,i2.2,1x,44x,4f13.9)
c		if(hh.ne.h.or.mm.ne.m.or.ss.ne.s)then
c			print 31,pauldate,hh,mm,ss,nname,h,m,s
c31			format(a10,3i2,', looking for ',a15,3i2)
c			go to 32
c		endif
c	endif			
	t=3600.*float(h)+60.*float(m)+float(s)
	if(k.ge.k1)t=t+86400.
	if(k.ge.k2)t=t+86400
	if(k.ge.k3)t=t+86400.
	if(k.ge.k4)t=t+86400
	if(k.ge.k5)t=t+86400.
	if(k.ge.k6)t=t+86400
	if(k.ge.k7)t=t+86400.
	if(k.ge.k8)t=t+86400
	if(k.ge.k9)t=t+86400.
	if(k.ge.k10)t=t+86400
	if(k.ge.k11)t=t+86400.
	if(k.ge.k12)t=t+86400
	if(k.ge.k13)t=t+86400
c	if(k.ge.k14)t=t+86400
	if(k.eq.11)t0=t
	if(k.gt.11.and.k.lt.kmin)go to 98
	if(iabs(k-8).le.1)go to 98		!three open-shutter frames, kill...
c***	
	dddate=ddate(1)
	ddoy=doy(1)
c	if(k.eq.10)ddoy='158'
c	if(k.ge.11)ddoy=doy(2)
	if(k.ge.k1)dddate=ddate(2)
	if(k.ge.k1)ddoy=doy(2)
	if(k.ge.k2)dddate=ddate(3)
	if(k.ge.k2)ddoy=doy(3)
	if(k.ge.k3)dddate=ddate(4)
	if(k.ge.k3)ddoy=doy(4)
	if(k.ge.k4)dddate=ddate(5)
	if(k.ge.k4)ddoy=doy(5)
	if(k.ge.k5)dddate=ddate(6)
	if(k.ge.k5)ddoy=doy(6)
	if(k.ge.k6)dddate=ddate(7)
	if(k.ge.k6)ddoy=doy(7)
	if(k.ge.k7)dddate=ddate(8)
	if(k.ge.k7)ddoy=doy(8)
	if(k.ge.k8)dddate=ddate(9)
	if(k.ge.k8)ddoy=doy(9)
	if(k.ge.k9)dddate=ddate(10)
	if(k.ge.k9)ddoy=doy(10)
	if(k.ge.k10)dddate=ddate(11)
	if(k.ge.k10)ddoy=doy(11)
	if(k.ge.k11)dddate=ddate(12)
	if(k.ge.k11)ddoy=doy(12)
	if(k.ge.k12)dddate=ddate(13)
	if(k.ge.k12)ddoy=doy(13)
	if(k.ge.k13)dddate=ddate(14)
	if(k.ge.k13)ddoy=doy(14)
c	if(k.ge.k14)dddate=ddate(15)
c	if(k.ge.k14)ddoy=doy(15)
c
	write(infile(36:38),'(a3)')ddoy
	write(infile(40:54),'(a15)')nname
	write(infile(55:60),'(3i2.2)')h,m,s
	write(outname(1:15),'(a15)')nname
	write(outname(16:21),'(3i2.2)')h,m,s
c	if(k.eq.k1)print *,infile
c	if(k.eq.k2)print *,infile
c	if(k.eq.k3)print *,infile
c	if(k.eq.k4)print *,infile
c	if(k.eq.k5)print *,infile
c	if(k.eq.k6)print *,infile
c	if(k.eq.k7)print *,infile
c	if(k.eq.k8)print *,infile
c	if(k.eq.k9)print *,infile
c	if(k.eq.k10)print *,infile
c	if(k.eq.k11)print *,infile
c	if(k.eq.k12)print *,infile
c	if(k.eq.k13)print *,infile
c	if(k.eq.k14)print *,infile
	open(11,file=infile,form='binary',readonly)
c	if(k.ge.11.and.k.lt.k4-4000)then
c	  close(11)
c	  go to 98
c	endif
	read(11)head,((iframe(i,j),i=1,318),j=1,64),trail1,q,trail2,tdum1,trail3,tdum2,trail4
c	if(k.eq.k1-1)print 20,k,dddate,h,m,s
c	if(k.eq.k2-1)print 20,k,dddate,h,m,s
c	if(k.eq.k3-1)print 20,k,dddate,h,m,s
c	if(k.eq.k4-1)print 20,k,dddate,h,m,s
c20	format('k = ',i6,'  "',a8,1x,i2.2,':',i2.2,':',i2.2,'"')
	close(11)
c	if(mod(k,1000).eq.0)print *,infile
c	if(k.ge.11)go to 98		!useful when making a quick run to check presence of frames
c	print*, infile
c	iqbad=0
c	do i=1,4
c		if(k.ge.11.and.dabs(q(i)-qq(i)).gt.1.0D-06)then
c			iqbad=iqbad+1
c		endif
c	enddo
c	if(iqbad.gt.0)then
c		do i=1,4
c			q(i)=qq(i)
c		enddo
c		print *,k,': Replaced bad quaternions on ',nname,h,m,s
c		iqcnt=iqcnt+1
c	endif		
	ped=0.
	dark=0.
	idarkcnt=0
	darksum=0.
	ipedcnt=0
	pedsum=0.
	pednaive=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
		    pednaive=pednaive+arg
		    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
		    pednaive=pednaive+arg
		    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
	pednaive=pednaive/127.
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,pednaive,', 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
			ped=-9999.
			dark=-9999.
		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
			dark=-9999.
		endif
		go to 98
	endif
c************************************************************
	if(k.gt.10)go to 96
		if(k.eq.1)print *,' forming pattern from first 10 frames'
		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
c***			else
			elseif(k.le.9)then
			    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))
c***		
		if(k.eq.10)then
		if(arg.lt.pattern(i,j)-10.)then
		print *, i,j,arg,pattern(i,j)
		pattern(i,j)=arg
		endif
		endif
c***		
		enddo
		enddo
		if(k.eq.10)print *,'Done forming initial pattern, ',icrcnt,' cosmic rays removed'
		go to 98
96	continue	
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
		if(iarg.eq.0)frame(i,j)=0.
		if(icrx.eq.1)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
			if(ikill(i,j).eq.1)frame(i,j)=0.	!KILL measles pixels
			if(icrx.eq.1)frameout(i,j)=frame(i,j)	!for removed measles/C.R. file
			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.ge.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)
c	
	targ=(t-t0-tbegin)/torbit
	iorbit=int(targ)+1
	if(targ.lt.0.)then
		targ=targ+1.
		iorbit=iorbit-1
	endif
	targ=amod(targ,1.0)
	iarg=1+int(1524.*targ)
	if(iarg.gt.1524)iarg=1524
	if(iarg.lt.iarglast)print 111,iorbit,k,infile
111	format('Beginning orbit #',I3.3,' at k = ',I6,', ',A64)
	iarglast=iarg
	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
c
c	cut=pixdif+float(meascnt)/20.					!cut at line between pixdif=10 and meascnt=200
c	if(cut.gt.10.)then
c		print 6,infile,meascnt,ineigh,pixdif
6		format(' file = ',a46,': Bad measles + grit, ',2i4,f10.1)
c		trash=trash+1
c		go to 98
c	endif
	do j=1,64
		do i=6,316
c			if(frame(i,j).gt.0.)frame(i,j)=frame(i,j)-0.025*float(meascnt)
			if(flatcorr(i,j).ne.0.)frame(i,j)=frame(i,j)/flatcorr(i,j)  !Flat Field Correction
		enddo
	enddo
c***		Old aurora condition follows	
c		if(pixsum.gt.3.0)then
c		write(19,19)dddate,h,hh,m,mm,s,ss,meascnt,pixsum,pixdif
c		print 19,dddate,h,hh,m,mm,s,ss,meascnt,pixsum,pixdif,k
c19		format('"',a8,1x,2i1,':',2i1,':',2i1,'"',2x,i6,2f10.2,i6)
c	endif		
	if(igraph.eq.1)write(12,12)dddate,h,m,s,ped-pedinit,dark,meascnt,pixsum,pixdif
	if(igraph.eq.1)write(21,33)outname,ddoy,ped,dark,q,meascnt,pixsum,pixdif
	if(mod(k,1000).eq.0)then
	print 12,dddate,h,m,s,ped-pedinit,dark,meascnt,pixsum,pixdif,k
12	format('"',a8,1x,i2.2,':',i2.2,':',i2.2,'"',2x,2f10.3,i6,2f10.2,i6)
33	format(a25,a3,2f10.3,4f13.9,i6,2f10.2)
	endif
c	write(out1file(21:41),'(a21)')infile(24:44)
c	write(out2file(21:41),'(a21)')infile(24:44)
c	write(out3file(21:41),'(a21)')infile(24:44)
c	if(icrx.eq.1)open (14,file=out1file)
c	if(icrx.eq.1)open (15,file=out2file)
c	if(ibin.eq.1)open (16,file=out3file,form='binary')
	if(icrx.eq.1)write(14,140)
	if(icrx.eq.1)write(15,140)
140	format('DSAA'/'318 64'/'0 318'/'0 64'/'-10 1000')
	if(icrx.eq.1)write(14,15)((framein(i,j),i=1,318),j=1,64)
	if(icrx.eq.1)write(15,15)((frameout(i,j),i=1,318),j=1,64)
	head(3)=+4
c								trail1,q,trail2,tdum1,trail3,tdum2,trail4
     	if(ibin.eq.1)write(16)head,((frame(i,j),i=1,318),j=1,64),trail1,q,trail2,ped,dark,pixdif,trail3,ipedcnt,idarkcnt,inegmeas,
     &	meascnt,ineigh,icrcnt,icrneig,trail4
     	if(icrx.eq.1)close(14)
	if(icrx.eq.1)close(15)
	if(ibin.eq.1)close(16)	
c	print*,k,'ped = ',ped, 'dark = ',dark,'measles = ',meascnt

98	continue
	enddo
	print*,'end'
c	go to 99
	close(10)
	close(30)
	if(igraph.eq.1)close(12)
	if(igraph.eq.1)close(21)
c	print *,'Start writing output files'
c	open (13,file='pattern2.grd')
c	print *,' writing out pattern2.grd'
c	write(13,13)
13	format('DSAA'/'318 64'/'0 318'/'0 64'/'-10 100')	
c	do j=1,64
c	write(13,15)(pattern(i,j),i=1,318)
15	format(20f9.3)
c	enddo
c	close(13)
c	
c	open (17,file='min2_4.dat')
c	do i=1,1524
c	write(17,17)i,minresponse(i)
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)
c
	close(20)	
	print *,'reached end,',trash,' discarded out of ',kmax-kmin+1, ', quaternions replaced on ',iqcnt
	stop 
99	continue
	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 		
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