! Started modifying this on 2003/06/26
! Original version should be on backup CD prior to this.

	program ped_dark2

	implicit none

	include		'openfile.h'

	integer		 nx
	parameter	(nx=1272)
	integer		 ny
	parameter	(ny= 256)
	integer		 nxy
	parameter	(nxy=nx*ny)

	character	cInfo*512
	integer		iframe(nx,ny)
	integer		h,hh,m,mm,s,ss,i,j,k,kmax,kmin,iarg,iitime,kdel,icolum,i1,i2,ioutcnt,ineigh
	integer		idarksum,idarkcnt,ipedcnt,ipedsum,ibadped,ibaddark,icrcnt,meascnt
	integer		framein(nx,ny),icrneig,inegmeas
	integer		indx(nxy),kskip
	real		ttime,t0,tlast,dark,darklast,ped,pedfirst,pedlast,arg,arg1
	real		frame(nx,ny)
	character	cFile*256	/'G:/20030227/c2frm_2003_058_hhmmss.nic'/		!data starting 2-27-03
	character	nname*15
	character	hhead(5)*10
	character	ddate(3)*8	/'02/27/03','02/28/03','03/01/03'/
	character	dddate*8
	integer		ihist(250)	/250*0/
	integer		i1		/228/
	integer		i2		/301/
	integer		ioutcnt		/  0/
	real		badsum		/0.0/
	real		pattern(nx,ny)	/nxy*0./
	real		fleast(nx,ny) 	/nxy*100000./
	integer		ipcnt(nx,ny)	/nxy*0/
	integer		n		/nxy/

	logical		bOK, bReadNic, bWriteNic

	open (12,file='least2.grd')
	read (12,'(A)')hhead
	do j=1,ny
	    read(12,'(20F9.1)')(fleast(i,j),i=1,nx)
	enddo
	close(12)

	open (12,file='temp.dat')

	icolum=1001
	kskip=16175
	kmin=2		!skip first frame, different exposure
	kmax=23484	!camera 2 all, into 3-1	!16174 add 2-28	!3658 for just 2-27
	!kmax=10000		
	kdel=kmax-kmin+1
	iitime=nint(6.83*float(kdel))

	open (10,file='namelist2.txt')

	read (10,'(39X,A15,6I1)')nname,h,hh,m,mm,s,ss
	t0 = 3600.*float(10*h+hh)+60.*float(10*m+mm)+float(10*s+ss)

	do k=2,kmin-1
	    read (10,'(39X,A15,6I1)') nname,h,hh,m,mm,s,ss
	end do

	do k=kmin,kmax

	    read (10,'(39X,A15,6I1)') nname,h,hh,m,mm,s,ss
	    ttime=3600.*float(10*h+hh)+60.*float(10*m+mm)+float(10*s+ss)

	    if (k-kmin .gt. 9 .and. k .lt. kskip) go to 98	!***start at "kskip"	

	    ttime = ttime-t0
	    if (k.ge. 3659) ttime=ttime+86400.
	    if (k.ge.16175) ttime=ttime+86400.

	    if (ttime-tlast .lt. 0.0) write (3,'(A)') 'file '//cfile(:itrim(cfile))//' out of order time'

	    write(cFile(13:27),'(A15)') nname
	    write(cFile(28:33),'(6I1)') h,hh,m,mm,s,ss
	    if(k.eq.3659.or.(k-kmin.gt.9.and.kskip.gt.3659.and.kskip.lt.16175))write(cFile(9:11),'(i3)')i1
	    !if(k.eq.16175.or.(k-kmin.gt.9.and.kskip.gt.16175))print *,cFile(:itrim(cFile))
	    if(k.eq.16175.or.(k-kmin.gt.9.and.kskip.gt.16175))write (cFile(9:11),'(i3)') i2

	    bOK = bReadNic(OPN__REOPEN+OPN__NOMESSAGE+OPN__STOP, cFile, nxy, nxx, nyy, iframe, cInfo)

	    ipedcnt	= 0
	    ipedsum	= 0
	    idarkcnt	= 0
	    idarksum	= 0
	    ibadped	= 0
	    ibaddark	= 0
	    icrcnt	= 0
	    icrneig	= 0
	    meascnt	= 0
	    ineigh	= 0
	    inegmeas	= 0

	    do j=1,ny

		do i=13,20
		    iarg = iframe(i,j)

		    if (iarg .ne. 0) then
			if(i.le.16)then
			    if(k-kmin .le. 9 .or. abs(float(iarg)-pedlast).lt.10.)then
				ipedsum=ipedsum+iarg
				ipedcnt=ipedcnt+1
			    endif
			endif
			if(i.ge.17)then
			    if(k-kmin .le. 9 .or. abs(float(iarg)-pedlast-darklast).lt.10.)then
				idarksum=idarksum+iarg
				idarkcnt=idarkcnt+1
			    endif
			endif
		    endif

		enddo

		do i=1265,nx
		    iarg = iframe(i,j)

		    if (iarg .ne. 0) then
			if(i.ge.1269)then
			    if(k-kmin .le. 9 .or.abs(float(iarg)-pedlast).lt.10.)then
				ipedsum=ipedsum+iarg
				ipedcnt=ipedcnt+1
			    endif
			endif
			if(i.le.1268)then
			    if(k-kmin .le. 9 .or. abs(float(iarg)-pedlast-darklast).lt.10.)then
				idarksum=idarksum+iarg
				idarkcnt=idarkcnt+1
			    endif
			endif
		    endif

		enddo

	    enddo

	    if(ipedcnt .lt.1850) ibadped =1
	    if(idarkcnt.lt.1000) ibaddark=1

	    ped =float(ipedsum)/float(ipedcnt)
	    dark=float(idarksum)/float(idarkcnt)-ped
	    dark=dark-0.0018*float(1965-idarkcnt)	!ad hoc dark current adjustment

	    !print *,ped,dark

	    if (k .eq. kmin) then
	    	write (*,'(A,/,A,F9.2,A,I5)') 'file '//infile(:itrim(infile)),'First pedestal=',ped,',plotting column #',icolum
		pedlast  = ped
		pedfirst = ped
		tlast    = ttime

		do j=1,ny
		    do i=17,1268
			if(iframe(i,j) .ne. 0) then
			    pattern(i,j) = float(iframe(i,j))-dark-ped
			    ipcnt  (i,j) = 1
			!else
			!    pattern(i,j) = 0.0
			!    ipcnt  (i,j) = 0
			end if
		    end do
		end do

	    else if (k-kmin .lt. 6 .or. k-kmin .eq. 9) then

		do j=1,ny
		    do i=17,1268

			if(iframe(i,j).ne.0)then
			    arg = float(iframe(i,j))-dark-ped

			    if (i .eq. 701 .and. j .eq. 47) then	!701,47 is real bad...

				pattern(i,j) = pattern(i,j)+arg
				ipcnt  (i,j) = ipcnt  (i,j)+1

			    else if (k .eq. kmin+1) then

				arg1 = arg-pattern(i,j)
				if (arg1 .lt. -22.2) then		!22.2 ADUs is 100 electrons
				    !print *,'1st one bad',k-1,i,j,pattern(i,j),arg
				    pattern(i,j) = arg
				else if (arg1 .le. +22.2) then
				    pattern(i,j) = pattern(i,j)+arg
				    ipcnt  (i,j) = ipcnt  (i,j)+1
				!else
				!   print *,'2nd one bad',k-1,i,j,arg,pattern(i,j)
				end if

			    else
				arg1=arg-pattern(i,j)/float(ipcnt(i,j))
				if (arg1 .le. +22.2) then
				    pattern(i,j) = pattern(i,j)+arg
				    ipcnt  (i,j) = ipcnt  (i,j)+1
				!else
				!    print *,'>=3rd one bad',k-1,i,j,arg,pattern(i,j)/float(ipcnt(i,j))
				end if

			    end if
			end if

		    end do
		end do

	    end if

	    if (k .eq. kmin+9) then
		do j=1,ny
		    do i=17,1268
			if (ipcnt(i,j) .gt. 0) pattern(i,j) = pattern(i,j)/ipcnt(i,j)
		    end do
		end do
		write (*,*) 'Done forming initial pattern'
	    endif

	    if(ibadped .eq. 1 .or. ibaddark .eq. 1 .or. dark .le. 1.17) then
		ioutcnt=ioutcnt+1
		write (5,'(A,A,2I5,F10.3)') ' file '//infile(:itrim(infile)),' Bad count or low dark ',ipedcnt,idarkcnt,dark

	    else
		pedlast=ped
		darklast=dark

		if (k-kmin .gt. 11) then				!skip initial frames with shutter closed

		    do j=1,ny
			do i=1,ny
			    iarg = iframe(i,j)
			    arg = float(iarg)-ped
			    if (arg .lt. -10.0 .and. iarg .ne. 0 .and. i .ge. 21 .and. i .le. 1264) then	!stamp out "negative measles"
				write (*,'(2A,2I5,A,I6)') 'file '//infile(:itrim(infile),': pixel ',i,j,',',iarg,' is less than pedestal'
				iarg=0
				inegmeas=inegmeas+1
			    endif
			    if(i.ge.16.and.i.le.1268)arg=arg-dark-pattern(i,j)
			    frame(i,j)=arg
			    frame(i,j)=frame(i,j)-fleast(i,j)		!remove bland background using "fleast(i,j)"
			    if(iarg.eq.0)frame(i,j)=0.
			    framein(i,j)=nint(frame(i,j))
			    !if(frame(i,j).lt.fleast(i,j))fleast(i,j)=frame(i,j)	!for making fleast file
			enddo
		    enddo

		    call IndexR4(1,nxy,1,nxy,frame,indx)
		    call measles(nx,ny,frame,indx,meascnt,ineigh,ihist)
		    if(meascnt.ge.5000) write (6,'(2A,2I5)') 'file '//infile(:itrim(infile),': Bad measles, ',meascnt,ineigh

		    call cosmicray(nx,ny,frame,indx,icrcnt,icrneig)

		    do j=1,ny
			do i=1,nx
			    if (frame(i,j) .ne. 0.0 .and. frame(i,j) .ne. -10.0) frame(i,j)=frame(i,j)+fleast(i,j)	!replace bland background
			    iframe(i,j) = nint(frame(i,j))	!This location for binary out...
			enddo	
		    enddo

		    dddate = ddate(1)
		    if(k.ge.3659)dddate=ddate(2)
		    if(k.gt.16175)dddate=ddate(3)

		    if (k-kmin .ge. 10) then
			write(12,('"',a8,1x,2i1,':',2i1,':',2i1,'"',2x,2f10.3,8i6)) dddate,h,hh,m,mm,s,ss,ped-pedfirst,dark,ipedcnt,idarkcnt,inegmeas,meascnt,ineigh,icrcnt,icrneig
			if(mod(k,100).eq.0)then
			    write (*,('"',a8,1x,2i1,':',2i1,':',2i1,'"',2x,2f10.3,8i6)) dddate,h,hh,m,mm,s,ss,ped-pedfirst,dark,ipedcnt,idarkcnt,inegmeas,meascnt,ineigh,icrcnt,icrneig,k
			endif

			bOK = bWriteNic(cFile+'crx', '', nx, ny, 2, iframe, cInfo)
		    endif

		end if

	    end if

98	    continue

	enddo


	close(10)
	close(12)
c	print *,'Start writing output files'
c	open (13,file='temp.grd')
c	write(13,130)kmax
c130	format('DSAA'/i6,' 256'/'0 400'/'0 256'/'-100 100')
c	open (13,file='pattern2.grd')
c	write(13,13)
c13	format('DSAA'/'1272 256'/'0 1272'/'0 256'/'-10 100')	
c	open (13,file='least2.grd')
c	write(13,13)
c13	format('DSAA'/'1272 256'/'0 1272'/'0 256'/'0 70')
c	do j=1,256
c	   write(13,'(20F9.1)')(pattern(i,j),i=1,nx)
c	   write(13,'(20F9.1)')(fleast(k,j),k=1,nx)
c	enddo
c	close(13)
c	do j=1,250
c	print 1000,100*(j-1),ihist(j)
c1000	format(2i10)
c	enddo	

	print *,'reached end,',ioutcnt,' discarded out of ',kmax-kmin+1

	stop 

	end
c*******
	subroutine cosmicray(nx,ny,frame,indx,icrcnt,icrneig)		!This version for 1x1 binned data only
	!Note*** this version doesn't try subtracting a "bland background"... might improve performance.	

	integer		indx (*)
	real		frame(nx,ny)
	real		dcut	/111./				!500-electron cut
	logical		bail

	n = nx*ny

	do m=1,n
	    iarg=1+n-m
	    j=(indx(iarg)-1)/nx
	    i=indx(iarg)-nx*j
	    j=j+1

	    ! What about partially covered column 21

	    if (21 .lt. i .and. i .lt. 1264 .and. 1 .lt. j .and. j .lt. ny) then begin

		if (frame(i,j) .ne.0.0 .and. frame(i,j) .ne. -10.0 .and. frame(i,j) .lt. dcut) return	!quit if below threshold

		! To be a cosmic ray, the pixel exceed the average, then further none of the 24 closest neighbors
		! can be greater, and the pixel must stand more than 500 electrons (111 ADUs) above the
		! average of the 24 closest neighbors,

		! First check if any neighbor or next-neighbor is more, then adjust background (24-pixel average)

		bail = .FALSE.

		idifcnt=0
		dsum=0.

		do jarg=max(1,j-2),min(j+2,ny)
		    do iarg=max(22,i-2),min(i+2,1263)
			if (iarg .ne. i .or. jarg .ne. j) then
			    arg = frame(iarg,jarg)
			    if (arg .ne. 0.0 .and. arg .ne. -10.0) then
				bail = bail .or. frame(i,j) .lt. arg	!must be biggest
				idifcnt=idifcnt+1
				dsum=dsum+arg
			    end if
			end if
		    end do
		end do

		if (.not. bail) then

		    backg=dsum/float(idifcnt)
		    diff=frame(i,j)-backg

		    bail = diff .lt. dcut	!must stand more than 500 electrons above average

		end if

		if (.not. bail) then

		    ! The overage sum cannot increase more than 3.75X when adding the <8 immediate neighbors,
		    ! "stars do, streaks don't", 
		    ! and the pixel must stand above at least 7*sqrt(variance) of the 24 other pixels.

		    sumsq=0.
		    isumcnt=0
		    difsum=0.

		    do jarg=max(1,j-2),min(j+2,ny)
			do iarg=max(22,i-2),min(i+2,1263)
			    if (iarg .ne. i .or. jarg .ne. j) then
				arg=frame(iarg,jarg)
				if (arg .ne. 0.0 .and. arg .ne. -10.0)then
				    arg=arg-backg
				    isumcnt=isumcnt+1
				    sumsq=sumsq+arg*arg
				    if (iabs(iarg-i).le.1.and.iabs(jarg-j).le.1) difsum=difsum+arg
				 end if
			    end if
			end do
		    end do

		    if (isumcnt .gt. 2) then		!Always true?

			variance=sumsq/float(isumcnt-1)
			sdev=sqrt(variance)			!Here's the only transcendental (could be avoided)...

			! Here  a cosmic ray" must meet both of the remaining criteria
			! The 7-sigma criterion is not sensitive to background, but the 3.75x criterion is...

			if (diff/sdev.gt.7.0 .and. difsum/diff .lt. 3.75) then

			    ! So it IS a cosmic ray...

			    frame(i,j) = -10.0			!Stamp it out (could use backg...)
			    icrcnt=icrcnt+1

			    ! Now do neighbors within the 5 x 5 block

			    do jarg=max(1,j-2),min(j+2,ny)
				do iarg=max(21,i-2),min(i+2,1264)!note extend neighbors to include cols 21, 1264.
				    arg=frame(iarg,jarg)
				    if(arg.ne.0..and.arg.ne.-10.0)then
					arg=arg-backg				!take away 24-pix backg
					if(arg.gt.dcut.and.arg.gt.0.1*diff.and.arg.gt.2.*sdev)then
					    !print *,i,j,ii,jj,frame(iarg,jarg),backg,sdev			    		    
					    frame(iarg,jarg) = -10.0
					    icrneig=icrneig+1
					endif
				    endif
				enddo
			    enddo

			end if

		    end if

		end if

	    end if

	enddo

	!print *,' CR m = ',m	
	return
	end
c*******
	subroutine measles(nx,ny,frame,indx,meascnt,ineigh,ihist)

	! This subroutine is a bare-bones version of a cosmic-ray spot remover which uses only the
	! ADUs departure of a pixel from its neighbors-average z0 to trigger its replacement/elimination,
	! accompanied by up to 3 neighbors if they fail the same threshold. The response must be above
	! a fixed value "dcut" and be at least double the 8-pixel average response z0.
 
	integer		ihist(*)
	integer		indx (*)
	real		frame(nx,ny)
	integer		nfix(-1:1,-1:1)
	real		z   (-1:1,-1:1)
	real		dcut		/222./				!1000-electron cut
	logical		bail

	n = nx*ny

	do m=1,n

	    iarg=1+n-m
	    j=(indx(iarg)-1)/nx
	    i=indx(iarg)-nx*j
	    j=j+1

	    ! bail if near edges
	    ! What about partially covered column 21

	    if (21 .lt. i .and. i .lt. 1264 .and. 1 .lt. j .and. j .lt. ny) then begin

		! Sanity check ??

		if (frame(i,j) .ne. 0. .and. frame(i,j) .lt. dcut) return	!quit if below threshold

		do jj=-1,1			! Pick up 3x3 block centered on i,j
		    do ii=-1,1
			z(ii,jj) = frame(i+ii,j+jj)
		    end do
		end do

		! bail if 3x3 group has empty pixel

		bail =	z(-1,-1) .eq. 0.0 .or. z(0,-1) .eq. 0.0 .or. z(1,-1) .eq. 0.0 .or.
     &			z(-1, 0) .eq. 0.0 .or. z(0, 0) .eq. 0.0 .or. z(1, 0) .eq. 0.0 .or.
     &			z(-1, 1) .eq. 0.0 .or. z(0, 1) .eq. 0.0 .or. z(1, 1) .eq. 0.0

		if (.not. bail) then

		    ! First check if any neighbor or next-neighbor is more,
		    ! then adjust background (8-pixel average)

		    call plane(z,z0)

		    ! Bail if not above threshold, or if not double the neighbors average

		    bail = z(0,0) .lt. dcut .or. z(0,0) .lt. 2.0*z0

		end if

		if (.not. bail) then

		    meascnt = meascnt+1

		    neigh=0

		    do jj=-1,1
			do ii=-1,1
			    if (ii .eq. 0 .and. jj .eq. 0) then		! Center pixel
				nfix(ii,jj) = 1
			    else
				if (z(ii,jj) .gt. 2.0*z0 .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
				else
				    nfix(ii,jj)=0
				end if
			    end if
			end do
		    end do

		    if (neigh .le. 3) then
			do jj=-1,1
			    do ii=-1,1
				if (nfix(ii,jj) .eq. 1) frame(i+ii,j+jj) = 0.0	!Put in zero...
			    end do
			end do

			if (neigh .gt. 0) ineigh = ineigh+1
		    else			!too many neighbors, just zero original pixel
			frame(i,j)=0.0		!put in zero...
		    end if

		end if

	    end if

	end do

	!print *,'measles m = ',m

	return
	end
c*******
	subroutine plane(z,z0)

	!Solves for a least-squares planar fit through the 8 z-values
	!surrounding [2,2], then reduces all 9 z's relative to the plane,
	!and returns with the reduced values.

	real	z(-1:1,-1:1)

	sumz=z(-1,0)+z(1,0)
	sumx=0.
	sumy=0.
	do i=-1,1
	    sumz=sumz+z( i,-1)+z(i,1)
	    sumx=sumx+z(-1, i)-z(1,i)
	    sumy=sumy+z( i,-1)-z(i,1)
	end do
	a=sumx/6.
	b=sumy/6.
	z0=sumz/8.
	do i=-1,1
	    do j=-1,1
		z(i,j)=z(i,j)-z0+a*float(i)+b*float(j)
	    enddo
	enddo

	return
	end