c
c name:		CRX_ped_dark_cam2.for
c purpose:	determines pedestal & dark current, finds measles
c inputs:	*.raw data files listed in "namelist2.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 files xxx.dat,txt present information for graphs/subseq of the various quantities
c calls:	subroutines  measles (+ plane), and INDEXX (2 versions) 
c mods:		April 2004, Andrew Buffington (UCSD; (858)-534-6630; abuffington@ucsd.edu)
c-
	program CRX_ped_dark_cam2
	implicit none
	integer*1 trail1(75),trail2a(94),trail2b(177),tdum1(12),trail3(21),tdum2(28),trail4(69)
	integer*2 head(3),iframe1(318,64),iframe2(318,64)
	integer*4 iframe(318,64),klast,jedge(180)
	integer*4 h,m,s,hh,mm,ss,i,j,k,kmax,kmin,iarg,ineigh,pedrecover,trash,ipcnt(318,64)
	integer*4 idarkcnt,ipedcnt,ibadped,ibaddark,ihist(250),icrcnt,meascnt,ikill(318,64),imed,med1,med2
	integer*4 icrneig,n,inegmeas,ipixsum,ipixdif,icrx,iavg,igraph,iffix,isunmoon,isat,isatcnt
	integer*4 indx(20352),k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,indxx(127),i10skip
	integer*4 mode,id,jd,jarg,iqbad,iqcnt,iorbit,iarglast,iavgcnt,itemp,iedgesave(180),indx180(180)
	real*4 dark,darksum,ped,pedlast,pedinit,pedinit2,pedsum
	real*4 frame(318,64),frameavg(318,64),frameout(318,64),badsum,pattern(318,64),arg,patcr,patcorr
	real*4 fframe(318,64),minresponse(1524),mingrit(1524),pixsum,pixdif,flatcorr(318,64)
	real*4 pednaive,darksv(127),patdark,darkdist(127,6000),darkmed(127),darknaive
	real*4 pedsv(127),darkscale,median,dmed,pedtest,ddark
	real*4 camtemp1,camtemp2,camtempsave,weight,weightsum,tempmean
	real*4 temp1,glaremap(318,64),glarecnt(318,64),edge(180),edgesave(180),edgedist(5265,180)
	real*8 q(4),qq(4),t,t0,torbit,targ,tprev
	character*48 out3file,out2file,out1file
	character infile*36,infile1*36, infile2*36	!semi-old naming convention
c	character infile*47,infile1*47, infile2*47	!new naming convention (stuff shifted right by 10)
	character doy(15)*3,ddoy*3  
c	character nname*15,outname*25
	character nname*14,outname*24
	character hhead(5)*10,dum*3
	character ddate(15)*8,dddate*8
	data dum /'  '/
c
c	2003
c
	data infile   /'G:/20030151/c2m2_2003_151_hhmmss.raw'/
	data infile1  /'G:/20030151/c2m2_2003_151_hhmmss.raw'/
	data infile2  /'G:/20030151/c2m2_2003_151_hhmmss.raw'/
c	data outname  /'c2frm_2003_DOY_hhmmss.nic'/		
c	data out1file /'G:/SEAN/buffer2/Crx/c2frm_2003_DOY_hhmmsscr1.grd'/	!cosmic-ray not removed file...
c	data out2file /'G:/SEAN/buffer2/Crx/c2frm_2003_DOY_hhmmsscrx.grd'/	!cosmic-ray removed file...
c	data out3file /'G:/SEAN/buffer2/Crx/c2avg_2003_DOY_hhmmsscrx.grd'/	!cosmic-ray removed binary file...
c
c	2004
c
c	data infile   /'G:/20040DOY/c2m2_2004_DOY_hhmmss.raw'/
c	data infile1   /'G:/20040DOY/c2m2_2004_DOY_hhmmss.raw'/
c	data infile2   /'G:/20040DOY/c2m2_2004_DOY_hhmmss.raw'/
c
c	New naming convention...
c
c	data infile    /'G:/SEAN/buffer2/Nic288/c2m2_2004_DOY_hhmmss.raw'/
c	data infile1   /'G:/SEAN/buffer2/Nic288/c2m2_2004_DOY_hhmmss.raw'/
c	data infile2   /'G:/SEAN/buffer2/Nic288/c2m2_2004_DOY_hhmmss.raw'/	!havn't got beyond here in convention chgs.
c	
c	data outname  /'c2frm_2004_DOY_hhmmss.nic'/		
c	data out1file /'G:/SEAN/buffer2/Crx/c2frm_2004_DOY_hhmmsscr1.grd'/	!cosmic-ray not removed file...
c	data out2file /'G:/SEAN/buffer2/Crx/c2frm_2004_DOY_hhmmsscrx.grd'/	!cosmic-ray removed file...
c	data out3file /'G:/SEAN/buffer2/Crx/c2avg_2004_DOY_hhmmsscrx.grd'/	!cosmic-ray removed binary file...
c	
	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',
	data ddate /'05/31/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',
c	data ddate /'04/19/04','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'/
c	data doy /'327','328','146','147','148','149','150','151','152','153','154','155','156','157','158'/
	data doy /'151','145','146','147','148','149','150','151','152','153','154','155','156','157','158'/
c	data doy /'110','145','146','147','148','149','150','151','152','153','154','155','156','157','158'/
	data badsum /0./
	data pattern /20352*0./
	data ipcnt /20352*0/
	data n /20352/
	data trash,iffix,iqcnt,isunmoon,isatcnt /5*0/
	data patcr /6.94/	!13.89 adu's = 1000 electrons
	data minresponse /1524*100000./
	data mingrit /1524*100000./
	data torbit /6096.1D0/	!value determined from "start_times.for"
	data iarglast /1524/
	data darkdist /762000*0./
	data frameavg /20352*0./
	data glaremap /20352*0./
	data glarecnt /20352*0./
	data jedge /13*64,2*63,62,2*61,2*60,2*59,2*58,57,2*56,2*55,2*54,2*53,2*52,2*51,2*50,2*49,3*48,2*47,2*46,
     &	3*45,2*44,3*43,2*42,3*41,2*40,3*39,3*38,3*37,3*36,4*35,3*34,4*33,4*32,4*31,4*30,5*29,5*28,5*27,7*26,8*25,
     &	10*24,42*23/
	data edgesave /180*0./
	data iedgesave /180*0/
	data weightsum /0./
	data edgedist /947700*0./
	data edge /180*0./
c
	icrx=0			!icrx=1 means make before and after files
	iavg=0			!iavg=1 means make average maps every 10 frames
	igraph=0		!igraph=1 means make GRAPHER output file
	mode = 4	
	id = 1272
	jd = 256
	id=id/mode
	jd=jd/mode
c			
c		!this section for reading in pattern
c
c	patdark=8.18373
c	open (13,file='c2pat_2003_326_4x4.grd',readonly)	!DOY327
c	print *,'reading in c2pat_2003_326_4x4.grd'		!DOY327
	open (13,file='c2pat_2003_144_4x4.grd',readonly)	!DOY151
	print *,'reading in c2pat_2003_144_4x4.grd'		!DOY151
c	open (13,file='c2pat_2004_092_4x4.grd',readonly)	!DOY110
c	print *,'reading in c2pat_2004_092_4x4.grd'		!DOY110
	read(13,1)hhead
1	format(a10)
	do j=1,64
	read(13,15)(pattern(i,j),i=1,318)
15	format(20f9.3)
	enddo
	read (13,16)patdark
16	format(f10.4)
	print *,'patdark = ',patdark	
	close(13)
c	open (10,file='namelist2_327.txt',readonly)		!DOY327
c	if(igraph.eq.1)open (12,file='cam2_327to328_CRX.dat')	!DOY327
	open (10,file='namelist2_151.txt',readonly)		!DOY151
	if(igraph.eq.1)open (12,file='cam2_151to151_CRX.dat')	!DOY151
c	open (10,file='namelist2_110.txt',readonly)		!DOY110
c	if(igraph.eq.1)open (12,file='cam2_110to110_CRX.dat')	!DOY110
c
	pedrecover=0
	pedinit=0.	!1014.
	pedlast=pedinit
	pedinit2=0.	!1020.
	icrcnt=0
	iavgcnt=0
	t0=+2014.5D0	!DOY144
c********************************************************************
c	kmax =  2905	!DOY327
	kmax =	3879	!DOY151
c	kmax =  1790	!DOY110
c	kmax =  6264	!DOY288
	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=270426	
	kmin=1
c	kmax=5365	!DOY151, all frames singly to the end...
	i10skip=0
	open(20,file='DOY151_del_glare2.dat')
	do k=1,kmax
	if(i10skip.gt.0)then
		read (10,10)nname,h,m,s
		i10skip=i10skip-1
		go to 98
	endif
	read (10,10)nname,h,m,s
c10	format(39x,a15,3i2.2)
10	format(39x,a14,3i2.2)
	if(k.lt.kmin)go to 98
c
c	if(k.gt.206)go to 98			!pick out glare range of interest for DOY151, 1st 
c	if(k.lt.1390.or.k.ge.1518)go to 98	!pick out glare range of interest for DOY151, 2nd
c	if(k.lt.1518.or.k.ge.1696)go to 98	!pick out glare range of interest for DOY151, 3rd
c	if(k.lt.2841.or.k.ge.3023)go to 98	!pick out glare range of interest for DOY151, 4th
	if((k.lt.1518.or.k.ge.1696).and.k.gt.206)go to 98	!pick out range for DOY151, 1st & 3rd
c	
	if(k.le.84)go to 98			!remove bad particles early in DOY151, 1st
	if(k.eq.1548)go to 98			!kill a bad frame, DOY151
c
c	if(k.lt.1021.or.k.gt.1196)go to 98	!pick out glare range of interest for DOY327, 1st
c	if(iabs(k-1131).le.2.or.iabs(k-1169).le.4)go to 98	!kill some bad frames, DOY327, 1st
c	if(k.lt.2495.or.k.ge.2623)go to 98	!pick out glare range of interest for DOY327, 2nd
c	if((k.lt.1021.or.k.gt.1196).and.(k.lt.2495.or.k.ge.2623))go to 98	!range for DOY327, 1st&2nd
c
c	if(k.lt.900.or.k.gt.1080)go to 98	!best-shot range for DOY110, from ephemeris: note Sun > 90 here...
c
c	Here remove frames which do not have a later frame to compare with
c	
c	if(k.eq.330.or.(k.ge.1348.and.k.le.1363).or.(k.ge.1797.and.k.le.1824))go to 98	!DOY327
c
	if(k.eq.35.or.k.eq.78)go to 98	!DOY151	
	if(k.eq.125.or.k.eq.170.or.k.eq.214.or.k.eq.257.or.k.eq.300.or.k.eq.344.or.k.eq.389)go to 98	!DOY151
	if(k.eq.434.or.k.eq.480.or.k.eq.525.or.k.eq.566.or.k.eq.608.or.k.eq.654.or.k.eq.700)go to 98	!DOY151
	if(k.eq.735.or.k.eq.767.or.k.eq.791.or.k.eq.821.or.k.eq.857.or.k.eq.902.or.k.eq.946)go to 98	!DOY151
	if(k.eq.989.or.k.eq.1029.or.k.eq.1074.or.k.eq.1121.or.k.eq.1168.or.k.eq.1215.or.k.eq.1262)go to 98	!DOY151
	if(k.eq.1309.or.k.eq.1355.or.k.eq.1403.or.k.eq.1448.or.k.eq.1489.or.k.eq.1531.or.k.eq.1576)go to 98	!DOY151
	if(k.eq.1622.or.k.eq.1666.or.k.eq.1710.or.k.eq.1752.or.k.eq.1795.or.k.eq.1838.or.k.eq.1877)go to 98	!DOY151
	if(k.eq.1920.or.k.eq.1966.or.k.eq.2012.or.k.eq.2046.or.k.eq.2066.or.k.eq.2087.or.k.eq.2113)go to 98	!DOY151
	if(k.eq.2123.or.k.eq.2159.or.k.eq.2224.or.k.eq.2251.or.k.eq.2288.or.k.eq.2324.or.k.eq.2367)go to 98	!DOY151
	if(k.eq.2412.or.k.eq.2456.or.k.eq.2497.or.k.eq.2540.or.k.eq.2586.or.k.eq.2633.or.k.eq.2681)go to 98	!DOY151
	if(k.eq.2729.or.k.eq.2776.or.k.eq.2822.or.k.eq.2867.or.k.eq.2913.or.k.eq.2945.or.k.eq.2964)go to 98	!DOY151
	if(k.eq.2987.or.k.eq.3028.or.k.eq.3074.or.k.eq.3119.or.k.eq.3163.or.k.eq.3207.or.k.eq.3249)go to 98	!DOY151
	if(k.eq.3289.or.k.eq.3322.or.k.eq.3348.or.k.eq.3378.or.k.eq.3418.or.k.eq.3464.or.k.eq.3509)go to 98	!DOY151
	if(k.eq.3554.or.k.eq.3593.or.k.eq.3636.or.k.eq.3637.or.k.eq.3682.or.k.eq.3729.or.k.eq.3768)go to 98	!DOY151
	if(k.eq.3800.or.k.eq.3830.or.k.eq.3871)go to 98	!DOY151
c
c	if(k.eq.6)go to 98	!DOY110
c
	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
c	write(infile1(13:27),'(a15)')nname
c	write(infile1(28:33),'(3i2.2)')h,m,s
	write(infile1(27:32),'(3i2.2)')h,m,s	!DOY151
c	write(outname(1:15),'(a15)')nname
c	write(outname(16:21),'(3i2.2)')h,m,s
	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
c	write(infile2(9:11),'(a3)')ddoy
c	write(nname(12:14),'(a3)')ddoy
c	write(infile2(13:27),'(a15)')nname
c	write(infile2(28:33),'(3i2.2)')hh,mm,ss
	write(infile2(27:32),'(3i2.2)')hh,mm,ss	!DOY151
c	if(k.ge.0)print *,k,dum,infile1,dum,infile2
	open(11,file=infile1,form='binary',readonly)
	read(11)head,((iframe1(i,j),i=1,318),j=1,64)	!,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,318),j=1,64)	!,trail1,q,trail2a,camtemp2,trail2b,tdum1,trail3,tdum2,trail4
	close(11)
c	if(trail1(15).eq.0)then		!open-shutter data only
c	if(camtemp.ne.999.99)then
c		if((t-tprev),gt.90.)then
c			print *,infile,',wind off 11 frames, time diff = ',t-tprev
c			i10skip=10
c			go to 98
c		endif
c	else
c		camtemp=camtempsave
c	endif
c	endif				!open-shutter data only
	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	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(iframe1(i,j).eq.-16389)isat=isat+1
	if(iframe2(i,j).eq.-16389)isat=isat+1
	enddo
	enddo
	if(isat.ne.0)then
		isatcnt=isatcnt+1
		print *,' frame ',k,': saturated pixels...'
		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,64
	do i=1,318
	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,318
	do j=1,64
		iarg=int4(iframe2(i,j))
		if(iarg.lt.0)iarg=iarg+66536
		iframe(i,j)=int4(iframe1(i,j))
		if(iframe(i,j).lt.0)iframe(i,j)=iframe(i,j)+65536
		iframe(i,j)=iframe(i,j)-iarg
	enddo
	enddo
c	
c**********************************!	Pedestal calculation
c
	do j=1,64
	if(iframe(4,j).ne.100)then
		iarg=iframe(4,j)
		arg=float(iarg)/0.75
		pednaive=pednaive+arg
		ipedcnt=ipedcnt+1
		pedsv(ipedcnt)=arg
	endif
	if(iframe(318,j).ne.100.and.j.ne.64)then
		iarg=iframe(318,j)
		arg=float(iarg)/0.75
		pednaive=pednaive+arg
		ipedcnt=ipedcnt+1
		pedsv(ipedcnt)=arg
	endif
	enddo
	pednaive=pednaive/127.
	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	
c
	do j=1,127
	    if(pedsv(indxx(j))-pedtest.gt.2.)go to 20	!exclude high/low pedestals above 144 electrons...
	    if(pedsv(indxx(j))-pedtest.gt.-2.)then
	    	ipedcnt=ipedcnt+1
	    	pedsum=pedsum+pedsv(indxx(j))
	    endif
	enddo
20	if(ipedcnt.lt.100)ibadped=1		!112 = 1 frame, scaled from eng. data
c	if(pedsv(indxx(1))-pedtest.lt.-10..and.pedsv(indxx(2))-pedtest.lt.-10.)then
c		crazyped1=pedsv(indxx(1))
c		crazyped2=pedsv(indxx(2))
c		print*,infile,': BAD PEDESTAL VALUES = ',crazyped1,crazyped2,pedtest
c		ibadped=2
c	endif
	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,64
	if(iframe(5,j).ne.100)then
		iarg=iframe(5,j)
		arg=float(iarg)/0.75-ped
		darksum=darksum+arg
	endif
	if(iframe(317,j).ne.100.and.j.ne.64)then
		iarg=iframe(317,j)
		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.100)then
		iarg=iframe(5,j)
		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
		if(k-kmin+1.le.6000.and.k-kmin+1.gt.0)darkdist(jarg,k-kmin+1)=ddark
	endif
	enddo
	do j=1,63
	if(iframe(317,j).ne.100)then
		iarg=iframe(317,j)
		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
		if(ddark.gt.patcr)then
		    idarkcnt=idarkcnt-1
		    darksum=darksum-arg
		    darkmed(jarg)=10000.
		    ddark=0.
		endif
		if(k-kmin+1.le.6000.and.k-kmin+1.gt.0)darkdist(jarg,k-kmin+1)=ddark
	endif
	enddo
	if(idarkcnt.le.10)then
		ibaddark=1
		go to 53
	endif
c	if(k-kmin+1.le.6000.and.k-kmin+1.gt.0)then
c	do j=1,127
c		if(darkmed(j).ne.0..and.darkmed(j).ne.10000.)darkdist(j,k-kmin+1)=darkmed(j)
c	enddo
c	endif
c
c	Get median:"interpolate within the histogram..."
c
	call indexxx(127,darkmed,indxx)
c	if(k-kmin+1.le.6000)then	
c	do j=1,127
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,127-imed
	if(darkmed(indxx(imed+j)).ne.dmed)then
		med2=imed+j
		go to 52
	endif
	enddo
52	continue
	median=darkmed(indxx(imed))+(float(imed-med1)/float(med2-med1)-0.5)/0.75	!pretty tricky...
	if(idarkcnt.lt.100)ibaddark=1
	dark=darksum/float(idarkcnt)
c	print 2000,k,darknaive,dark,median,idarkcnt	!,dark/median
2000	format(i6,3f10.2,i6)
c	dark=median	
c
c	Now correct time-plotted covered pixels by best dark current scale factor	
c
c	if(k-kmin+1.le.6000.and.k-kmin+1.gt.0)then
c	patcorr=(dark/patdark)**1.5
c	do j=1,127
c		darkdist(j,k-kmin+1)=darkdist(j,k-kmin+1)-darksv(j)*patcorr
c	enddo
c	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
		dark=median
	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,a37)
			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,a37)
		endif
		ped=-9999.
		dark=-9999.
		if(k-kmin+1.le.6000.and.k-kmin+1.gt.0)then
		do j=1,127
			darkdist(j,k-kmin+1)=0.
		enddo
		endif
		go to 98
	endif
c
c	if(k.ge.1)go to 97	!$$$$$	bail quick	
c
	pixsum=0.
	ipixsum=0
	pixdif=0.
	ipixdif=0	
	patcorr=dark/patdark	!was raised to the 1.5 power for engineering data
c	patcorr=0.	!*******************
	do j=1,64
	do i=1,318
		iarg=iframe(i,j)
		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
		fframe(i,j)=abs(arg)
c		if(iarg.eq.0)frame(i,j)=0.
	enddo
	enddo
	temp1=0.
	itemp=0
	do i=14,180
	j=jedge(i)
			if(frame(i,j).lt.125.)then
				temp1=temp1+frame(i,j)
				itemp=itemp+1
				edgesave(i)=edgesave(i)+frame(i,j)
				iedgesave(i)=iedgesave(i)+1
				edge(i)=frame(i,j)
				edgedist(k,i)=edge(i)
			endif
	enddo
	if(itemp.lt.30)then
		print 1002,k,itemp,infile1
1002		format('BAD GLARE',i6,', only',i3,' below 125ADU,'1x,a37)
		go to 98
	endif
	call index180(180,edge,indx180)
	imed=ifix(float(itemp)/2.0)
	tempmean=edge(indx180(imed))
	do i=1,180
c		edgedist(k,i)=edge(indx180(i))
	enddo
	temp1=temp1/float(itemp)
c	if(mod(k,1000).eq.0)print *,k
c	if(k.ne.0)go to 98	!bail quick
	weight=abs(temp1/29.)		!the triangular eclipse depth below gives slightly less glare...
c	print *,k,ped,dark,meascnt,itemp,weight
c	weight=1.
c	
	call INDEXX(n,fframe,indx)
	call measles(n,frame,indx,meascnt,ineigh,ihist,ikill)
	iavgcnt=iavgcnt+1
	do j=1,64
		do i=1,318
			if(iframe1(i,j).eq.0.and.iframe2(i,j).eq.0)frame(i,j)=0.	!impose ROI
			if(ikill(i,j).eq.1)frame(i,j)=0.	!KILL measles pixels
			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.ge.6.and.i.lt.317)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.ge.6.and.i.lt.317)then
					ipixdif=ipixdif+1
					pixdif=pixdif+abs(frame(i,j-1)-frame(i,j))
				endif
			endif
		enddo	
	enddo
c
	pixsum=pixsum/float(ipixsum)
	pixdif=pixdif/float(ipixdif)
	if(meascnt.gt.20)go to 98	!this cut was at 100 for DOY151-2 and 4
	if(k.lt.1518)then		!triangular eclipse depth, not quite right.
		weight=1.-float(iabs(k-90))/116.
	else
		weight=1.-float(1585-k)/110.
		if(weight.gt.1.)weight=1.-float(k-1585)/116.
	endif
	write(20,56)k,temp1/29,tempmean,ped-pedinit,dark,meascnt,weight
56	format(i6,4f10.2,i6,f10.2)
	do i=1,318
	do j=1,64
		if(iframe1(i,j).ne.0)then
		if(k.ge.1518)then	!if(k.lt.10000)then	!if(k.ge.1518)then	!special for DOY151
			if(frame(i,j).lt.15..and.frame(i,j).gt.-75.)then	!75 for DOY151, 15 for DOY327
				glaremap(i,j)=glaremap(i,j)-frame(i,j)*weight
				glarecnt(i,j)=glarecnt(i,j)+weight
			endif
		else
			if(frame(i,j).gt.-15..and.frame(i,j).lt.75.)then
				glaremap(i,j)=glaremap(i,j)+frame(i,j)*weight
				glarecnt(i,j)=glarecnt(i,j)+weight
			endif
		endif
		endif
	enddo
	enddo
c	
c	if(k.ge.1)go to 98	
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
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,64
c		do i=6,316
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)write(12,120)dddate,h,m,s,ped-pedinit,dark,meascnt,pixsum,pixdif,camtemp1-camtemp2	!exp((camtemp+10.)/8.7)
c	if(mod(k,1000).eq.0)then
	print 12,dddate,h,m,s,ped-pedinit,dark,meascnt,pixsum,pixdif,k
12	format('"',a8,1x,i2.2,':',i2.2,':',i2.2,'"',2x,2f10.3,i6,2f10.2,i7)
120	format('"',a8,1x,i2.2,':',i2.2,':',i2.2,'"',2x,2f10.3,i6,3f10.2)
33	format(a25,a3,2f10.3,4f13.9,i6,2f10.2)
c	endif
	write(out1file(32:41),'(a10)')infile1(24:33)
	write(out2file(32:41),'(a10)')infile1(24:33)
	write(out3file(32:41),'(a10)')infile1(24:33)
	if(iavg.eq.1.and.iavgcnt.eq.10)then
	open (16,file=out3file)
	write(16,140)
	write(16,150)((frameavg(i,j),i=1,318),j=1,64)
	close(16)
	iavgcnt=0
	do j=1,64
	do i=1,318
	frameavg(i,j)=0.
	enddo
	enddo
	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'/'318 64'/'0 318'/'0 64'/'-50 50')
	if(icrx.eq.1)write(15,150)((frameout(i,j),i=1,318),j=1,64)
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)
c	
c	open (17,file='min2_4.dat')
c	do i=1,1524
c	write(17,17)i,minresponse(i)
c	enddo
c	close(17)	
c	open (17,file='mingrit2_4.dat')
c	do i=1,1524
c	write(17,17)i,mingrit(i)
c	enddo
c	close(17)
c
	close(20)
	print *,'Start writing output files:'
	print *,isatcnt,' discarded as sun/moon based on saturated pixel(s)'
	print *,isunmoon,' discarded as sun/moon based on raw/crazy pedestal'
	print *,trash,' others ped/dark discarded out of ',kmax-kmin+1
	if(igraph.eq.1)then
	print *,'It aint done until the fat lady sings...'
	open(10,file='edgedist.grd')
	write(10,13)
13	format('DSAA'/'5265 180'/'0 5265'/'0 1600'/'0 125')
	do j=1,180
	do i=1,5265
	if(edgedist(i,j).gt.999.99)edgedist(i,j)=999.99
	if(edgedist(i,j).lt.-99.99)edgedist(i,j)=-99.99
	enddo
	write(10,130)(edgedist(i,j),i=1,5265)	
130	format(20f7.2)	
	enddo
	close(10)
	endif
c
c	form, write output glare map
c
c	if(igraph.eq.100)then	
	open(10,file='d:\sean\eclipse\glare2_151_1_el_depth.grd')
	write(10,100)
100	format('DSAA'/'318 64'/'0 318'/'0 64'/'-1 30')
	do j=1,64
	do i=1,318
		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,318),j=1,64)
110	format(20f9.2)	
	close(10)
	open(11,file='edge.dat')
	do i=1,180
	if(iedgesave(i).gt.0)edgesave(i)=edgesave(i)/iedgesave(i)
	write(11,11)i,edgesave(i)
11	format(i4,f10.2)
	enddo	
c	
c	endif
	print *,'done'
	stop 
	end
c********************************************************************************************
	subroutine measles(n,frame,indx,meascnt,ineigh,ihist,ikill)
	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(20352),ikill(318,64)
	real*4 frame(318,64),z(3,3),z0,z0sv
	real*4 dcut,sumzsq,devcut
	
	data dcut /13.88/				!13.88 is a 1000-electron cut
	data devcut /400./				!20 standard deviation cut
	
	meascnt=0

	do m=1,n
		iarg=1+n-m
		j=(indx(iarg)-1)/318
		i=indx(iarg)-318*j
		j=j+1
		if(i.le.6.or.i.ge.316.or.j.lt.2.or.j.gt.63)go to 98	!bail if near edges
				
c		if(frame(i,j).ne.0.and.frame(i,j).ne.-10..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(20352),N,I,J,L,IR,INDXT	!Absoft Compiler dislikes INDX(N), etc.
      real*4 ARRIN(20352),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(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            
c*********************************************************************************
      SUBROUTINE INDEX180(N,ARRIN,INDX)	!from Numerical Receipes...
      implicit none
      integer*4 INDX(180),N,I,J,L,IR,INDXT	!Absoft Compiler dislikes INDX(N), etc.
      real*4 ARRIN(180),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