c
c name:		Pattern_Get_cam2.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-
	program Pattern_Get_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,ii,jj,iarg,pedrecover,trash,idarkcnt,ipedcnt,ibadped
	integer*4 iffix,isunmoon,imed,nmap,iorbitkeep,ineigh,nmin,icount,isat,isatcnt
	integer*4 k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,indxx(127),ibaddark
	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(318,64),framemin(318,64),frame1(318,64),pattern(318,64),arg,patcr
	real*4 pednaive,patdark,patcorr,darknaive,crazyped1,crazyped2,frameout(318,64)
	real*4 darkmed(127),pedsv(127),darksv(127),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:/Buffer2/NicDOY/c2frm_2003_DOY_hhmmss.nic'/
	data outname  /'c2pattern_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,isatcnt /4*0/
	data patcr /6.94/	!13.89 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 framemin /20352*66000./
	data frame1 /20352*0./
	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 = 4	
	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='pattern2h_144.grd',readonly)
	print *,'reading in pattern2h_144.grd'
	read(13,1)hhead
1	format(a10)	
	do j=1,64
	read(13,13)(pattern(i,j),i=1,318)
13	format(20f9.3)
	enddo
	close(13)
c***	
c	
	open (10,file='namelist2.txt',readonly)
c
	pedinit=1016.
	pedinit2=1020.
	pedlast=pedinit
	pedrecover=0
c********************************************************************
	kmax =  270258	!No DOY158 on this file presently, skip out at orbit184, but need 1st one to get map
	k1=10087
	k2=31042
	k3=52012
	k4=72962
	k5=93902
	k6=114852
	k7=135736
	k8=156661
	k9=177590
	k10=198503
	k11=219415
	k12=240332
	k13=261242
	k14=290000	!no DOY158 is present	
	kmin=1
c	kmax=148300
	t0=+2014.5D0	!DOY144
	do k=1,kmax
	read (10,10)nname,h,m,s
10	format(39x,a15,3i2.2)
c	if(k.lt.kmin)go to 98	!for kmin=kmax, one frame treatment
	if(k.lt.kmin.and.nmap.ge.2)go to 98
	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
	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'/'318 64'/'0 318'/'0 64'/'-10 100')			
			do j=1,64
				write(20,21)(framemin(i,j),i=1,318)
c				do i=1,318
c					if(framemin(i,j).lt.-10.or.framemin(i,j).gt.200.)print *,i,j,framemin(i,j)
c				enddo
21				format(20f10.3)
				do i=1,318
					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,64
			do i=1,318
				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,63
			do i=2,317
			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.5.)icount=icount+1
23			continue
			enddo	
			enddo
c					
			write(outname(16:18),'(i3.3)')iorbit-1
			print *,'writing previous orbits map = ',outname,': ',icount,' >5 ADU pattern pixels found'
			open (20,file=outname)
			write(20,22)
22			format('DSAA'/'318 64'/'0 318'/'0 64'/'-25 25')			
			do j=1,64
c				write(20,24)((framemin(i,j)-frame1(i,j)),i=1,318)
				write(20,24)(frameout(i,j),i=1,318)
24				format(20f9.3)				
				do i=1,318
					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***	
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
c
c	This particular frame is out of sequence: found because of discrepant dark current
c	
	if(ddoy.eq.'156'.and.h.eq.08.and.m.eq.30.and.s.eq.55)print *,infile,': Killed, k =',k
	if(ddoy.eq.'156'.and.h.eq.08.and.m.eq.30.and.s.eq.55)go to 98
c
c	Frames with stripes not extending to ped/dark region
c	
	if(ddoy.eq.'145'.and.h.eq.13.and.m.eq.44.and.s.eq.39)print *,infile,': Killed, k =',k
	if(ddoy.eq.'145'.and.h.eq.13.and.m.eq.44.and.s.eq.39)go to 98
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.05.and.s.eq.31)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.05.and.s.eq.31)go to 98
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.08.and.s.eq.11)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.08.and.s.eq.11)go to 98
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.13.and.s.eq.07)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.13.and.s.eq.07)go to 98
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.16.and.s.eq.15)print *,infile,': Killed, k =',k
	if(ddoy.eq.'151'.and.h.eq.14.and.m.eq.16.and.s.eq.15)go to 98
c	if(k.gt.kmin)print *,k,iorbit,infile
	open(11,file=infile,form='binary',readonly)
	read(11)head,((iframe(i,j),i=1,318),j=1,64),trail1,q,trail2,tdum1,trail3,tdum2,trail4
	close(11)
c
c	check for Moon nearby (may trigger on planets), by finding saturated 4x4 binned pixels
c
	isat=0
	do i=6,316
	do j=1,64
	if(iframe(i,j).eq.-16389)isat=isat+1
	enddo
	enddo
	if(isat.ne.0)then
		isatcnt=isatcnt+1
		go to 98
	endif
c	
	ped=0.
	dark=0.
	idarkcnt=0
	ibaddark=0
	darksum=0.
	ipedcnt=0
	ibadped=0
	pedsum=0.
	pednaive=0.
c	
c**********************************!	Pedestal calculation
c
	do j=1,64
	if(iframe(4,j).ne.0)then
		iarg=int4(iframe(4,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)/0.75
		pednaive=pednaive+arg
		ipedcnt=ipedcnt+1
		pedsv(ipedcnt)=arg
		endif
	if(iframe(318,j).ne.0.and.j.ne.64)then
		iarg=int4(iframe(318,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)/0.75
		pednaive=pednaive+arg
		ipedcnt=ipedcnt+1
		pedsv(ipedcnt)=arg
	endif
	enddo
	pednaive=pednaive/127.
	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(127,pedsv,indxx)
	pedtest=pedsv(indxx(63))		!here compare with the naive pedestal median
c	
c	Here exclude "white measles", bad bands, and high pedestals above 144 electrons...
c	
	do j=1,127
	    if(pedsv(indxx(j))-pedtest.gt.2.)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.115)ibadped=1		!112=scaled from eng. data
	if(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.1)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,64
	if(iframe(5,j).ne.0)then
		iarg=int4(iframe(5,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)/0.75-ped
		darksum=darksum+arg
	endif
	if(iframe(317,j).ne.0.and.j.ne.64)then
		iarg=int4(iframe(317,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)/0.75-ped
		darksum=darksum+arg
	endif
	enddo
	darknaive=darksum/127.
	darkscale=darknaive/patdark	!naive dark current yields scale factor to "patdark"
c	note for engineering data this was raised to the 1.5 power: seems not needed here	
	if(darkscale.gt.1.0)darkscale=1.0	!this assumes "high" pattern is used
c
c	Now go back through using scale factor to filter CR hits from dark calculation
c 	
	darksum=0.
	do j=1,64
	if(iframe(5,j).ne.0)then
		iarg=int4(iframe(5,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)/0.75-ped	
		idarkcnt=idarkcnt+1
		darksum=darksum+arg
		jarg=1+2*(j-1)
		if(k.eq.1)darksv(jarg)=pattern(5,j)
		darkmed(jarg)=arg
		ddark=arg-pattern(5,j)*darkscale	!+patmore(i,j))*darkscale
		if(ddark.gt.patcr)then	!500 e-: see also just below...
		    idarkcnt=idarkcnt-1
		    darksum=darksum-arg
		    darkmed(jarg)=10000.
		    ddark=0.
		endif
	endif
	enddo
	do j=1,63
	if(iframe(317,j).ne.0)then
		iarg=int4(iframe(317,j))
		if(iarg.lt.0)iarg=iarg+65536
		arg=float(iarg)/0.75-ped
		idarkcnt=idarkcnt+1
		darksum=darksum+arg
    	        jarg=2+2*(j-1)
		if(k.eq.1)darksv(jarg)=pattern(317,j)
		darkmed(jarg)=arg
		ddark=arg-pattern(317,j)*darkscale	!+patmore(i,j))*darkscale
		if(ddark.gt.patcr)then
		    idarkcnt=idarkcnt-1
		    darksum=darksum-arg
		    darkmed(jarg)=10000.
		    ddark=0.
		endif
	endif
	enddo
c
c	Get median:"interpolate within the histogram..."
c
	call indexxx(127,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,127-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
c
c************************************************************
	if(ibaddark.eq.0.or.k.eq.1)then
		pedlast=ped
		pedrecover=0
		dark=median
	else	
		trash=trash+1
c			print 1001,k,pednaive-pedinit,ped-pedinit,ipedcnt,dark,median,idarkcnt,infile
1001			format('BAD DARK ',i6,', ped='2f6.2,', ipedcnt=',i3,' dark = ',2f5.2,', darkcnt=',i3,1x,a43)
		ped=-9999.
		dark=-9999.
		go to 98
	endif
c	
c	if(k.ge.1)go to 97	!$$$$$	bail quick	
c
	sumall=0.
	isumall=0
	patcorr=dark/patdark
	do j=1,64
	iless=0
	do i=1,318
	    frame(i,j)=66000.
	    if(iframe(i,j).ne.0)then
		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(i.gt.4.and.i.lt.318)arg=arg-pattern(i,j)*patcorr
		frame(i,j)=arg
		sumall=sumall+arg
		isumall=isumall+1
		if(frame(i,j).lt.framemin(i,j))framemin(i,j)=frame(i,j)
		if(iorbit.gt.1.and.frame(i,j)-frame1(i,j).lt.-10.)iless=iless+1
	    endif
	enddo
	if(iless.ge.10)print *,infile,': ***WARNING, ',iless,' lower pixels (< 10) in row ',j
	enddo
c	
97	continue
98	continue
	enddo
	print *,'reached end, 6 on the list were hand killed'
	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
	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(127),N,I,J,L,IR,INDXT	!Absoft Compiler dislikes INDX(N), etc.
      real*4 ARRIN(127),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