c
c name:		ENG_pattern_cam3.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:		May 2004, Andrew Buffington (UCSD; (858)-534-6630; abuffington@ucsd.edu)
c		Please Note that "headroom" and small-scale FF are turned OFF when camera 3 in 1x1 mode as here...
c		Another consequence of this is that columns 21 and 1264 are NOT stamped out as would be
c		the case if the SSFF were imposed. Thus these have a full amount of pedestal and dark current,
c		but only an unknown and variable-with-row fraction of sky between 0.5 and unity. Yuck!
c-
	program ENG_pattern_cam3
	implicit none
	integer*1 trail1(75),trail2a(94),trail2b(177),tdum1(12),trail3(21),tdum2(28),trail4(69)
	integer*2 iframe(1272,256),head(3)
	integer*4 h,m,s,i,j,k,kmax,kmin,iarg,ineigh,ptrash,dtrash,ipcnt(1272,256)
	integer*4 idarkcnt,ipedcnt,ibadped,ibadark,ihist(256),icrcnt,meascnt,ikill(1272,256)
	integer*4 icrneig,n,inegmeas,ipixsum,ipixdif,igraph,iffix,isunmoon,imed
	integer*4 indx(325632),indxx(2044),k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14
	integer*4 mode,id,jd,jarg,med1,med2,isat,ipatt
	real*4 dark,darklast,darksum,ped,ped1,ped2,pediff,pedlast,pedinit,pedright,pedsum,ddark
	real*4 frame(1272,256),badsum,pattern(1272,256),arg,arg1,patcr
c	real*4 pednaive,patdark,patcorr,patternbin(636,128),darkdist(2044,13030),darkmed(2044)
	real*4 pednaive,patdark,patcorr,patternbin(636,128),darkdist(2044,4000),darkmed(2044)
c	real*4 pednaive,patdark,patcorr,patternbin(636,128),darkdist(2044,8000),darkmed(2044)
	real*4 pedsv(2044),darksv(2044),median,dmed,pedtest,pixsum,pixdif
	real*4 camtemp
	real*8 q(4)
	character infile*43,infilesv(10)*43
	character doy(15)*3,ddoy*3 
	character nname*15,outname*25
	character ddate(15)*8,dddate*8
	data infile   /'D:/Buffer3/NicDOY/c3frm_2003_DOY_hhmmss.nic'/
	data outname  /'c3frm_2003_DOY_hhmmss.nic'/		
	
c	data ihist /256*0/
c	data ddate /'02/04/03','02/05/03','03/16/03','03/20/03','03/25/03','03/31/03','04/01/03','04/03/03','04/05/03',
c     &	'04/16/03','04/24/03','05/12/03','05/22/03','05/23/03','05/24/03'/
c	data doy /'035','036','075','079','084','090','091','093','095','106','114','132','142','143','144'/
c	data ddate /'06/07/03','06/11/03','06/12/03','06/14/03','06/21/03','07/05/03','07/12/03','07/24/03','08/02/03',
c     &	'08/08/03','08/15/03','08/18/03','08/20/03','08/22/03','09/02/03'/
c	data doy /'158','162','163','165','172','186','193','205','214','220','227','230','232','234','245'/
	data ddate /'09/06/03','09/13/03','09/20/03','10/06/03','10/07/03','10/08/03','10/23/03','11/06/03',
     &	'11/07/03','11/15/03','11/22/03','11/29/03','12/10/03','12/20/03','12/26/03'/
	data doy /'249','256','263','279','280','281','296','310','311','319','326','333','344','354','360'/
c	data ddate /'12/29/03','12/31/03','01/08/04','01/15/04','01/16/04','01/17/04','01/24/04','02/05/04',
c     &	'02/14/04','02/22/04','02/24/04','03/07/04','03/17/04','03/24/04','04/01/04'/
c	data doy /'363','365','008','015','016','017','024','036','045','053','055','067','077','084','092'/
c	data doy /'363','365','373','380','381','382','389','401','410','418','420','432','442','449','457'/
c	data ddate /'04/08/04','04/24/04','05/01/04','01/15/04','01/16/04','01/17/04','01/24/04','02/05/04',
c     &	'02/14/04','02/22/04','02/24/04','03/07/04','03/17/04','03/24/04','04/01/04'/
c	data doy /'099','115','122','015','016','017','024','036','045','053','055','067','077','084','092'/
c	data doy /'464','480','487','380','381','382','389','401','410','418','420','432','442','449','457'/
	data badsum /0./
	data pattern /325632*0./
	data ipcnt /325632*0/
	data n /325632/
	data ptrash,dtrash,iffix,isunmoon /4*0/
	data patcr /44.44/	!222.22 adu's = 1000 electrons for 1 x 1 binning
	data patdark /0./
	data patternbin /81408*0./
c	data darkdist /26633320*0./
	data darkdist /8176000*0./
c	data darkdist /16352000*0./
	data darksv /2044*0./
c
	ipatt=0			!ipatt=1 means output the pattern *.grd files
	igraph=0		!igraph=1 means make GRAPHER output file and *.txt file for indexing...
	mode = 1	
	id = 1272
	jd = 256
	id=id/mode
	jd=jd/mode
	open (10,file='namelist3engc.txt',readonly)
c	if(igraph.eq.1)open (12,file='cam3_363to457_ped.dat')	!for pedestal run
	if(igraph.eq.1)open (12,file='cam3_249to310_CRX.dat')	!for CRX run
c
	pedinit=1046.
	icrcnt=0
	pedlast=pedinit		!these presuppose that the chosen pattern frames have these values
	pedright=0.0		!pedright is the amount by which the right-hand-side pedestal exceeds the left-hand one
********************************************************************
c			5	4	3	2	1	Chronological order of sections
	kmax=	1919	!516	!3096	!3479	!3177	!7889	!10 extra added in here for the initial pattern
	k1=	169	!172	!306	!169	!281	!337
	k2=	342	!345	!467	!342	!500	!846
	k3=	521	!588	!588	!521	!731	!1152
	k4=	722	!742	!742	!722	!936	!1394
	k5=	894	!912	!912	!894	!1214	!1625
	k6=	1597	!1082	!1082	!1597	!1477	!2746
	k7=	1770	!1252	!1252	!1770	!1721	!3515
	k8=	1920	!1938	!1938	!1920	!1911	!3980
	k9=	1920	!2102	!2102	!2084	!2069	!4440
	k10=	1920	!2263	!2263	!2261	!2373	!4644
	k11=	1920	!2431	!2431	!2433	!2530	!5830
	k12=	1920	!2601	!2601	!2603	!2675	!6404
	k13=	1920	!2766	!2766	!2771	!2842	!6577
	k14=	1920	!2933	!2933	!2932	!3008	!6931
c	kmax=10
c
	kmin=1597
	do k=1,kmax
	if(k.eq.11)then			!this presupposes that the data begin at these 2nd values
		pedlast=pedinit
		darklast=patdark
	endif
c
	read (10,10)nname,h,m,s
10	format(39x,a15,3i2.2)
c
	if(k.ge.11.and.k.lt.kmin)go to 98
	dddate=ddate(1)
	ddoy=doy(1)
	if(k.le.10)dddate=ddate(7)	!these put in appropriate DOY when pattern is from a later section
	if(k.le.10)ddoy=doy(7)		!this appropriate for first section of data
	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)
	if(k.ge.k14)dddate=ddate(15)
	if(k.ge.k14)ddoy=doy(15)
	write(infile(15:17),'(a3)')ddoy
	write(infile(19:33),'(a15)')nname
	write(infile(34:39),'(3i2.2)')h,m,s
	if(k.le.10)infilesv(k)=infile
	write(outname(1:15),'(a15)')nname
	write(outname(16:21),'(3i2.2)')h,m,s
c	print *,infile
	open(11,file=infile,form='binary',readonly)
	read(11)head,((iframe(i,j),i=1,1272),j=1,256),trail1,q,trail2a,camtemp,trail2b,tdum1,trail3,tdum2,trail4
	close(11)
	ped=0.
	ped1=0.
	ped2=0.
	dark=0.
	idarkcnt=0
	darksum=0.
	ipedcnt=0
	ibadped=0
	pedsum=0.
	pednaive=0.
	ibadark=0
	if(k.gt.10)icrcnt=0
	icrneig=0
	meascnt=0
	ineigh=0
	inegmeas=0
	isat=0
	do j=1,256
	do i=1,1272
	ikill(i,j)=0
	enddo
	enddo
c**********************************	Here begin calculation of electronic pedestal offset
	do j=1,256
	do i=13,16
	    if(iframe(i,j).ne.0)then
		iarg=int4(iframe(i,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)		!/0.75	!"headroom" restored: this occurs 3 more times below
		if(k.gt.10)arg=arg-pattern(i,j)	!UGH: a pattern even in the pedestal! Do not scale with dark current
		pednaive=pednaive+arg
		ped1=ped1+arg
		ipedcnt=ipedcnt+1
		pedsv(ipedcnt)=arg
	    endif
	enddo
	enddo
	do j=1,255
	do i=1269,1272
	    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(k.gt.10)arg=arg-pattern(i,j)
		if(k.gt.10)arg=arg-pedright
		pednaive=pednaive+arg
		ped2=ped2+arg
		ipedcnt=ipedcnt+1
		pedsv(ipedcnt)=arg
	    endif
	enddo
	enddo
	pednaive=pednaive/2044.
	if(ipedcnt.ne.2044.or.pednaive-pedinit.gt.50.)then
	if(pednaive-pedinit.le.50.)print *,infile,': warning, incomplete pedestal roster, icount = ',ipedcnt
c	print *,infile,' crazy raw pedestal = ',pednaive-pedinit,ipedcnt,k-10
		isunmoon=isunmoon+1
		go to 98
	endif
	ped1=ped1/1024.
	ped2=ped2/1020.
	pediff=ped2-ped1
c	print *,k-10,pednaive,ped1,ped2
c
c	Here update the differential correction to the right-hand-side pedestal
c	Use an exponential running average over 10 frames 
c	Exclude changes abs(pediff) > 4 ADUs keeps crazies out, maybe
c
	if(k.gt.10.and.abs(pediff).lt.4)pedright=pedright+0.1*pediff
	ipedcnt=0
	call indexxx(2044,pedsv,indxx)
	pedtest=pedsv(indxx(1022))		!here compare with the naive pedestal median
	do j=1,2044
	    if(pedsv(indxx(j))-pedtest.gt.4)go to 20		!bail at > 18 electrons...same for all cameras
	    ipedcnt=ipedcnt+1
	    pedsum=pedsum+pedsv(indxx(j))
	enddo
20	if(ipedcnt.lt.1800)ibadped=1
	ped=pedlast
	if(ipedcnt.gt.0)ped=pedsum/float(ipedcnt)
	print *,ipedcnt
c	do i=1,2044		!*********************************************************ped dist
c		darkdist(i,k)=pedsv(indxx(i))-pedtest
c	enddo
c	if(k.ge.11)go to 97	!*********************************************************ped dist
c*********************************************	Begin dark current calculation	
	do j=1,256
	do i=17,20
	if(iframe(i,j).ne.0)then
		iarg=int4(iframe(i,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)		!/0.75	!"headroom" restored: this occurs 3 more times below
   		jarg=i-16+8*(j-1)
		darkmed(jarg)=arg
		darksum=darksum+arg
		if(k.ge.11)then		!note here the CR removal differs for first 10; OK here, no CRs found!
			ddark=arg-pattern(i,j)*darklast/patdark-ped
			if(ddark.gt.2.*patcr)then		!CR cut for the pattern, = 400 e-
				idarkcnt=idarkcnt+1
				darkmed(jarg)=1000.
				darksum=darksum-arg
			endif
		endif
	endif
	enddo
	enddo
	do j=1,255
	do i=1265,1268
	if(iframe(i,j).ne.0)then
		iarg=int4(iframe(i,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)		!/0.75
		jarg=i-1260+8*(j-1)
		darkmed(jarg)=arg
		darksum=darksum+arg
		if(k.ge.11)then
			ddark=arg-pattern(i,j)*darklast/patdark-ped
			if(ddark.gt.2.*patcr)then		!CR cut for the pattern, = 400 e-
				idarkcnt=idarkcnt+1
				darkmed(jarg)=1000.
				darksum=darksum-arg
			endif
		endif
	endif
	enddo
	enddo
	if(idarkcnt.ge.684)ibadark=1	!we find imed < 680 sometimes fails to have a different value < the reduced median
	do j=1,2044
		if(k.le.10)darksv(j)=darksv(j)+(darkmed(j)-ped)/10.
		if(darkmed(j).ne.1000.)darkdist(j,k)=darkmed(j)-ped
	enddo
c
c	Get median:"interpolate within the histogram..."
c
	call indexxx(2044,darkmed,indxx)
c	do j=1,2044
c		darkdist(j,k)=darkmed(indxx(j))-ped	!now is sorted
c	enddo
	med1=0
	med2=0
	imed=ifix((2044.-float(idarkcnt))/2.0)
	dmed=darkmed(indxx(imed))
	do j=1,imed-1
	if(darkmed(indxx(imed-j)).ne.dmed)then
		med1=imed-j+1
		go to 51
	endif
	enddo
51	continue	
	do j=1,2044-imed
	if(darkmed(indxx(imed+j)).ne.dmed)then
		med2=imed+j
		go to 52
	endif
	enddo
52	continue
	if(ibadped.eq.0.and.ibadark.eq.0.and.(med1.eq.0.or.med2.eq.0))then
		print *,'Bad median determination ',ped-pedinit,med1,imed,med2,infile,k-10
		ibadark=1
	endif
	if(ibadped.eq.0.and.ibadark.eq.1)then
		dtrash=dtrash+1
		do j=1,2044
			darkdist(j,k)=0.
		enddo
		go to 98
	endif	
	median=darkmed(indxx(imed))-0.5+float(imed-med1)/float(med2-med1)-ped
	dark=darksum/float(2044-idarkcnt)-ped
	dark=median
	do j=1,2044
		if(k.ge.11.and.darkdist(j,k).eq.1000.)darkdist(j,k)=0.
		if(k.ge.11.and.darkdist(j,k).ne.0.)darkdist(j,k)=darkdist(j,k)-darksv(j)*(dark/patdark)
c		if(k.ge.11)darkdist(j,k)=darkdist(j,k)-darksv(j)
	enddo
c		
	if(k.lt.11)print *,k,infile,pednaive,ped,pediff,pedright,ipedcnt,dark,idarkcnt
c
	if(ibadped.eq.0.or.k.le.10)then
		pedlast=ped
		darklast=dark
	else	
c		print*,'Bad ped or dark ',ped-pedinit,dark,ipedcnt,2044-idarkcnt,infile,k-10
		ptrash=ptrash+1
c		print*,'BAD PEDESTAL',k,'ped= ',ped,pednaive,', ipedcnt= ',ipedcnt,infile
		ped=-9999.
		dark=-9999.
		do j=1,2044
			darkdist(j,k)=0.
		enddo
		go to 98
	endif
c************************************************************
	if(k.ge.11)go to 96
		if(k.eq.1)print *,' forming pattern from the chosen ten frames'
		infilesv(k)=infile
c		print *, infile, ped, dark, icrcnt
		patdark=patdark+dark
		do j=1,256
		do i=5,1272
		if(iframe(i,j).ne.0.and.(i.le.1267.or.j.ne.256))then
			iarg=int4(iframe(i,j))
			if(iarg.lt.0)iarg=iarg+65536
			arg=float(iarg)-ped		!/0.75
			if(k.eq.1)then
				pattern(i,j)=arg
				ipcnt(i,j)=1
			else
			    arg1=arg-pattern(i,j)/float(ipcnt(i,j))
				if(k.eq.2)then
					if(arg1.lt.-0.5*patcr)then	!2nd frame rules if less by 100 e- or more
						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.+0.5*patcr)then	!1st frame rules if 2nd more by 100e- or more
						icrcnt=icrcnt+1
c						print *,'2nd one bad',k,i,j,arg,pattern(i,j)
						go to 50
					endif
				else
					if(arg1.gt.+0.5*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)then
			patdark=patdark/10.
			print *,'Done forming initial pattern, ',icrcnt,' total cosmic rays removed, patdark = ',patdark
			if(ipatt.eq.1)then
			open (13,file='c3pat_2003_296_eng.grd')
			print *,' writing out c3pat_2003_296_eng.grd'
			write(13,13)
13			format('DSAA'/'1272 256'/'0 1272'/'0 256'/'-10 100')	
			do j=1,256
				write(13,15)(pattern(i,j),i=1,1272)
15				format(20f9.3)
			enddo
			write(13,16)patdark,(infilesv(i),i=1,10)
			close(13)
16			format(f10.4,/,10(a43,/))			
			do i=1,1272
				iarg=1+(i-1)/2
				do j=1,256
					jarg=1+(j-1)/2
					patternbin(iarg,jarg)=patternbin(iarg,jarg)+pattern(i,j)/4.
				enddo
			enddo
			open (13,file='c3pat_2003_296_2x2.grd')
			print *,' writing out c3pat_2003_296_2x2.grd'
			write(13,130)
130			format('DSAA'/'636 128'/'0 1272'/'0 256'/'-10 100')	
			do j=1,128
				write(13,15)(patternbin(i,j),i=1,636)
			enddo
			write(13,16)patdark,(infilesv(i),i=1,10)
			close(13)
			endif
		endif
		go to 98
96	continue
c
	if(k.ge.1)go to 97	!$$$$$	bail quick	
c
	pixsum=0.
	ipixsum=0
	pixdif=0.
	ipixdif=0
	patcorr=dark/patdark
	do j=1,256
	do i=1,1272
		iarg=int4(float(iframe(i,j)))
		if(iarg.lt.0)iarg=iarg+65536
		if(iarg.eq.65535)isat=isat+1
		arg=float(iarg)	!/0.75
		arg=arg-ped
		if(i.ge.17.and.i.le.1268)arg=arg-pattern(i,j)*patcorr	!scale pattern by dark current
		frame(i,j)=arg
		if(iarg.eq.0)frame(i,j)=0.
	enddo
	enddo
c	
	call INDEXX(n,frame,indx)
	call measles(n,frame,indx,meascnt,ineigh,ihist,ikill)
	do j=1,256
		do i=1,1272
			if(ikill(i,j).eq.1)frame(i,j)=0.	!KILL measles pixels
			if(i.gt.22.and.i.lt.1263.and.j.gt.1)then
			    if(frame(i,j).ne.0..and.frame(i,j).lt.5000.)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.22.and.i.lt.1263)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.gt.21.and.i.lt.1263)then
					ipixdif=ipixdif+1
					pixdif=pixdif+abs(frame(i,j-1)-frame(i,j))
				endif
			    endif
			endif
		enddo	
	enddo
c
	pixsum=pixsum/float(ipixsum)
	pixdif=pixdif/float(ipixdif)
	if(isat.gt.20)print 12,dddate,h,m,s,ped-pedinit,dark,meascnt,pixsum,pixdif,isat,k-10
	if(isat.gt.20)go to 98
c	
c97	if(igraph.eq.1)write(12,120)dddate,h,m,s,pedright,pediff
97	if(igraph.eq.1)write(12,120)dddate,h,m,s,ped-pedinit,dark,meascnt,pixsum,pixdif,exp((camtemp+40.)/8.7),isat
c	if(igraph.eq.1)write(21,33)outname,ddoy,ped,dark,q,meascnt,pixsum,pixdif
c	if(mod(k,10).eq.0)then
c	print 12,dddate,h,m,s,pedright,pediff,k-10
	print 12,dddate,h,m,s,ped-pedinit,dark,meascnt,pixsum,pixdif,isat,k-10
12	format('"',a8,1x,i2.2,':',i2.2,':',i2.2,'"',2x,2f10.3,i6,2f10.2,2i7)
120	format('"',a8,1x,i2.2,':',i2.2,':',i2.2,'"',2x,2f10.3,i6,3f10.2,i7)
33	format(a25,a3,2f10.3,4f13.9,i6,2f10.2)
c	endif
c	print*,k,'ped = ',ped, 'dark = ',dark,'measles = ',meascnt	

98	continue
	enddo
	close(10)
	if(igraph.eq.1)close(12)
c	if(igraph.eq.1)close(21)
	print *,'Start writing output files'
	close(20)	
	print *,'reached end'
	print *,isunmoon,' discarded as sun/moon based on raw/crazy pedestal'
	print *,ptrash,',',dtrash,' others ped, dark discarded out of ',kmax-kmin+1
	if(igraph.eq.0)stop	!#################################################
	print *,'It aint done until the fat lady sings...'
	open(10,file='darkdist36all.grd')
	write(10,1300)
1300	format('DSAA'/'4000 2044'/'0 4000'/'0 2044'/'-10 100')
	do j=1,2044
	do i=1,4000
	if(darkdist(j,i).gt.999.99)darkdist(j,i)=999.99
	if(darkdist(j,i).lt.-99.99)darkdist(j,i)=-99.99
	enddo
	write(10,1301)(darkdist(j,i),i=1,4000)	
1301	format(20f7.2)	
	enddo
	close(10)
	print *,'Now its done!'
	stop 
	end
c********************************************************************************************
	subroutine measles(n,frame,indx,meascnt,ineigh,ihist,ikill)	!this version for 1x1 binning
	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(325632),ikill(1272,256)
	real*4 frame(1272,256),z(3,3),z0,z0sv
	real*4 dcut,sumzsq,devcut
	
	data dcut /27.76/				!55.52 is a 1000-electron cut: we have 500 here.
	data devcut /25./				!5 standard deviation cut, since this is sigma**2
	
	meascnt=0

	do m=1,n
		iarg=1+n-m
		j=(indx(iarg)-1)/1272
		i=indx(iarg)-1272*j
		j=j+1
		if(i.le.21.or.i.ge.1264.or.j.lt.2.or.j.gt.255)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(325632),N,I,J,L,IR,INDXT	!Absoft Compiler dislikes INDX(N), etc.
      real*4 ARRIN(325632),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
c*********************************************************************************
      SUBROUTINE INDEXXX(N,ARRIN,INDX)	!from Numerical Receipes...
      implicit none
      integer*4 INDX(2044),N,I,J,L,IR,INDXT	!Absoft Compiler dislikes INDX(N), etc.
      real*4 ARRIN(2044),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