c
c name:		Pattern_Get_cam3.for
c purpose:	determines pedestal & dark current, finds Least Pattern for each orbit
c inputs:	*.nic data files listed in "namelistCRX.txt"
c outputs:	min pattern *.grd files 
c calls:	subroutine INDEXX , PLANE
c mods:		March 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 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 Pattern_Get_cam3
	implicit none
	integer*1 trail1(75),trail2(275),tdum1(12),trail3(21),tdum2(28),trail4(69)
	integer*2 iframe(636,128),head(3)
	integer*4 h,m,s,i,j,k,kmax,kmin,ii,jj,iarg,pedrecover,trash,idarkcnt,ipedcnt,ibadped
	integer*4 iffix,isunmoon,ibadorbit,imed,nmap,iorbitkeep,ineigh,nmin,icount
	integer*4 k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,indxx(511)
	integer*4 mode,id,jd,jarg,iorbit,iarglast,med1,med2,isumall,iless,iisv,jjsv
	real*4 dark,darksum,ped,pedlast,pedinit,pedinit2,pedsum,ddark,sumall,z(3,3),z0,dcut
	real*4 frame(636,128),framemin(636,128),frame1(636,128),pattern(636,128),arg,patcr
	real*4 pednaive,patdark,patcorr,darknaive,crazyped1,crazyped2,frameout(636,128)
	real*4 darkmed(510),pedsv(510),darksv(510),darkscale,median,dmed,pedtest
	real*8 q(4),t,t0,torbit,targ
	character infile*43
	character doy(15)*3,ddoy*3  
	character nname*15,outname*22
	character hhead(5)*10
	character ddate(15)*8,dddate*8
	data infile   /'D:/Buffer3/NicDOY/c3frm_2003_DOY_hhmmss.nic'/
	data outname  /'c3pattern_orbit###.grd'/		
	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 trash,iffix,isunmoon,ibadorbit /4*0/
	data patcr /11.10/	!55.52 adu's = 1000 electrons for 2 x 2 binning
	data torbit /6096.297D0/	!value determined from quaternion analysis (see orbitcompare.for & badquarts.for)
	data iarglast /1524/
	data patdark /8.68373/	!value from eng. data DOY144, last orbit
	data framemin /81408*66000./
	data nmap /0/
	data iorbitkeep / 0 /	!0 means do them all
	data nmin / 0 /		!0 means do them all starting at the beginning
	data dcut /10./		!defines a minimum deviation for a pixel (one of eight) to be a discrepant neighbor
	mode = 2	
	id = 1272
	jd = 256
	id=id/mode
	jd=jd/mode
c
c***		!this section for reading in pattern rather than making it here.
	patdark=8.18373
	open (13,file='pattern3l_144.grd',readonly)
	print *,'reading in pattern3l_144.grd'
	read(13,1)hhead
1	format(a10)	
	do j=1,128
	read(13,13)(pattern(i,j),i=1,636)
13	format(20f9.3)
	enddo
	close(13)
c***	
c	
	open (10,file='namelist3.txt',readonly)
c
	pedinit=1046.
	pedinit2=1051.
	pedrecover=0
c********************************************************************
	kmax =  275061	!275252	!275061 marks end of orbit #183, 275255 is total # on list
	k1=10289
	k2=31657
	k3=53025
	k4=74355
	k5=95670
	k6=116941
	k7=138215
	k8=159486
	k9=180769
	k10=202008
	k11=223281
	k12=244555
	k13=265882
	k14=275234	!19 frames included in DOY 158 (after major calibration)
	kmin=1
	t0=+2014.5D0	!DOY144
c	kmax=6100
	do k=1,kmax
	if(k.eq.3)pedlast=pedinit
c
	read (10,10)nname,h,m,s
10	format(39x,a15,3i2.2)
	if(k.le.2)go to 98	!skip first two frames, an artifact from early pattern days...
	if(k.gt.2.and.k.lt.kmin)go to 98
c	
	t=3600.D0*dfloat(h)+60.D0*dfloat(m)+dfloat(s)
	if(k.ge.k1)t=t+86400.D0
	if(k.ge.k2)t=t+86400.D0
	if(k.ge.k3)t=t+86400.D0
	if(k.ge.k4)t=t+86400.D0
	if(k.ge.k5)t=t+86400.D0
	if(k.ge.k6)t=t+86400.D0
	if(k.ge.k7)t=t+86400.D0
	if(k.ge.k8)t=t+86400.D0
	if(k.ge.k9)t=t+86400.D0
	if(k.ge.k10)t=t+86400.D0
	if(k.ge.k11)t=t+86400.D0
	if(k.ge.k12)t=t+86400.D0
	if(k.ge.k13)t=t+86400.D0
	if(k.ge.k14)t=t+86400.D0
c
c***	Here bail if the original (when forming the mingrit array) grittiness > 40 indicates bad contamination
c	Note this section must be commented out when forming the mingrit and min response arrays
c	
	targ=(t-t0)/torbit
	iorbit=int(targ)+1
	if(targ.lt.0.D0)then
		targ=targ+1.
		iorbit=iorbit-1
	endif
	iorbit=iorbit-7	!make numbering start at 1 on noon, DOY144
	targ=dmod(targ,1.0D0)
	iarg=1+int(1524.*targ)
	if(iarg.gt.1524)iarg=1524
	if(iarg.lt.iarglast)then
		print *,'Beginning orbit #',iorbit,', at k = ',k,', ',nname,h,m,s
		nmap=nmap+1
		if(nmap.eq.2)then
			write(outname(16:18),'(i3.3)')iorbit-1
			print *,'writing first map = ',outname
			open (20,file=outname)
			write(20,20)
20			format('DSAA'/'636 128'/'0 636'/'0 128'/'-10 100')			
			do j=1,128
				write(20,21)(framemin(i,j),i=1,636)
				do i=1,636
c					if(framemin(i,j).lt.-10.or.framemin(i,j).gt.200.)print *,i,j,framemin(i,j)
				enddo
21				format(20f9.3)
				do i=1,636
					frame1(i,j)=framemin(i,j)
					framemin(i,j)=66000.
				enddo			
			enddo
			close(20)
		endif
		if(nmap.gt.2.and.(iorbitkeep.eq.0.or.iorbit-1.eq.iorbitkeep))then
c		
			do j=1,128
			do i=1,636
				frameout(i,j)=0.
			enddo
			enddo
			icount=0
c			
c	Here use the plane program to find deviations of a pixel from its immediate neighbors, call that the "pattern"
c
			do j=2,127
			do i=2,635
			do ii=1,3
			do jj=1,3
				if(framemin(i-2+ii,j-2+jj).eq.66000..or.frame1(i-2+ii,j-2+jj).eq.66000.)go to 23 !bail if 3x3 has any empty pixel
				z(ii,jj)=framemin(i-2+ii,j-2+jj)-frame1(i-2+ii,j-2+jj)	!always compare with the first map
			enddo
			enddo
			ineigh=0
			do ii=1,3
			do jj=1,3
			if(z(ii,jj).gt.z(2,2))ineigh=ineigh-1
			if(z(ii,jj).lt.z(2,2))ineigh=ineigh+1
			enddo
			enddo
			if(iabs(ineigh).ne.8)go to 23	!here insist that it be biggest/smallest of the nine
			call plane(z,z0)
c
c	Here take a shot at improving a pattern value when just one screwy neighbor is present
c			
			ineigh=0
			do ii=1,3
			do jj=1,3
			if(abs(z(ii,jj)).gt.dcut.and.abs(z(ii,jj)/z(2,2)).gt.0.1)then
				ineigh=ineigh+1
				iisv=ii
				jjsv=jj
			endif
			enddo
			enddo
			if(ineigh.gt.1)go to 23
			if(ineigh.eq.1)z(2,2)=z(2,2)-z(iisv,jjsv)/8.
			frameout(i,j)=z(2,2)
			if(abs(z(2,2)).gt.10.)icount=icount+1
23			continue
			enddo	
			enddo
c					
			write(outname(16:18),'(i3.3)')iorbit-1
			print *,'writing subsequent map = ',outname,': ',icount,' >10 pattern pixels found'
			open (20,file=outname)
			write(20,22)
22			format('DSAA'/'636 128'/'0 636'/'0 128'/'-25 25')			
			do j=1,128
c				write(20,21)((framemin(i,j)-frame1(i,j)),i=1,636)
				write(20,21)(frameout(i,j),i=1,636)
				do i=1,636
					framemin(i,j)=66000.
				enddo			
			enddo
			close(20)
		endif
	endif
	iarglast=iarg
	if(iorbitkeep.ne.0.and.iorbit.ne.1.and.iorbit.ne.iorbitkeep)go to 98	!bail except for selected orbit
	if(nmap.gt.1.and.nmin.ne.0.and.nmap.lt.nmin)go to 98			!quick bail when limited range wanted
c
c	Here bail out (don't even read the frame) if orbital position is known to lie within bad sun
c	A cut on grittinexx at 40 was imposed when forming the "min" arrays: this here is at 866
c	The corresponding min value for avg response is 690, both on the trailing edge after a sun transit
c	
	if(k.gt.2.and.iarg.ge.283.and.iarg.le.866)then	!bail out
		ibadorbit=ibadorbit+1
		go to 98
	endif
c***	
c	
	dddate=ddate(1)
	ddoy=doy(1)
	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)
c
	write(infile(15:17),'(a3)')ddoy
	write(infile(19:33),'(a15)')nname
	write(infile(34:39),'(3i2.2)')h,m,s
	if(ddoy.eq.'145'.and.h.eq.13.and.m.eq.44.and.s.eq.07)print *,infile,': Killed, k =',k
	if(ddoy.eq.'145'.and.h.eq.13.and.m.eq.44.and.s.eq.07)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.28.and.s.eq.19)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.28.and.s.eq.19)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.39.and.s.eq.27)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.39.and.s.eq.27)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.42.and.s.eq.23)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.42.and.s.eq.23)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.51.and.s.eq.11)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.51.and.s.eq.11)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.52.and.s.eq.39)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.52.and.s.eq.39)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.52.and.s.eq.51)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.52.and.s.eq.51)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.54.and.s.eq.43)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.54.and.s.eq.43)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.56.and.s.eq.23)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.56.and.s.eq.23)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.58.and.s.eq.19)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.58.and.s.eq.19)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.15)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.15)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.19)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.19)go to 98
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.47)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.47)go to 98
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.00.and.s.eq.15)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.00.and.s.eq.15)go to 98
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.03)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.03)go to 98
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.07)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.07)go to 98
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.15)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.15)go to 98
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.39)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.39)go to 98
	print *,k,iorbit,infile
	open(11,file=infile,form='binary',readonly)
	read(11)head,((iframe(i,j),i=1,636),j=1,128),trail1,q,trail2,tdum1,trail3,tdum2,trail4
	close(11)
c
	ped=0.
	dark=0.
	idarkcnt=0
	darksum=0.
	ipedcnt=0
	ibadped=0
	pedsum=0.
	pednaive=0.
c	
c**********************************!	Pedestal calculation
c
	do j=1,128
	do i=7,8
	    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
		pednaive=pednaive+arg
		ipedcnt=ipedcnt+1
		pedsv(ipedcnt)=arg
	    endif
	enddo
	enddo
	do j=1,127
	do i=635,636
	    if(iframe(i,j).ne.0)then
		iarg=int4(iframe(i,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)		!/0.75
		pednaive=pednaive+arg
		ipedcnt=ipedcnt+1
		pedsv(ipedcnt)=arg
	    endif
	enddo
	enddo
	pednaive=pednaive/510.
	if(pednaive-pedinit.gt.50.)then	!This number seems to work OK for late May...
		isunmoon=isunmoon+1
		go to 98
	endif
	ipedcnt=0
	call indexxx(510,pedsv,indxx)
	pedtest=pedsv(indxx(255))		!here compare with the naive pedestal median
c	
c	Here exclude "white measles" and high pedestals above 72 electrons...
c	
	do j=1,510
	    if(pedsv(indxx(j))-pedtest.gt.4.)go to 32
	    if(pedsv(indxx(j))-pedtest.gt.-10.)then
	    	ipedcnt=ipedcnt+1
	    	pedsum=pedsum+pedsv(indxx(j))
	    endif	    
	enddo
32	if(ipedcnt.lt.450)ibadped=1		!scaled from eng. data
	if(ibadped.eq.0.and.pedsv(indxx(1))-pedtest.lt.-10..and.pedsv(indxx(2))-pedtest.lt.-10.)then
		crazyped1=pedsv(indxx(1))
		crazyped2=pedsv(indxx(2))
		print*,infile,': BAD PEDESTAL VALUES = ',crazyped1,crazyped2,pedtest
		ibadped=2
	endif
	ped=pedlast				!pedlast is a final resort value when everyone fails
	if(ipedcnt.gt.0)ped=pedsum/float(ipedcnt)
	if(ibadped.eq.0.or.k.le.3)then
		pedlast=ped
		pedrecover=0
	else	
		trash=trash+1
c		if(ibadped.eq.1)print*,infile,': BAD PEDESTAL, peds = ',ped,', ',pedtest,', ipedcnt= ',ipedcnt
		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
		ped=-9999.
		dark=-9999.
		go to 98
	endif
c	
c**********************************!	begin dark current calculation
c
	do j=1,128
	do i=9,10
	if(iframe(i,j).ne.0)then
		iarg=int4(iframe(i,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)-ped		!/0.75	!"headroom" restored: this occurs more times below
		darksum=darksum+arg
	endif
	enddo
	enddo
	do j=1,127
	do i=633,634
	if(iframe(i,j).ne.0)then
		iarg=int4(iframe(i,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)-ped		!/0.75
		darksum=darksum+arg
	endif
	enddo
	enddo
	darknaive=darksum/510.
	darkscale=0.5*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=int4(iframe(i,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)-ped		!/0.75	
		idarkcnt=idarkcnt+1
		darksum=darksum+arg
		jarg=i-8+4*(j-1)
		if(k.eq.3)darksv(jarg)=pattern(i,j)
		darkmed(jarg)=arg
		if(k.ge.3)then
		  ddark=arg-pattern(i,j)*darkscale
		  if(ddark.gt.patcr)then	!200 e-
		    idarkcnt=idarkcnt-1
		    darksum=darksum-arg
		    darkmed(jarg)=10000.
		  endif
		endif
	endif
	enddo
	enddo
	do j=1,127
	do i=633,634
	if(iframe(i,j).ne.0)then
		iarg=int4(iframe(i,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)-ped		!/0.75
		idarkcnt=idarkcnt+1
		darksum=darksum+arg
    	        jarg=i-630+4*(j-1)
    	        if(k.eq.3)darksv(jarg)=pattern(i,j)
		darkmed(jarg)=arg
		if(k.ge.3)then
		  ddark=arg-pattern(i,j)*darkscale
		  if(ddark.gt.patcr)then
		    idarkcnt=idarkcnt-1
		    darksum=darksum-arg
		    darkmed(jarg)=10000.
		  endif
		endif
	endif
	enddo
	enddo
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)
	dark=darksum/float(idarkcnt)
c	print *,k,infile,ped,ipedcnt,darknaive,median,idarkcnt,dark/median
	dark=0.76*median		!ad hoc scale change in dark current for 2x2 binned data
	dark=dark-0.75			!ad hoc shift because Mark Cooke's 2x2 binning must differ from mine...
c
c************************************************************
	if(k.le.2)go to 98
96	continue
c	
c	if(k.ge.1)go to 97	!$$$$$	bail quick	
c
	sumall=0.
	isumall=0
	patcorr=(dark/patdark)
	do j=1,128
	iless=0
	do i=1,636
		iarg=int4(float(iframe(i,j)))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)	!/0.75
c		if(i.eq.11.or.i.eq.632)arg=arg*2.0	!see note at top...
		if(iarg.ne.0)arg=arg-ped
		if(iarg.ne.0.and.i.ge.9.and.i.le.634)arg=arg-pattern(i,j)*patcorr
		frame(i,j)=arg
		if(iarg.ne.0)sumall=sumall+arg
		if(iarg.ne.0)isumall=isumall+1
		if(iarg.eq.0)frame(i,j)=66000.
		if(frame(i,j).lt.framemin(i,j))framemin(i,j)=frame(i,j)
		if(iarg.ne.0.and.iorbit.gt.1.and.frame(i,j)-frame1(i,j).lt.-20.)iless=iless+1
	enddo
	if(iless.ge.20)print *,infile,': ***WARNING, ',iless,' lower pixels (< 20) in row ',j
	enddo
c	
97	continue
98	continue
	enddo
	print *,'reached end, 18 on the list were hand killed'
	print *,ibadorbit,' discarded due to bad-sun orbital location between (283&866)/1524, inclusive'
	print *,isunmoon,' more discarded as sun/moon based on raw/crazy pedestal'
	print *,trash,' others ped/dark discarded out of ',kmax-kmin+1
	print *,'done'
	stop 
	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 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