program Cam3Average
	implicit none
c
c	This program averages together camera 3 maps
c	This is a square 180 x 180 degrees 
c
	character*34 polething
	character*46 outfile1,outfile2,outfile4
	character*47 outfile3(52)
	character*58 filename1,filename2,blankk
	character*62 filename3,blank3
	character*10 head(5)
	character*12 XY
	character*3  DOY
	character*1 i100,i10,i1
	integer*4 i,ii,imax,j,k,idoy,nmap(360,360),icount(360,360,52),iarg,imap,iweek,resmap(200,52),iarg1,iidoy,iiyr,iiweek
	real*4 frame(360,360),xsq,arg,arg1,arg2,arg3,arg4,eps,omega,phi,dSMEI
	real*4 polemap(360,360),amap(360,360),corr32,corr3,cmap(360,360,52),rDis,rDismin,rLat,rLng,PA	!,zodi0,c,s,betabar
	data omega /78.25/
	data rDismin / 0.983292 /
	data blankk /'                                                          '/
	data blank3 /'                                                              '/
	data icount /6739200*0/
	data cmap  /6739200*0./
	data resmap /10400*0/
	data imax /1947/
c	
c	imap = 0 means make the average map, don't output weekly averages
c	imap = 1 means read the average map, do make & output weekly averages
c	imap = 2 means read the average map, but make & output only the one "iiweek" average
c
	data imap /1/
c
	data iiweek / 39 /	
c
c	data outfile1 /'g:\GegenscheinXStars\c3avg_2003through2008.grd'/	!use this when imap = 0
c	data outfile1 /'g:\GegenscheinXStars\c3avg_2003-8rolledoff.grd'/	!use this when imap = 1
c
c	Note here we bring in the rolled off average map ...
c	When this program is run with imap = 0, don't forget to remake the rolled off map with Sun_Fisheye.for
c	before running this program with imap = 1...
c	
	data polething /'g:\GegenscheinXStars\polething.grd'/
	data outfile2 /'g:\GegenscheinXStars\n3avg_2003through2008.grd'/
	data outfile4 /'g:\GegenscheinXStars\resmap2003through2008.grd'/
c	
	data outfile3(1)  /'g:\GegenscheinXStars\cmap01_2003through2008.grd'/
	data outfile3(2)  /'g:\GegenscheinXStars\cmap02_2003through2008.grd'/
	data outfile3(3)  /'g:\GegenscheinXStars\cmap03_2003through2008.grd'/
	data outfile3(4)  /'g:\GegenscheinXStars\cmap04_2003through2008.grd'/
	data outfile3(5)  /'g:\GegenscheinXStars\cmap05_2003through2008.grd'/
	data outfile3(6)  /'g:\GegenscheinXStars\cmap06_2003through2008.grd'/
	data outfile3(7)  /'g:\GegenscheinXStars\cmap07_2003through2008.grd'/
	data outfile3(8)  /'g:\GegenscheinXStars\cmap08_2003through2008.grd'/
	data outfile3(9)  /'g:\GegenscheinXStars\cmap09_2003through2008.grd'/
	data outfile3(10) /'g:\GegenscheinXStars\cmap10_2003through2008.grd'/
	data outfile3(11) /'g:\GegenscheinXStars\cmap11_2003through2008.grd'/
	data outfile3(12) /'g:\GegenscheinXStars\cmap12_2003through2008.grd'/
	data outfile3(13) /'g:\GegenscheinXStars\cmap13_2003through2008.grd'/
	data outfile3(14) /'g:\GegenscheinXStars\cmap14_2003through2008.grd'/
	data outfile3(15) /'g:\GegenscheinXStars\cmap15_2003through2008.grd'/
	data outfile3(16) /'g:\GegenscheinXStars\cmap16_2003through2008.grd'/
	data outfile3(17) /'g:\GegenscheinXStars\cmap17_2003through2008.grd'/
	data outfile3(18) /'g:\GegenscheinXStars\cmap18_2003through2008.grd'/
	data outfile3(19) /'g:\GegenscheinXStars\cmap19_2003through2008.grd'/
	data outfile3(20) /'g:\GegenscheinXStars\cmap20_2003through2008.grd'/
	data outfile3(21) /'g:\GegenscheinXStars\cmap21_2003through2008.grd'/
	data outfile3(22) /'g:\GegenscheinXStars\cmap22_2003through2008.grd'/
	data outfile3(23) /'g:\GegenscheinXStars\cmap23_2003through2008.grd'/
	data outfile3(24) /'g:\GegenscheinXStars\cmap24_2003through2008.grd'/
	data outfile3(25) /'g:\GegenscheinXStars\cmap25_2003through2008.grd'/
	data outfile3(26) /'g:\GegenscheinXStars\cmap26_2003through2008.grd'/
	data outfile3(27) /'g:\GegenscheinXStars\cmap27_2003through2008.grd'/
	data outfile3(28) /'g:\GegenscheinXStars\cmap28_2003through2008.grd'/
	data outfile3(29) /'g:\GegenscheinXStars\cmap29_2003through2008.grd'/
	data outfile3(30) /'g:\GegenscheinXStars\cmap30_2003through2008.grd'/
	data outfile3(31) /'g:\GegenscheinXStars\cmap31_2003through2008.grd'/
	data outfile3(32) /'g:\GegenscheinXStars\cmap32_2003through2008.grd'/
	data outfile3(33) /'g:\GegenscheinXStars\cmap33_2003through2008.grd'/
	data outfile3(34) /'g:\GegenscheinXStars\cmap34_2003through2008.grd'/
	data outfile3(35) /'g:\GegenscheinXStars\cmap35_2003through2008.grd'/
	data outfile3(36) /'g:\GegenscheinXStars\cmap36_2003through2008.grd'/
	data outfile3(37) /'g:\GegenscheinXStars\cmap37_2003through2008.grd'/
	data outfile3(38) /'g:\GegenscheinXStars\cmap38_2003through2008.grd'/
	data outfile3(39) /'g:\GegenscheinXStars\cmap39_2003through2008.grd'/
	data outfile3(40) /'g:\GegenscheinXStars\cmap40_2003through2008.grd'/
	data outfile3(41) /'g:\GegenscheinXStars\cmap41_2003through2008.grd'/
	data outfile3(42) /'g:\GegenscheinXStars\cmap42_2003through2008.grd'/
	data outfile3(43) /'g:\GegenscheinXStars\cmap43_2003through2008.grd'/
	data outfile3(44) /'g:\GegenscheinXStars\cmap44_2003through2008.grd'/
	data outfile3(45) /'g:\GegenscheinXStars\cmap45_2003through2008.grd'/
	data outfile3(46) /'g:\GegenscheinXStars\cmap46_2003through2008.grd'/
	data outfile3(47) /'g:\GegenscheinXStars\cmap47_2003through2008.grd'/
	data outfile3(48) /'g:\GegenscheinXStars\cmap48_2003through2008.grd'/
	data outfile3(49) /'g:\GegenscheinXStars\cmap49_2003through2008.grd'/
	data outfile3(50) /'g:\GegenscheinXStars\cmap50_2003through2008.grd'/
	data outfile3(51) /'g:\GegenscheinXStars\cmap51_2003through2008.grd'/
	data outfile3(52) /'g:\GegenscheinXStars\cmap52_2003through2008.grd'/
c
c	Read in "polething" map, made by program polething.for
c	
	print *,'reading in pole map ',polething
	open(12,file=polething,readonly)
		read(12,1)head
		do k=1,360
			read(12,120)(polemap(j,k),j=1,360)
120			format(20f8.2)			
		enddo
		close(12)
c
c	Read in or zero average map
c
	if(imap.eq.0)then
		outfile1 ='g:\GegenscheinXStars\c3avg_2003through2008.grd'
		do k=1,360
		do j=1,360
			amap(j,k)=0.
			nmap(j,k)=0
		enddo
		enddo
	else
		outfile1 ='g:\GegenscheinXStars\c3avg_2003-8rolledoff.grd'
		print *,'reading in average map ',outfile1
		open(12,file=outfile1)
		read(12,1)head
		do k=1,360
			read(12,11)(amap(j,k),j=1,360)
		enddo
		close(12)
	endif
c
	open (20,file='G:\GegenscheinXStars\zodi_dist.txt')	!this was made by program zodi_dist_filemake.for
	open (10,file='G:\GegenscheinXStars\filenames123.txt',readonly)
	do i=1,imax
	if(mod(i,100).eq.0)print *,i
		read(10,10)filename1,filename2,filename3
10		format(6x,a58,9x,a58,9x,a62)
		if(filename3.eq.blank3)go to 99
c		
c		if(i.ne.999)go to 99			!***************pick out a single frame**************
c
c	Stamp out bad sections
c		
c		if(i.ge. 658.and.i.le. 654)go to 99	!Kill DOYs   47-53, 2005 grooves+white,  recovered, bad pattern
c		if(i.ge. 943.and.i.le. 945)go to 99	!Kill DOYs   26-28, 2006 suddenly white, recovered, bad pattern
c		if(i.ge.1415.and.i.le.1419)go to 99	!Kill DOYs 235-240, 2007 bad groove,     recovered, bad pattern
c		if(i.eq.1481)go to 99			!Kill DOY      303, 2007 mask-out in midst, now recovered
c		if(i.ge.1710.and.i.le.1722)go to 99	!Kill DOYs 206-218, 2008 "bad patterns", recovered
c
		if(i.ge. 933.and.i.le. 937)go to 99	!Kill DOYs   13-17, 2006 suddenly black
		if(i.eq. 950)go to 99			!Kill DOY       33, 2006 bad grooves
		if(i.eq.1034)go to 99			!Kill DOY      155, 2006 very little sky coverage
		if(i.ge.1302.and.i.le.1305)go to 99	!Kill DOYs   73-76, 2007 some black
		if(i.eq.1435)go to 99			!Kill DOY      256, 2007 suddenly white
		if(i.ge.1581.and.i.le.1591)go to 99	!Kill DOYs   43-53, 2008 very black
		if(i.eq.1641.or .i.eq.1642)go to 99	!Kill DOYs 121,122, 2008 "hot"
		if(i.ge.1737.and.i.le.1748)go to 99	!Kill DOYs 233-244, 2008 black
c		
		write(XY(1:12),'(a12)')filename3(47:58)
		write(DOY(1:3),'(a3)')XY(3:5)
c	Decoding DOY and evaluate empirical DOY-dependent offset for camera 3		
		write(i100,'(a1)')XY(3:3)
		write(i10,'(a1)')XY(4:4)
		write(i1,'(a1)')XY(5:5)
		idoy=100*(ichar(i100)-48)+10*(ichar(i10)-48)+ichar(i1)-48
		iweek=1+idoy/7
			if(imap.eq.2.and.iweek.ne.iiweek) go to 99	!does 1 week at a time
		if(iweek.gt.52)iweek=52
		phi=float(idoy-1)*360./365.25
		if(phi.gt.360.)phi=360.
		arg=phi-omega+90.
		if(arg.gt.180.)arg=arg-360.
		arg=abs(arg)
		xsq=float(idoy)*float(idoy)
		corr3 =6.7*cos(float(idoy)/55.+2.4+0.00002*xsq)+1.73-0.03*float(idoy)	!flattens cam 3 alone
		corr32=7.8*cos(float(idoy)/55.+2.4+0.00002*xsq)-2.73-0.03*float(idoy)	!matches cams 2 & 3
c
c	Include ad hoc mask-in corrections: ***************depends upon sequence number in list********************
c
		if(i.gt. 972.and.i.le. 996)corr32=corr32     -0.5*float(i-966)
		if(i.gt.1136.and.i.le.1250)corr32=corr32    +0.22*float(i-1136)
		if(i.gt.1250.and.i.le.1325)corr32=corr32+28.-0.04*float(i-1250)
		if(i.gt.1434.and.i.le.1655)then
			if(i.lt.1581.or.i.gt.1591)corr32=corr32+52.-0.15*float(i-1435)
		endif
		if(i.ge.1750)then
			if(i.le.1901.or.i.ge.1915)corr32=corr32+70.
		endif
c
c	Read files, get rDis
c
		do ii=1,400
c	write(11,11)adoy,rDis,rLng,rLat,dSMEI,iyr-2003
c11	format(1x,a3,f9.6,1x,3f10.4,1x,i3)
			read(20,20)iidoy,rDis,dSMEI,iiyr
20			format(i4,f9.6,21x,f10.4,i4)
			if(iidoy.eq.idoy)go to 21
		enddo
21		continue
c		print *,filename3,ii,iidoy,dSMEI
		iiyr=iiyr+2003
		rDis=rDis/rDismin
c
		open(11,file=filename3,readonly)
3		read(11,1)head
1		format(a10)			
		do k=1,360
			read(11,11)(frame(j,k),j=1,360)
11			format(20f8.1)				
		enddo
		close(11)
c
c	Form Average Map
c
		if(imap.eq.0)then
		     do k=1,360
			rLat=float(k-180)/2.
			arg2=cosd(rLat)
			do j=1,360
				rLng=float(j-180)/2.
				arg1=cosd(rLng)
				eps=acosd(arg1*arg2)-20.
				if(eps.lt.0.)eps=0.
				PA=atan2d(sind(rLng),tand(rLat))
				arg3=0.
				if(arg.lt.45..or.arg.gt.135.)then
				   arg3=arg3+200.*exp(-eps/10.)*exp(-abs(arg2)/10.) !exponential in both elongation & latitude
				   arg3=arg3*cosd(2.*arg)
				   if(arg.gt.135.)arg3=1.5*arg3
				endif
				arg3=arg3-250.*exp(-eps/10.)*(-0.3+dSMEI/15.)*((cosd(PA))**2)	!Distance Sun-to-SMEI axis adjustment, feathered to South
				arg3=arg3-60.0*exp(-eps/15.)*exp(-abs(float(iarg-41)/4.))	!Adjustment when Sun closest: Fall mostly feathered
				arg3=arg3-40.0*exp(-eps/5.)*exp(-abs(float(iarg-41)/4.))	!Adjustment when Sun closest: Fall mostly feathered
				if(k.ge.190)arg3=0.
				if(k.ge.170.and.k.lt.190)arg3=arg3*(float(190-k)/20.)		!feather along latitude: mostly to South
				arg3=arg3-35.*exp(-eps/35.)*exp(-abs(float(iarg-11)/4.))	!Adjustment when Sun closest: Spring not feathered
				arg3=arg3-10.*exp(-eps/5.)*exp(-abs(float(iarg-9)/4.))		!Adjustment when Sun closest: Spring not feathered
				arg3=arg3-5.0*exp(-eps/40.)*exp(-abs(float(iarg-41)/4.))	!Adjustment when Sun closest: This adds some in the North
				arg3=arg3+40.*exp(-eps/5.)*exp(-abs(float(iarg-41)/4.))		!Adjustment when Sun closest: This reverses 2nd term above
				arg3=arg3-0.0	!was -4.0
				arg4=frame(j,k)-corr32-arg3-polemap(j,k)			!Take away corr32, arg3 & polemap
				if(frame(j,k).ne.0..and.abs(arg4).lt.100.)then
					amap(j,k)=amap(j,k)+arg4
					nmap(j,k)=nmap(j,k)+1
				endif
			enddo
		     enddo
		else
c
c	Count & average maps contributing each week (DOYs 265,6 are added to week 52)
c
			iarg=iweek
			if(iarg.gt.52)iarg=52
			do k=1,360
				rLat=float(k-180)/2.
				arg2=cosd(rLat)
c				betabar=1.5*(sqrt(1.+(rLat/1.5)**2)-1.)
				do j=1,360
					rLng=float(j-180)/2.
					arg1=cosd(rLng)
					eps=acosd(arg1*arg2)-20.
					if(eps.lt.0.)eps=0.
					PA=atan2d(sind(rLng),tand(rLat))
c					c=cosd(eps)
c					s=sind(eps)
c					zodi0=(65.+120.*c-185.*c*c+65.*c*c*c)*(10.**(-sind(betabar/(0.009*(eps+60.)))))*(s**(-2.3))
					arg3=0.
					if(arg.lt.45..or.arg.gt.135.)then
					   arg3=arg3+200.*exp(-eps/10.)*exp(-abs(arg2)/10.) !exponential in both elongation & latitude
					   arg3=arg3*cosd(2.*arg)
					   if(arg.gt.135.)arg3=1.5*arg3
					endif
c					
c					arg3=arg3-zodi0*(1.-rDis**0.2)
c					arg3=arg3-400.*exp(-eps/10.)*(rDis**0.2)*((cosd(PA))**2)	  !Distance-to-Sun adjustment, feathered to South
c					arg3=arg3-200.*exp(-eps/10.)*(1.-sqrt(cosd(phi)))*((cosd(PA))**2) !Distance-to-Sun adjustment, feathered to South
c
					arg3=arg3-250.*exp(-eps/10.)*(-0.3+dSMEI/15.)*((cosd(PA))**2)	!Distance-Sun-to-SMEI axis adjustment, feathered to South
					arg3=arg3-60.0*exp(-eps/15.)*exp(-abs(float(iarg-41)/4.))	!Adjustment when Sun closest: Fall mostly feathered
					arg3=arg3-40.0*exp(-eps/5.)*exp(-abs(float(iarg-41)/4.))	!Adjustment when Sun closest: Fall mostly feathered
					if(k.ge.190)arg3=0.
					if(k.ge.170.and.k.lt.190)arg3=arg3*(float(190-k)/20.)		!feather along latitude: mostly to South
					arg3=arg3-35.*exp(-eps/35.)*exp(-abs(float(iarg-11)/4.))	!Adjustment when Sun closest: Spring not feathered
					arg3=arg3-10.*exp(-eps/5.)*exp(-abs(float(iarg-9)/4.))		!Adjustment when Sun closest: Spring not feathered
					arg3=arg3-5.0*exp(-eps/40.)*exp(-abs(float(iarg-41)/4.))	!Adjustment when Sun closest: This adds some in the North
					arg3=arg3+40.*exp(-eps/5.)*exp(-abs(float(iarg-41)/4.))		!Adjustment when Sun closest: This reverses 2nd term above
					arg4=frame(j,k)-corr32-amap(j,k)-arg3-polemap(j,k)		!Take away polemap, residue map, & the other stuff
					if(frame(j,k).ne.0..and.abs(arg4).lt.100.)then
						icount(j,k,iarg)=icount(j,k,iarg)+1
						cmap(j,k,iarg)=cmap(j,k,iarg)+arg4
						iarg1=nint(arg4+100.5)
						if(iarg1.lt.1)iarg1=1
						if(iarg1.gt.200)iarg1=200
						resmap(iarg1,iarg)=resmap(iarg1,iarg)+1
					endif
				enddo
			enddo
		endif
99	continue
	enddo
c
	if(imap.eq.0)then	
		do k=1,360
		do j=1,360
			if(nmap(j,k).gt.0)amap(j,k)=amap(j,k)/float(nmap(j,k))
			if(nmap(j,k).lt.10)amap(j,k)=0.
		enddo
		enddo
c
		print *,'writing ',outfile1		
		open(12,file=outfile1)
		open(13,file=outfile2)
		write(12,12)
		write(13,13)
12		format('DSAA'/'360 360'/'-90 90'/'-90 90'/'-100 100')
13		format('DSAA'/'360 360'/'-90 90'/'-90 90'/'0 2000')
		do k=1,360
			write(12,11)(amap(j,k),j=1,360)
			write(13,14)(nmap(j,k),j=1,360)
14			format(20i5)		
		enddo
		close(12)
		close(13)
	else
		if(imap.eq.1)then
		open(15,file=outfile4)
		write(15,15)
15		format('DSAA'/'200 52'/'-100 100'/'0 52'/'0 100000')
		endif		
		do iarg=1,52
		if(imap.eq.2.and.iarg.ne.iiweek)go to 101
		if(imap.eq.1)write(15,150)(resmap(j,iarg),j=1,200)
150		format(20i7)		
		do k=1,360
		do j=1,360
			if(icount(j,k,iarg).gt.1)cmap(j,k,iarg)=cmap(j,k,iarg)/float(icount(j,k,iarg))
			if(icount(j,k,iarg).lt.10)cmap(j,k,iarg)=0.
			if(abs(cmap(j,k,iarg)).gt.9999.)cmap(j,k,iarg)=9999.
		enddo
		enddo
c	
		open (12,file=outfile3(iarg))
		write(12,12)
		do k=1,360
			write(12,11)(cmap(j,k,iarg),j=1,360)
		enddo
		close(12)
101		enddo
		if(imap.eq.1)close(15)
	endif
100	stop 'done'
	end