c
c name:		CRX_ped_dark_cam3.for
c purpose:	determines pedestal & dark current, finds black/white measles & cosmic rays
c inputs:	*.fts data files listed in "namelistCRX.txt"
c outputs:	before and after *.grd files 
c    "		The output file xxx.dat presents information for graphs of the various quantities
c calls:	subroutines  measles (+ plane), and INDEXX 
c mods:		April 2005, 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 2x2 mode as here...
c		Another consequence of this is that columns 11 and 631 are NOT half 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 CRX_ped_dark_cam3
	implicit none
	integer*1 headfts(5760),iframe1(2,636,128),iframe2(2,636,128)
c	integer*2 head(3)
	integer*4 iframe(636,128),klast,itriangle,makemap
	integer*4 h,m,s,hh,mm,ss,i,j,k,nb,kmax,kmin,iarg,ineigh,pedrecover,trash,ipcnt(636,128)
	integer*4 idarkcnt,ipedcnt,ibadped,ibaddark,ihist(250),icrcnt,meascnt,ikill(636,128),iavgcnt,isat
	integer*4 icrneig,n,inegmeas,ipixsum,ipixdif,icrx,iavg,igraph,iffix,isunmoon,isatcnt,imed
	integer*4 indx(81408),indxx(511)		!this last one a +1: strange
	integer*4 mode,id,jd,jarg,med1,med2,iarg1,iarg2,iarg3,iarg4
	real*4 dark,darksum,ped,pedlast,pedinit,pedinit2,pedsum,ddark,trianglesum,glare
	real*4 frame(636,128),frameavg(636,128),frameout(636,128),badsum,pattern(636,128),arg,patcr,patcorr
	real*4 pixsum,pixdif
	real*4 pednaive,darksv(510),patdark,darkmed(510),darknaive
	real*4 pedsv(510),darkscale,median,dmed,pedtest
	real*4 weight,weightsum
	real*4 glaremap(636,128),glarecnt(636,128)
c	real*8 q(4),qq(4),t,t0,torbit,targ,tprev
	character*47 out2file,out1file
	character infile*45,infile1*45,infile2*45	!new file naming convention
	character nname*14				!new file naming convention
	character hhead(5)*10,dum*1
	data dum /' '/
c
	data infile    /'D:/Sean/Eclipse/Data/c3m1_2005_098_hhmmss.raw'/
	data infile1   /'D:/Sean/Eclipse/Data/c3m1_2005_098_hhmmss.raw'/
	data infile2   /'D:/Sean/Eclipse/Data/c3m1_2005_098_hhmmss.raw'/
c
	data out1file  /'D:/Sean/Eclipse/Crx/c3m1_2005_098_hhmmsscr1.grd'/	!cosmic-ray not removed file...
	data out2file  /'D:/Sean/Eclipse/Crx/c3m1_2005_098_hhmmsscrx.grd'/	!cosmic-ray removed file...
c	
	data ihist /250*0/
	data badsum /0./
	data pattern /81408*0./
	data ipcnt /81408*0/
	data n /81408/
	data trash,iffix,isunmoon,isatcnt /4*0/
	data patcr /55.52/	!55.52 adu's = 1000 electrons for 2 x 2 binning
	data weightsum /0./
	data iframe1,iframe2 /325632*0/
	data glaremap,glarecnt /162816*0./
c
	icrx=0			!icrx=1 means make before and after files
	igraph=0		!igraph=1 means make GRAPHER output file and *.txt file for indexing...
	mode = 2
	makemap=1	
	id = 1272
	jd = 256
	id=id/mode
	jd=jd/mode
c	open (13,file='D:\Sean\Patterns\c3cal_2005_061_050915_1.grd',readonly)
c	print *,'reading in D:\Sean\Patterns\c3cal_2005_061_050915_1.grd'
	open (13,file='D:\Sean\Patterns\c3cal_2005_097_063859_1.grd',readonly)
	print *,'reading in D:\Sean\Patterns\c3cal_2005_097_063859_1.grd'
	read(13,1)hhead
1	format(a10)
	do j=1,128
	read(13,15)(pattern(i,j),i=1,636)
15	format(20f9.3)
	enddo
	read(13,1500)patdark
1500	format(f10.4)
	print *,' patdark = ',patdark
	close(13)
	open (10,file='D:\Sean\Eclipse\Data\namelist3.txt',readonly)
c	if(igraph.eq.1)open (12,file='cam3_eclipse2005_CRX.dat')
	if(igraph.eq.1)open (21,file='cam3_eclipse2005_CRX.txt')
	pedinit=0.	!0.	!1070.
	pedinit2=0.
	pedrecover=0
	icrcnt=0
	iavgcnt=0
c	open(20,file='DOY098_del_glare3.dat')
c
c********************************************************************
	kmin=1
c	kmax=651
c	kmax=1304		!DOY098, 2005, the whole file
	kmax=300		!DOY098, 2005, 1st early	
c	kmin=301		!DOY098, 2005, 1st late
c	kmax=600		!DOY098, 2005, 1st late
c	kmin=1126		!DOY098, 2005, 2nd early
c	kmax=1199		!DOY098, 2005, 2nd early	!1196?
	do k=1,kmax
	if(k.eq.1)pedlast=pedinit
c
	read (10,10)nname,h,m,s
10	format(39x,a14,3i2.2)	
	if(k.lt.kmin)go to 98
	klast=k
c
	write(infile1(22:35),'(a14)')nname	
	write(infile1(36:41),'(3i2.2)')h,m,s
c	
	ss=s+36
	mm=m+41
	hh=h+1
	if(ss.ge.60)then
	ss=ss-60
	mm=mm+1
	endif
	if(mm.ge.60)then
	mm=mm-60
	hh=hh+1
	endif
	if(hh.ge.24)then
	hh=hh-24
	endif
c	
	write(infile2(22:35),'(a14)')nname	
	write(infile2(36:41),'(3i2.2)')hh,mm,ss
c	print *,k,dum,infile1,dum,infile2
	open(11,file=infile1,form='binary',readonly)
	read(11)headfts,(((iframe1(nb,i,j),nb=1,2),i=1,636),j=1,128)
	close(11)
	open(11,file=infile2,form='binary',readonly)
	read(11)headfts,(((iframe2(nb,i,j),nb=1,2),i=1,636),j=1,128)
	close(11)
c	if(k.ge.1)go to 98		!useful when making a quick run to check presence of frames
c
c
c	check for Moon nearby (may trigger on planets), by finding saturated pixels
c
	isat=0
c	do i=11,632
c	do j=1,128
c	if(iframe1(i,j).lt.0.and.iframe1(i,j).gt.-10)isat=isat+1
c	if(iframe2(i,j).lt.0.and.iframe1(i,j).gt.-10)isat=isat+1
c	enddo
c	enddo
	if(isat.ne.0)then
		isatcnt=isatcnt+1
		print *,' frame ',k,': saturated pixels...',isat
		go to 98
	endif
c
	ped=0.
	dark=0.
	idarkcnt=0
	ibaddark=0
	darksum=0.
	ipedcnt=0
	ibadped=0
	pedsum=0.
	pednaive=0.
	icrcnt=0
	icrneig=0
	meascnt=0
	ineigh=0
	inegmeas=0
	do j=1,128
	do i=1,636
	ikill(i,j)=0
	enddo
	enddo
c
c	Here load up the difference (one orbit to the next) between the two frames
c	
	do i=1,636
	do j=1,128
		iarg1=int4(iframe1(1,i,j))
		iarg1=iarg1+128				!this undoes Paul's subtraction of 32768 from the 2-byte number
		iarg2=int4(iframe1(2,i,j))
		if(iarg2.lt.0)iarg2=iarg2+256		!this fixes bit #8 in the less-significant byte
		iarg3=256*iarg1+iarg2
		iarg1=int4(iframe2(1,i,j))
		iarg1=iarg1+128				!this undoes Paul's subtraction of 32768 from the 2-byte number
		iarg2=int4(iframe2(2,i,j))
		if(iarg2.lt.0)iarg2=iarg2+256		!this fixes bit #8 in the less-significant byte
		iarg4=256*iarg1+iarg2
c
		iframe(i,j)=iarg3-iarg4			!this subtracts the frames
		iframe(i,j)=iframe(i,j)/2		!after DOY055, 2005, 1700h UT
c		if(iarg4.gt.2500.or.iarg3.gt.2500)iframe(i,j)=0
	enddo
	enddo
c	
c**********************************!	Pedestal calculation
c
	do j=1,128
	do i=7,8
	    if(iframe(i,j).ne.0)then	!the "100" here is an artifact of distinguishing zero from absent (non-FOV)
		iarg=iframe(i,j)
		arg=float(iarg)		!/0.75	!"headroom" restored: this occurs 3 more times below
		pednaive=pednaive+arg
		ipedcnt=ipedcnt+1
		pedsv(ipedcnt)=arg
	    endif
	enddo
	do i=635,636
	    if(iframe(i,j).ne.0.and.j.ne.128)then
		iarg=iframe(i,j)
		arg=float(iarg)		!/0.75
		pednaive=pednaive+arg
		ipedcnt=ipedcnt+1
		pedsv(ipedcnt)=arg
	    endif
	enddo
	enddo
	pednaive=pednaive/510.
	ipedcnt=0
	call indexxx(510,pedsv,indxx)
c	print 2000,(i,pedsv(indxx(i)),i=1,510)
2000	format(i4,1x,f10.1)
	pedtest=pedsv(indxx(255))		!here compare with the naive pedestal median
c
c	Here exclude "white measles", bad bands, and high pedestals	
c
	do j=1,510
c	    if(pedsv(indxx(j))-pedtest.gt.4.)go to 20	!exclude high pedestals above 72 electrons...before 4-05
c	    if(pedsv(indxx(j))-pedtest.gt.-4.)then
	    if(pedsv(indxx(j))-pedtest.gt.20.)go to 20	!exclude high pedestals above 180 electrons...after 4-05
	    if(pedsv(indxx(j))-pedtest.gt.-10.)then
	    	ipedcnt=ipedcnt+1
	    	pedsum=pedsum+pedsv(indxx(j))
	    endif
	enddo
c20	if(ipedcnt.lt.450)ibadped=1		!scaled from eng. data
20	if(ipedcnt.lt.320)ibadped=1		!200=DOY288: 320 would be OK for 2005)
	ped=pedlast				!pedlast is a final resort value when everyone fails
	if(ipedcnt.gt.0)ped=pedsum/float(ipedcnt)
c	print *,k,pednaive,ped,ipedcnt
c	
c**********************************!	begin dark current calculation
c
	do j=1,128
	do i=9,10
	if(iframe(i,j).ne.0)then
		iarg=iframe(i,j)
		arg=float(iarg)-ped		!/0.75	!"headroom" restored: this occurs more times below
		darksum=darksum+arg
	endif
	enddo
	do i=633,634
	if(iframe(i,j).ne.0.and.j.ne.128)then
		iarg=iframe(i,j)
		arg=float(iarg)-ped		!/0.75
		darksum=darksum+arg
	endif
	enddo
	enddo
	darknaive=darksum/510.
	darkscale=0.75*darknaive/patdark	!naive dark current yields scale factor to "patdark"
	if(darkscale.gt.2.5)darkscale=2.5
c
c	Now go back through using scale factor to filter CR hits from dark calculation
c 	
	darksum=0.
	do j=1,128
	do i=9,10
	if(iframe(i,j).ne.0)then
		iarg=iframe(i,j)
		arg=float(iarg)-ped		!/0.75	
		idarkcnt=idarkcnt+1
		darksum=darksum+arg
		jarg=i-8+4*(j-1)
		if(k.eq.1)darksv(jarg)=pattern(i,j)
		darkmed(jarg)=arg
		  ddark=arg-pattern(i,j)*darkscale
		  if(ddark.gt.patcr)then	!1000 e-
		    idarkcnt=idarkcnt-1
		    darksum=darksum-arg
		    darkmed(jarg)=10000.
		    ddark=0.
		  endif
	endif
	enddo
	enddo
	do j=1,127
	do i=633,634
	if(iframe(i,j).ne.0)then
		iarg=iframe(i,j)
		arg=float(iarg)-ped		!/0.75
		idarkcnt=idarkcnt+1
		darksum=darksum+arg
    	        jarg=i-630+4*(j-1)	!this is weird & probably wrong...
    	        if(k.eq.1)darksv(jarg)=pattern(i,j)
		darkmed(jarg)=arg
		  ddark=arg-pattern(i,j)*darkscale
		  if(ddark.gt.patcr)then
		    idarkcnt=idarkcnt-1
		    darksum=darksum-arg
		    darkmed(jarg)=10000.
		    ddark=0.
		  endif
	endif
	enddo
	enddo
	if(idarkcnt.le.10)then	!hopeless...
		ibaddark=1
		go to 53
	endif
c
c	Get median:"interpolate within the histogram..."
c
	call indexxx(510,darkmed,indxx)
	med1=0
	med2=0
	imed=ifix(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,510-imed
	if(darkmed(indxx(imed+j)).ne.dmed)then
		med2=imed+j
		go to 52
	endif
	enddo
52	continue
	median=darkmed(indxx(imed))-0.5+float(imed-med1)/float(med2-med1)
	if(idarkcnt.lt.100)ibaddark=1
	dark=darksum/float(idarkcnt)
c	print *,k,infile,ped,ipedcnt,darknaive,median,idarkcnt,dark/median
	dark=0.85*median		!ad hoc scale change in dark current for 2x2 binned data
c	dark=dark-0.75			!ad hoc shift because Mark Cooke's 2x2 binning must differ from mine...
c
c	Now correct time-plotted covered pixels by best dark current scale factor	
c
c	print *,k,pednaive,ped,ipedcnt,dark,idarkcnt
c
53	if((ibadped.eq.0.and.ibaddark.eq.0).or.k.eq.1)then
		pedlast=ped
		pedrecover=0
	else	
		trash=trash+1
		if(ibadped.eq.1)then
			print 1000,k,pednaive-pedinit,ped-pedinit,ipedcnt,dark,median,idarkcnt,infile1
1000			format('BAD PED  ',i6,', ped='2f6.2,', ipedcnt=',i3,' dark = ',2f7.2,', darkcnt=',i3,1x,a47)
			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
c				print*, 'Attempting Pedestal Recovery '
			endif
		else
			print 1001,k,pednaive-pedinit,ped-pedinit,ipedcnt,dark,median,idarkcnt,infile1
1001			format('BAD DARK ',i6,', ped='2f6.2,', ipedcnt=',i3,' dark = ',2f7.2,', darkcnt=',i3,1x,a47)
		endif
		ped=-9999.
		dark=-9999.
		go to 98
	endif
	weight=1.-float(iabs(k-129))/120.		!DOY098, 2005, 1st early
	if(k.gt.300)weight=1.-float(iabs(k-429))/120.	!DOY098, 2005, 1st late
	if(k.gt.600)weight=1.-float(iabs(k-1197))/71.	!DOY098, 2005, 2nd early
c	
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,128
	do i=1,636
		iarg=iframe(i,j)
		arg=float(iarg)	!/0.75
c		if(i.eq.11.or.i.eq.632)arg=arg*2.0	!see note at top...
		arg=arg-ped
		if(i.ge.9.and.i.le.634)arg=arg-pattern(i,j)*patcorr	
c		if(i.ge.9.and.i.le.634)arg=arg-dark*(1.+float(j-64)/100.)	!subtract an avg dark current, with a slope
		frame(i,j)=arg
		if(iarg.eq.0)frame(i,j)=0.
	enddo
	enddo
c
c	Calculate the average response in the left-hand triangle towards the Sun but beyond the FOV
c	Closely similar results when smaller triangle used, but noise increases.
c	Indication of a shift, 10% factor occurring sometime in the ratio of triangle to single corner pixel.
c	For now, this probably indicates an upper limit to systematic error here.
c
	trianglesum=0.
	itriangle=0
	do i=11,16
	do j=112+i,128
	trianglesum=trianglesum+frame(i,j)
	itriangle=itriangle+1
	enddo
	enddo
	glare=trianglesum/float(itriangle)
	glare=glare/16.65	!put close to correction scale in order to join to camera 2's map	!15.5 or 258. (16.65)?
c
	call INDEXX(n,frame,indx)
	call measles(n,frame,indx,meascnt,ineigh,ihist,ikill)
	iavgcnt=iavgcnt+1
	do j=1,128
		do i=1,636
			if(frame(i,j).gt.32768.)print *,i,j,frame(i,j)
c			if(ikill(i,j).eq.1)frame(i,j)=-100.	!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
				if(iavg.eq.1)frameavg(i,j)=frameavg(i,j)+frame(i,j)/10.	
				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.11.and.i.lt.633)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.10.and.i.lt.633)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	write(20,56)k,weight,glare,ped-pedinit,darkscale,dark,pixsum,meascnt/10,ipedcnt,idarkcnt
56	format(i6,6f10.2,3i6)
	glare=weight
	do i=1,636
	do j=1,128
c		if(iframe1(1,i,j).ne.-128)then
		if(glare.gt.0.)then
		    if(k.ge.301.and.k.le.600)then
c			if(frame(i,j).lt.15..and.frame(i,j).gt.-75.)then	!75 for DOY151, 15 for DOY327
			if(frame(i,j).lt.10..and.frame(i,j).gt.-100.)then
				glaremap(i,j)=glaremap(i,j)-frame(i,j)*glare
				glarecnt(i,j)=glarecnt(i,j)+glare
			endif
		    else
c			if(frame(i,j).gt.-15..and.frame(i,j).lt.75.)then
			if(frame(i,j).gt.-10..and.frame(i,j).lt.100.)then
				glaremap(i,j)=glaremap(i,j)+frame(i,j)*glare
				glarecnt(i,j)=glarecnt(i,j)+glare
			endif
		    endif
		endif
c		endif
	enddo
	enddo
c	
97	continue
	if(igraph.eq.1)write(21,33)nname,h,m,s,ped,dark,meascnt,pixsum,pixdif,glare
	if(mod(k,10).eq.0)then
	print 12,nname,h,m,s,ped-pedinit,dark,meascnt,pixsum,pixdif,glare,k
12	format(a14,3i2.2,2x,2f10.3,i6,3f10.2,i7)
33	format(a14,3i2.2,2x,2f10.3,i6,3f10.2)
	endif
	write(out1file(31:40),'(a10)')infile1(24:33)
	write(out2file(31:40),'(a10)')infile1(24:33)
	if(iavg.eq.1.and.iavgcnt.eq.10)then
c	open (16,file=out3file)
c	write(16,140)
c	write(16,150)((frameavg(i,j),i=1,636),j=1,128)
c	close(16)
	iavgcnt=0
	do j=1,128
	do i=1,636
	frameavg(i,j)=0.
	enddo
	enddo
	endif	
c	if(icrx.eq.1)open (14,file=out1file)
	if(icrx.eq.1)open (15,file=out2file)
c	if(icrx.eq.1)write(14,140)
	if(icrx.eq.1)write(15,140)
140	format('DSAA'/'636 128'/'0 636'/'0 128'/'-50 50')
	if(icrx.eq.1)write(15,150)((frameout(i,j),i=1,636),j=1,128)
150	format(20f9.2)
c    	if(icrx.eq.1)close(14)
	if(icrx.eq.1)close(15)
c	print*,k,'ped = ',ped, 'dark = ',dark,'measles = ',meascnt	

98	continue
	enddo
	close(10)
	if(igraph.eq.1)close(12)
	if(igraph.eq.1)close(21)
c	
	close(20)	
	print *,'reached end'
	print *,isatcnt,' discarded as sun/moon based on saturated pixel(s)'
	print *,isunmoon,' more discarded as sun/moon based on raw/crazy pedestal'
	print *,trash,' others ped/dark discarded out of ',kmax-kmin+1
c
c	form, write output glare map
c
	if(makemap.eq.1)then
	open(10,file='D:\Sean\Eclipse\Glare3_2005_1st_e.grd')
	write(10,100)
100	format('DSAA'/'636 128'/'0 636'/'0 128'/'-1 30')
	do j=1,128
	do i=1,636
		if(glarecnt(i,j).ne.0.)glaremap(i,j)=glaremap(i,j)/glarecnt(i,j)
	enddo
	enddo
	write(10,110)((glaremap(i,j),i=1,636),j=1,128)
110	format(20f9.2)	
	close(10)
	endif
	print *,'done'
	stop 
	end
c********************************************************************************************
	subroutine measles(n,frame,indx,meascnt,ineigh,ihist,ikill)	!this version for 2x2 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	**********************************note this version rejects both positive and negative measles...
c	to accomplish this, we need to sort by abs(response) rather than the usual just-response...
c	Also note if used, the "if below threshold" needs to be changed from 99 to 98, or the negatives get missed.
c 		
	implicit none
	integer*4 i,j,ii,jj,meascnt,neigh,ineigh,nfix(3,3),ihist(250),iarg,n,m
	integer*4 indx(81408),ikill(636,128)
	real*4 frame(636,128),z(3,3),z0,z0sv
	real*4 dcut,sumzsq,devcut
	
c	data dcut /55.52/				!55.52 is a 1000-electron cut for 2x2 binning
	data dcut /10./	
c	data devcut /25./				!5 standard deviation cut, since this is sigma**2
	data devcut /9./
	
	meascnt=0

	do m=1,n
		iarg=1+n-m
		j=(indx(iarg)-1)/636
		i=indx(iarg)-636*j
		j=j+1
		if(i.le.12.or.i.ge.631.or.j.lt.2.or.j.gt.127)go to 98	!bail if near edges
c		if(frame(i,j).ne.0.and.frame(i,j).ne.-10..and.abs(frame(i,j)).lt.dcut)go to 98	!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(abs(z(2,2)).lt.dcut)go to 98		!bail if not above threshold
		do jj=1,3
		do ii=1,3
			if(z(2,2).gt.0..and.z(ii,jj).gt.z(2,2))go to 98
			if(z(2,2).lt.0..and.z(ii,jj).lt.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.abs(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 N
      integer*4 INDX(N),I,J,L,IR,INDXT	!Absoft Compiler dislikes INDX(N), etc.
      real*4 ARRIN(N),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(510),N,I,J,L,IR,INDXT	!Absoft Compiler dislikes INDX(N), etc.
      real*4 ARRIN(510),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