	Program cam3glaremap
	implicit none
c
c	This program combines 4 camera 3 glare maps to make one
c	
	real*4 glaremap(4,636,128),sum,glaremapall(636,128),glaremapallpolar(636,128)
	real*4 x,y,r,rm,theta,arg,arg1,glaremapmodel(636,128)
	integer*4 i,j,k,isum,jr,itheta,mapmakeall,ii,jj
	character*11 head(5),headsv(5)
	data glaremapall /81408*0./
	data glaremapallpolar /81408*0./
	data glaremapmodel /81408*0./
	mapmakeall=0	!1 means remake "all" map
c	
c
1	format(a11)
2	format(20f9.2)
	if(mapmakeall.eq.1)then	
	    do k=1,4
		print *,k
		if(k.eq.1)open(10,file='d:\sean\eclipse\glare3_288_2nd_e.grd',readonly)
		if(k.eq.2)open(10,file='d:\sean\eclipse\glare3_288_2nd_l.grd',readonly)
		if(k.eq.3)open(10,file='d:\sean\eclipse\glare3_288_3rd_e.grd',readonly)	
		if(k.eq.4)open(10,file='d:\sean\eclipse\glare3_288_3rd_l.grd',readonly)
		read(10,1)head
		print *,head
			read(10,2)((glaremap(k,i,j),i=1,636),j=1,128)
		close(10)
	    enddo
	    print *,'input files read in'
c
	    do i=1,636
	    do j=1,128
		sum=0.
		isum=0
		do k=1,4
			if(k.le.2)then
				sum=sum+glaremap(k,i,j)
				isum=isum+1
			else
				sum=sum+1.37*glaremap(k,i,j)
				isum=isum+1
			endif
		enddo
		glaremapall(i,j)=sum/float(isum)	
    	    enddo
	    enddo
c	
	    open (11,file = 'd:\sean\eclipse\glare3_288_all.grd')
	    write(11,1)(head(k),k=1,5)
	    write(11,2)((glaremapall(i,j),i=1,636),j=1,128)
	    close(11)
c	    
	else
	    open(10,file='d:\sean\eclipse\glare3_288_all.grd',readonly)
	    read(10,1)head
	    read(10,2) ((glaremapall(i,j),i=1,636),j=1,128)
	    close(10)
	endif
	do i=1,5
		headsv(i)=head(i)
	enddo
c	
	do j=1,128
	do i=1,636
c	glaremapallsmoo(i,j)=glaremapall(i,j)
		x=float(i-320)
		y=float(625-j)
		theta=atan2d(x,y)
		itheta=nint(317.+theta*10.)
		r=sqrt(x*x+y*y)
		rm=r-575.
		jr=nint(rm)
		if(glaremapall(i,j).ne.0.)then
			if(jr.lt.1)jr=1
			if(jr.gt.128)jr=128
			if(itheta.lt.1)print *,i,j,theta,itheta
			if(itheta.lt.1)itheta=1
			if(itheta.gt.636)print *,i,j,theta,itheta
			if(itheta.gt.636)itheta=636
			glaremapallpolar(itheta,jr)=glaremapall(i,j)	!*(r/600.)
c
c	Here form the model, subtract...at the end, remove the squares and covered pixels
c			
			arg=6.0+4.0*((theta-2.)/25.)**2
			if(jr.gt.12)arg=arg+1.5
c			if(jr.gt.14)arg=arg+0.5
c			if(jr.gt.16)arg=arg+0.5
			arg1=13.0*(exp(-((theta+27.0-rm/22.)**2)/22.))*exp(-(((rm-22.)/2.)**2)/260.)
			if(jr.le.10)arg1=arg1/1.15
			arg=arg+arg1
			arg1=6.5*(exp(-((theta-25.5)**2)/22.))*exp(-(((rm-27.)/2.)**2)/260.)
			if(jr.le.10)arg1=arg1/1.15
			arg=arg+arg1
			if(theta.lt.-28.0)arg=arg*exp((theta+28.0)/10.)
			if(theta.gt.+28.0)arg=arg*exp((28.0-theta)/10.)
			glaremapallpolar(itheta,jr)=glaremapallpolar(itheta,jr)-arg
			if(i.gt.10.and.i.lt.633)glaremapmodel(i,j)=arg
			if(j.le.50.and.(i.lt.50.or.i.gt.590))glaremapmodel(i,j)=0.
		endif
	enddo
	enddo
c	
	print *,'glaremap formed, beginning output'
	open (11,file = 'd:\sean\eclipse\glare3_288_allpolardif.grd')
	head(3)='-31.8 +31.8'
	head(4)='0 64'
	write(11,1)(head(k),k=1,5)
	write(11,2)((glaremapallpolar(i,j),i=1,636),j=1,128)
	close(11)
c
c	fill holes in glare map
c
	do j=2,127
	do i=12,631
		if(glaremapmodel(i,j).eq.0.)then
			sum=0.
			isum=0
			do ii=i-1,i+1
			do jj=j-1,j+1
				if(glaremapmodel(ii,jj).ne.0.)then
					sum=sum+glaremapmodel(ii,jj)
					isum=isum+1
				endif
			enddo
			enddo
			if(isum.eq.8)then
				glaremapmodel(i,j)=sum/float(isum)
			endif
		endif
	enddo
	enddo
c		
	open (11,file = 'd:\sean\eclipse\glare3_2004_288_CCD_model.grd')
	write(11,1)(headsv(k),k=1,5)
	write(11,2)((glaremapmodel(i,j),i=1,636),j=1,128)
	close(11)
c	
	stop 'done'
	end

