c
c name:		CRX_ped_dark_cam3.for
c purpose:	determines pedestal & dark current, finds black/white measles & cosmic rays
c inputs:	*.nic data files listed in "namelistCRX.txt"
c outputs:	before and after *.grd files 
c    "		Note here that the pedestal, dark, and CR information incorporated into the "trailer"
c    "		The output file xxx.dat presents information for graphs of the various quantities
c calls:	subroutines  measles (+ plane), and INDEXX 
c mods:		April 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 CRX_ped_dark_cam3
	implicit none
	integer*1 trail1(75),trail2a(94),trail2b(177),tdum1(12),trail3(21),tdum2(28),trail4(69)
	integer*2 iframe1(636,128),iframe2(636,128),head(3)
	integer*4 iframe(636,128),klast,itriangle,makemap
	integer*4 h,m,s,hh,mm,ss,i,j,k,kmax,kmin,iarg,iarg1,iarg2,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),k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,indxx(511)	!this last one a +1: strange
	integer*4 mode,id,jd,jarg,iqbad,iqcnt,iorbit,iarglast,med1,med2
	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 minresponse(1524),mingrit(1524),pixsum,pixdif
	real*4 pednaive,darksv(510),patdark,darkdist(510,6000),darkmed(510),darknaive
	real*4 pedsv(510),darkscale,median,dmed,pedtest
	real*4 camtemp1,camtemp2,camtempsave,weight,weightsum,tempmean
	real*4 temp1,glaremap(636,128),glarecnt(636,128)
	real*8 q(4),qq(4),t,t0,torbit,targ,tprev
	character*48 out3file,out2file,out1file
c	character infile*37,infile1*37,infile2*37	!old file naming convention
	character infile*47,infile1*47,infile2*47	!new file naming convention, shifted over 10
	character doy(15)*3,ddoy*3  
c	character nname*15,outname*25			!old file naming convention
	character nname*14,outname*24			!new file naming convention
	character hhead(5)*10,dum*1
	character ddate(15)*8,dddate*8
	data dum /' '/
c
c	2003
c	
c	data infile   /'G:/20030DOY/c3frm_2003_DOY_hhmmss.raw'/
c	data infile1   /'G:/20030DOY/c3frm_2003_DOY_hhmmss.raw'/
c	data infile2   /'G:/20030DOY/c3frm_2003_DOY_hhmmss.raw'/
c	data outname  /'c3frm_2003_DOY_hhmmss.nic'/		
c	data out1file /'G:/SEAN/buffer3/Crx/c3frm_2003_DOY_hhmmsscr1.grd'/	!cosmic-ray not removed file...
c	data out2file /'G:/SEAN/buffer3/Crx/c3frm_2003_DOY_hhmmsscrx.grd'/	!cosmic-ray removed file...
c	data out3file /'G:/SEAN/buffer3/Crx/c3avg_2003_DOY_hhmmsscrx.grd'/	!cosmic-ray removed binary file...
c
c	2004
c
c	old file naming convention
c
c	data infile   /'G:/20040DOY/c3frm_2004_DOY_hhmmss.raw'/
c	data infile1   /'G:/20040DOY/c3frm_2004_DOY_hhmmss.raw'/
c	data infile2   /'G:/20040DOY/c3frm_2004_DOY_hhmmss.raw'/
c	data outname  /'c3frm_2004_DOY_hhmmss.nic'/
c
c	new file naming convention
c
	data infile    /'D:/Sean/Buffer3/Nic288/c3m1_2004_DOY_hhmmss.raw'/
	data infile1   /'D:/Sean/Buffer3/Nic288/c3m1_2004_DOY_hhmmss.raw'/
	data infile2   /'D:/Sean/Buffer3/Nic288/c3m1_2004_DOY_hhmmss.raw'/
	data outname   /'c3m1_2004_DOY_hhmmss.nic'/	
c
	data out1file /'D:/SEAN/buffer3/Crx/c3frm_2004_DOY_hhmmsscr1.grd'/	!cosmic-ray not removed file...
	data out2file /'D:/SEAN/buffer3/Crx/c3frm_2004_DOY_hhmmsscrx.grd'/	!cosmic-ray removed file...
	data out3file /'D:/SEAN/buffer3/Crx/c3avg_2004_DOY_hhmmsscrx.grd'/	!cosmic-ray removed binary file...
	
	data ihist /250*0/
c	data ddate /'11/23/03','11/24/03','05/26/03','05/27/03','05/28/03','05/29/03','05/30/03','05/31/03','06/01/03',
c	data ddate /'05/31/03','11/24/03','05/26/03','05/27/03','05/28/03','05/29/03','05/30/03','05/31/03','06/01/03',
c	data ddate /'04/19/04','11/24/03','05/26/03','05/27/03','05/28/03','05/29/03','05/30/03','05/31/03','06/01/03',
	data ddate /'10/13/04','10/14/04','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'/
c	data doy /'327','328','146','147','148','149','150','151','152','153','154','155','156','157','158'/
c	data doy /'151','328','146','147','148','149','150','151','152','153','154','155','156','157','158'/
c	data doy /'110','328','146','147','148','149','150','151','152','153','154','155','156','157','158'/
	data doy /'287','288','146','147','148','149','150','151','152','153','154','155','156','157','158'/
	data badsum /0./
	data pattern /81408*0./
	data ipcnt /81408*0/
	data n /81408/
	data trash,iffix,iqcnt,isunmoon,isatcnt /5*0/
	data patcr /55.52/	!55.52 adu's = 1000 electrons for 2 x 2 binning
	data minresponse /1524*100000./
	data mingrit /1524*100000./
	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 darkdist /3060000*0./
	data weightsum /0./
	data iframe1,iframe2 /162816*0/
	data glaremap,glarecnt /162816*0.0/
c
	icrx=0			!icrx=1 means make before and after files
	iavg=0			!ibin=1 means make binary CRX files
	igraph=0		!igraph=1 means make GRAPHER output file and *.txt file for indexing...
	mode = 1
	makemap=1	
	id = 1272
	jd = 256
	id=id/mode
	jd=jd/mode
	patdark=119.447		!DOY287, 2004
	open (13,file='c3patt_2004_287_045015_2x2.grd',readonly)
	print *,'reading in c3patt_2004_287_045015_2x2.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
c	read(13,1500)patdark	!bummer, have to type in "patdark" above each time the pattern changes, mired in *.fts file...
1500	format(f10.4)
	print *,' patdark = ',patdark
	close(13)
c	open (10,file='namelist3_327.txt',readonly)
c	if(igraph.eq.1)open (12,file='cam3_327to328_CRX.dat')
c	open (10,file='namelist3_151.txt',readonly)
c	if(igraph.eq.1)open (12,file='cam3_151to151_CRX.dat')
c	open (10,file='namelist3_110.txt',readonly)
c	if(igraph.eq.1)open (12,file='cam3_110to110_CRX.dat')
	open (10,file='namelist3_288.txt',readonly)
	if(igraph.eq.1)open (12,file='cam3_287to288_CRX.dat')
	if(igraph.eq.1)open (21,file='cam3_287to288_CRX.txt')
	pedinit=0.	!0.	!1066.
	pedinit2=0.
	pedrecover=0
	icrcnt=0
	iavgcnt=0
c	open(20,file='DOY288_del_glare3.dat')
	
	t0=+2014.5D0	!DOY144
c********************************************************************
c	kmax =  2832	!DOY327
c	kmax =  2962	!DOY151
c	kmax =  1549	!DOY110, a bit left at the end of this, but spotty coverage...
	kmax =  6267	!DOY288, the whole file
c	k1=10289
	k1=871		!DOY288
	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	!list includes 19 frames after major calibration in DOY 158
	kmin=1
c	kmax=200
c	kmin=202	!DOY288, 2004: earlier - 1st eclipse
c	kmax=366	!DOY288, 2004: earlier - 1st eclipse
c	kmin=1550	!DOY288, 2004: earlier - 2nd eclipse
c	kmax=1745	!DOY288, 2004: earlier - 2nd eclipse
c	kmin=2886	!DOY288, 2004: earlier - 3rd eclipse
c	kmax=3088	!DOY288, 2004: earlier - 3rd eclipse
c	kmin=1726	!DOY288, 2004: later - 1st eclipse
c	kmax=1860	!DOY288, 2004: later - 1st eclipse
c	kmin=3074	!DOY288, 2004: later - 2nd eclipse
c	kmax=3269	!DOY288, 2004: later - 2nd eclipse
	kmin=4410	!DOY288, 2004: later - 3rd eclipse
	kmax=4612	!DOY288, 2004: later - 3rd eclipse
c	kmax=4740	!DOY288, 2004: whole file differenced
	do k=1,kmax
	if(k.eq.1)pedlast=pedinit
c
	read (10,10)nname,h,m,s
c10	format(39x,a15,3i2.2)	!old file naming convention
10	format(39x,a14,3i2.2)	!new file naming convention
	if(k.lt.kmin)go to 98
c	if(k.eq.4.or.k.eq.243.or.k.eq.330.or.k.eq.1331.or.k.eq.1334.or.k.eq.1336.or.k.eq.1337)go to 98	     !DOY327
c	if(k.eq.1339.or.k.eq.1340.or.k.eq.1343.or.k.eq.1346.or.k.eq.1347.or.k.eq.1349.or.k.eq.1350)go to 98	!DOY327
c	if(k.eq.1352.or.k.eq.1354.or.k.eq.1356.or.k.eq.1358.or.k.eq.1363.or.k.eq.1371.or.k.eq.1382)go to 98	!DOY327
c	if(iabs(k-1764).le.13)go to 98	!DOY327
c	if(k.eq.248.or.k.eq.447.or.k.eq.451.or.k.eq.455.or.k.eq.458.or.k.eq.1123.or.k.eq.1127)go to 98	!DOY151
c	if(k.eq.1129.or.k.eq.1133.or.k.eq.1142)go to 98	!DOY151
c	the following for DOY110, 2004
c	if(k.eq.6.or.k.eq.27.or.k.eq.48.or.k.eq.53.or.k.eq.60.or.k.eq.944.or.k.eq.1507.or.k.eq.1543.or.k.eq.1546)go to 98
c	if(k.lt.895.or.k.gt.1075)go to 98	!DOY110, 2004
	if(k.eq.1542.or.k.eq.3144.or.k.eq.3213.or.k.eq.3503.or.k.eq.3622)go to 98	!DOY288, 2004
	klast=k
	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
	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
c	write(infile1(9:11),'(a3)')ddoy		!old naming convention
c	write(infile1(13:27),'(a15)')nname	!old naming convention
c	write(infile1(28:33),'(3i2.2)')h,m,s	!old naming convention
c	write(outname(1:15),'(a15)')nname	!old naming convention
c	write(outname(16:21),'(3i2.2)')h,m,s	!old naming convention
c
c	write(infile1(20:22),'(a3)')ddoy	!new naming convention	!skip this step for DOY288,2004
	write(infile1(24:37),'(a14)')nname	!new naming convention
	write(infile1(38:43),'(3i2.2)')h,m,s	!old naming convention
	write(outname(1:13),'(a14)')nname	!new naming convention
	write(outname(15:20),'(3i2.2)')h,m,s	!new naming convention
c	
	ddoy=doy(2)	!special for DOY288, 2004
	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
	ddoy=doy(2)
	endif
	write(nname(11:13),'(a3)')ddoy
c	
c	write(infile2(9:11),'(a3)')ddoy		!old naming convention
c	write(infile2(13:27),'(a15)')nname	!old naming convention
c	write(infile2(28:33),'(3i2.2)')hh,mm,ss	!old naming convention
c	
c	write(infile2(20:22),'(a3)')ddoy	!new naming convention	!skip this step for DOY288,2004
	write(infile2(24:37),'(a14)')nname	!new naming convention
	write(infile2(38:43),'(3i2.2)')hh,mm,ss	!new naming convention
c	print *,k,dum,infile1,dum,infile2
	open(11,file=infile1,form='binary',readonly)
	read(11)head,((iframe1(i,j),i=1,636),j=1,128),trail1,q,trail2a,camtemp1,trail2b,tdum1,trail3,tdum2,trail4
	close(11)
	open(11,file=infile2,form='binary',readonly)
	read(11)head,((iframe2(i,j),i=1,636),j=1,128),trail1,q,trail2a,camtemp2,trail2b,tdum1,trail3,tdum2,trail4
	close(11)
	camtempsave=camtemp1
	tprev=t
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 2x2 binned pixels
c
	isat=0
	do i=11,632
	do j=1,128
	if(iframe1(i,j).lt.0.and.iframe1(i,j).gt.-10)isat=isat+1
	if(iframe2(i,j).lt.0.and.iframe1(i,j).gt.-10)isat=isat+1
	enddo
	enddo
	if(isat.ne.0)then
		isatcnt=isatcnt+1
		print *,' frame ',k,': saturated pixels...',isat
		go to 98
	endif
c	iqbad=0
c	do i=1,4
c		if(k.ge.11.and.dabs(q(i)-qq(i)).gt.1.0D-06)then
c			iqbad=iqbad+1
c		endif
c	enddo
c	if(iqbad.gt.0)then
c		do i=1,4
c			q(i)=qq(i)
c		enddo
c		print *,k,': Replaced bad quaternions on ',nname,h,m,s
c		iqcnt=iqcnt+1
c	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(i,j))
		iarg2=int4(iframe2(i,j))
		if(iarg1.lt.0)iarg1=iarg1+66536
		if(iarg2.lt.0)iarg2=iarg2+65536
c		iframe(i,j)=iarg1-iarg2		!frame 1 - frame 2
		iframe(i,j)=iarg2-iarg1		!frame 2 - frame 1
		if(iarg1.gt.2500.or.iarg2.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
		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	open (30,file='temp1.dat')
c	write(30,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...
c	    if(pedsv(indxx(j))-pedtest.gt.-4.)then
	    if(pedsv(indxx(j))-pedtest.gt.10.)go to 20	!exclude high pedestals above 180 electrons...
	    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		!DOY288 (320 works for individual frames)
	ped=pedlast				!pedlast is a final resort value when everyone fails
	if(ipedcnt.gt.0)ped=pedsum/float(ipedcnt)
c	print *,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" (was 0.5 before, see also below)
	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	!200 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(k-kmin+1.le.6000.and.k-kmin+1.gt.0)then
	do j=1,510
		if(darkmed(j).ne.0..and.darkmed(j).ne.10000.)darkdist(j,k-kmin+1)=darkmed(j)
	enddo
	endif
	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)
c	print 2000,(i,darkmed(indxx(i)),i=1,510)
c	if(k-kmin+1.le.6000)then	
c	do j=1,510
c		darkdist(j,k-kmin+1)=darkmed(indxx(j))-ped
c	enddo
c	endif
	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 (was 0.76, see also above)
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
	if(k-kmin+1.le.6000.and.k-kmin+1.gt.0)then
	do j=1,510
		darkdist(j,k-kmin+1)=darkdist(j,k-kmin+1)-darksv(j)*dark/patdark
	enddo
	endif	
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.
		if(k-kmin+1.le.6000.and.k-kmin+1.gt.0)then
		do j=1,510
			darkdist(j,k-kmin+1)=0.
		enddo
		endif
		go to 98
	endif
c	weight=1.
c	weight=1.-float(iabs(k-299))/67.	!DOY288, 2004: earlier - 1st eclipse
	weight=1.-float(iabs(k-1647))/98.	!DOY288, 2004: earlier - 2nd eclipse
c	weight=1.-float(iabs(k-2987))/99.	!DOY288, 2004: earlier - 3rd eclipse
c	weight=1.-float(iabs(k-1793))/67.	!DOY288, 2004: later - 1st eclipse
c	weight=1.-float(iabs(k-3171))/98.	!DOY288, 2004: later - 2nd eclipse
c	weight=1.-float(iabs(k-4511))/99.	!DOY288, 2004: later - 3rd eclipse
c	
c	if(k.ge.1)go to 97	!$$$$$	bail quick	
c
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
	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	
	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
c			if(ikill(i,j).eq.1)frame(i,j)=-100.	!KILL measles pixels (otherwise neighbors average is put in)
			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
c	if(meascnt.gt.100)go to 98
	pixsum=pixsum/float(ipixsum)
	pixdif=pixdif/float(ipixdif)
c	write(20,56)k,glare,ped-pedinit,darkscale,dark,pixsum,meascnt/10,ipedcnt,idarkcnt
56	format(i6,5f10.2,3i6)
	do i=1,636
	do j=1,128
		if(iframe1(i,j).ne.0)then
		if(glare.gt.0.)then
c			if(frame(i,j).gt.-2..and.frame(i,j).lt.75.)then			!75 for DOY151, 15 for DOY327
			if(frame(i,j).gt.-10..and.frame(i,j).lt.100.)then		!75 for DOY151, 15 for DOY327
				glaremap(i,j)=glaremap(i,j)+frame(i,j)*glare
				glarecnt(i,j)=glarecnt(i,j)+glare
			endif
		endif
		endif
	enddo
	enddo
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.D0*targ)
	if(iarg.gt.1524)iarg=1524
	if(iarg.lt.iarglast)print *,'Beginning orbit #',iorbit,', at k = ',k,', ',infile1
	iarglast=iarg
c	pixsum=pixsum-minresponse(iarg)					!for use when minresponse table read in
c	pixdif=pixdif-mingrit(iarg)					!for use when mingrit table read in
c	if(minresponse(iarg).gt.pixsum) minresponse(iarg)=pixsum	!for use in forming minresponse table
c	if(mingrit(iarg).gt.pixdif) mingrit(iarg)=pixdif		!for use in forming mingrit table
c	pixsum=pixsum-0.025*float(meascnt)				!for use when minresponse table read in, BUT Cam 3?
c	
c	cut=pixdif+float(meascnt)/20.					!cut at line between pixdif=10 and meascnt=200
c	if(cut.gt.10.)then
c		print 6,infile,meascnt,ineigh,pixdif
6		format(' file = ',a46,': Bad measles + grit, ',2i4,f10.1)
c		trash=trash+1
c		go to 98
c	endif
c	do j=1,128
c		do i=11,632
c			if(frame(i,j).gt.0.)frame(i,j)=frame(i,j)-0.025*float(meascnt)
c			if(flatcorr(i,j).ne.0.)frame(i,j)=frame(i,j)/flatcorr(i,j)  !Flat Field Correction
c		enddo
c	enddo
c***		Old aurora condition follows	
c		if(pixsum.gt.3.0)then
c		write(19,19)dddate,h,hh,m,mm,s,ss,meascnt,pixsum,pixdif
c		print 19,dddate,h,hh,m,mm,s,ss,meascnt,pixsum,pixdif,k
c19		format('"',a8,1x,2i1,':',2i1,':',2i1,'"',2x,i6,2f10.2,i6)
c	endif		
97	if(igraph.eq.1.and.k.gt.2)write(12,120)dddate,h,m,s,ped-pedinit,dark,meascnt,pixsum,pixdif,glare,camtemp1
	if(igraph.eq.1.and.k.gt.2)write(21,33)outname,ddoy,ped,dark,q,meascnt,pixsum,pixdif,glare
c	if(mod(k,100).eq.0)then
c	print 12,dddate,h,m,s,ped-pedinit,dark,meascnt,pixsum,pixdif,glare,k
12	format('"',a8,1x,i2.2,':',i2.2,':',i2.2,'"',2x,2f10.3,i6,3f10.2,i7)
120	format('"',a8,1x,i2.2,':',i2.2,':',i2.2,'"',2x,2f10.3,i6,4f10.2)
33	format(a25,a3,2f10.3,4f13.9,i6,3f10.2)
c	endif
c	write(out1file(32:41),'(a10)')infile1(24:33)
c	write(out2file(32:41),'(a10)')infile1(24:33)
c	write(out3file(32:41),'(a10)')infile1(24:33)
c	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)
c	iavgcnt=0
c	do j=1,128
c	do i=1,636
c	frameavg(i,j)=0.
c	enddo
c	enddo
c	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	
c	print *,'Start writing output files'
c	open (13,file='pattern3.grd')
c	print *,' writing out pattern3.grd'
c	write(13,13)
c13	format('DSAA'/'636 128'/'0 636'/'0 128'/'-10 100')	
c	do j=1,128
c	write(13,15)(pattern(i,j),i=1,636)
c15	format(20f9.3)
c	enddo
c	close(13)
c	
c	open (17,file='min3_2.dat')
c	do i=1,1524
c	write(17,17)i,minresponse(i)
c	enddo
c	close(17)	
c	open (17,file='mingrit3_2.dat')
c	do i=1,1524
c	write(17,17)i,mingrit(i)
c	enddo
c	close(17)
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	print *,'quaternions replaced on ',iqcnt
c	do i=121,250
c	print *,i-120,ihist(i)
c	enddo
	if(igraph.eq.1)then
	print *,'It aint done until the fat lady sings...'
	open(10,file='darkdist3ecl.grd')
	write(10,13)
13	format('DSAA'/'6000 510'/'0 6000'/'0 2044'/'-5 50')
	do j=1,510
	do i=1,6000
	if(darkdist(j,i).gt.999.99)darkdist(j,i)=999.99
	if(darkdist(j,i).lt.-99.99)darkdist(j,i)=-99.99
	enddo
	write(10,130)(darkdist(j,i),i=1,6000)	
130	format(20f7.2)	
	enddo
	close(10)
	endif
c
c	form, write output glare map
c
	if(makemap.eq.1)then
	open(10,file='d:\sean\eclipse\glare3_288_3rd_l.grd')		!glaremapall.grd')	!eclipse\glare3_288_2nd_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.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 INDX(81408),N,I,J,L,IR,INDXT	!Absoft Compiler dislikes INDX(N), etc.
      real*4 ARRIN(81408),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