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),trail2a(94),trail2b(177),tdum1(12),trail3(21),tdum2(28),trail4(69)
	integer*2 iframe(636,128),head(3), framekill
	integer*4 h,m,s,i,j,k,kmax,kmin,ii,jj,iarg,pedrecover,trash,idarkcnt,ipedcnt,ibadped,iargsave
	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,shutcount,akill
	real*4 ped1,ped2,pedright,pediff
	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),frame2(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,camtemp,camtempsave
	real*8 q(4),t,tprev,t0,torbit,targ
	character infile*52
	character doy(15)*3,ddoy*3
	character nname*15,outname*26
	character hhead(5)*10
	character ddate(15)*8,dddate*8
	data infile   /'../../zone/frames/cam3_DOY/c3frm_2003_DOY_hhmmss.nic'/
	data outname  /'c3pattern_296_orbit###.grd'/
	data ddate /'10/23/03','10/24/03','10/25/03','10/26/03','10/27/03','10/28/03','10/29/03','10/30/03','10/31/03',
     &	'11/01/03','11/02/03','11/03/03','11/04/03','11/05/03','11/06/03'/
	data doy /'296','297','298','299','300','301','302','303','304','305','306','307','308','309','310'/
	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.130D0/	!value determined from quaternion analysis (see orbitcompare.for & badquarts.for)
	data iarglast /1524/
	data framemin /81408*66000./
	data frame2   /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
	shutcount=0
	akill=0
	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.
	open (13,file='c3pat_2003_296_2x2.grd',readonly)
	print *,'reading in c3pat_2003_296_2x2.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
	read(13,159)patdark
159	format(f10.4)
	print*, 'Patdark: ',patdark
	close(13)
c***
	open (10,file='cam3_list_296to310.txt',readonly)
c
	pedinit=1046.
	pedinit2=1051.
	pedrecover=0
	framekill=0
c********************************************************************
cc	kmax =  275061	!275252	!275061 marks end of orbit #183, 275255 is total # on list
	kmax =  277152  !marks the end of orbit #209 !!278230 is the last frame in the list...
	k1=16525
	k2=37677
	k3=56249
	k4=76899
	k5=98174
	k6=118242
	k7=137929
	k8=138059
	k9=156153
	k10=177001
	k11=196303
	k12=217387
	k13=238380
	k14=259228
	kmin=1
c***	!258887 is the last frame of orbit 195, following orbits have 1x1 and shutter closed frames for c3frm_2003_309/310
	!I think the program needs the 1st frame of the next orbit to print out the preceeding orbit's pattern
	kmax=258888
	t0=18780.6D0

	do k=1,kmax
	if(k.eq.1)pedlast=pedinit
c
	read (10,10)nname,h,m,s
10	format(a15,3i2.2)
	if(k.le.4)go to 98	!skip first two frames, an artifact from early pattern days...
	if(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
c	if(targ.lt.0.D0)then
c		targ=targ+1.
c		iorbit=iorbit-1
c	endif
c	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(20:22),'(i3.3)')iorbit-1
			print *,'writing first map = ',outname
			open (20,file=outname)
			open (201,file='partial_map.grd')
			write(20,20)
			write(201,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)
				write(201,21)(frame2(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)
			close(201)
		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(20:22),'(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
	iargsave=iarg
	if(k.gt.k6.and.iarg.ge.975.and.iarg.le.1325)then	!bail out after the blast, [1000,1300] would probably work...
		ibadorbit=ibadorbit+1
		go to 98
	endif

	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(24:26),'(a3)')ddoy
	write(infile(39:41),'(a3)')ddoy
	write(infile(43:48),'(3i2.2)')h,m,s
c	if(ddoy.eq.'145'.and.h.eq.13.and.m.eq.44.and.s.eq.07)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'145'.and.h.eq.13.and.m.eq.44.and.s.eq.07)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.28.and.s.eq.19)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.28.and.s.eq.19)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.39.and.s.eq.27)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.39.and.s.eq.27)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.42.and.s.eq.23)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.42.and.s.eq.23)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.51.and.s.eq.11)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.51.and.s.eq.11)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.52.and.s.eq.39)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.52.and.s.eq.39)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.52.and.s.eq.51)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.52.and.s.eq.51)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.54.and.s.eq.43)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.54.and.s.eq.43)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.56.and.s.eq.23)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.56.and.s.eq.23)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.58.and.s.eq.19)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.58.and.s.eq.19)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.15)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.15)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.19)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.19)go to 98
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.47)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.13.and.m.eq.59.and.s.eq.47)go to 98
c	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.00.and.s.eq.15)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.00.and.s.eq.15)go to 98
c	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.03)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.03)go to 98
c	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.07)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.07)go to 98
c	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.15)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.15)go to 98
c	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.39)print *,infile,': Killed, k =',k
c	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.02.and.s.eq.39)go to 98
	if(mod(k,500).eq.0)  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,trail2a,camtemp,trail2b,tdum1,trail3,tdum2,trail4
	close(11)

c*******!Added by Aaron, not sure if this is needed since it also in the CRX program**********
	!Check for shutter closed data
	if(trail1(15).eq.0)then
	  shutcount=shutcount+1
	  if(k.eq.1) tprev=t-100.0
	  if(camtemp.gt.0.)then
c	    print*, 'Warning, bad camtemp:',camtemp,' for:  ',infile
	    camtemp=camtempsave
	  elseif((t-tprev).gt.90.) then
	      framekill=1
c	      print*, 'Killed Shutter Closed Frame:  ',infile
	      akill=akill+1
	      go to 98
	  endif
	elseif(framekill.eq.1)then
	  shutcount=shutcount+1
	  framekill=0
	  akill=akill+1
c	  print*, 'Killed Shutter Open Frame:    ',infile
	  go to 98
	endif
	framekill=0
	camtempsave=camtemp
	tprev=t
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
		arg=arg-pattern(i,j)
		pednaive=pednaive+arg
		ped1=ped1+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
		arg=arg-pattern(i,j)
		arg=arg-pedright
		pednaive=pednaive+arg
		ped2=ped2+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
	
	ped1=ped1/256.
	ped2=ped2/254.
	pediff=ped2-ped1
	if(abs(pediff).lt.4)  pedright=pedright+0.1*pediff

	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"
c	print *,k,darknaive,darkscale
	if(darkscale.gt.1.0)darkscale=1.0
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.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.
		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.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.
		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		!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(nmap.eq.1.and.iargsave.lt.975.or.iargsave.gt.1325)then
		  if(frame(i,j).lt.frame2(i,j)) frame2(i,j)=frame(i,j)
		endif
		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, ?? on the list were hand killed'
	print *,ibadorbit,' discarded due to bad-sun orbital location between (cutout)/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*, akill,' killed by Aaron'
	print*, shutcount,'shutter closed frames'
	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