C+
C NAME:
C	orbit_average_bin
C PURPOSE:
C	This program will take composite images of many orbits and average them together
C CATEGORY:
C	SMEI data analysis
C OUTPUTS:
C	A grid file that represents the average between many orbits
C PROCEDURE:
C	Read in the orbits to be used, average them together
C MODIFICATION HISTORY:
C	July 2003, Aaron Smith (UCSD/CASS)
C-
	PROGRAM orbit_average_bin
	implicit none

	integer*4	i, j, k, s, hour, date, numorbstart,numorbend,numorb, ic, mode
	integer*4	m,mmax,iarg,iiarg,jarg
	real*4 		delta,variance,amean,arg
	real*4		eqorbit(600,600)
	real*4		eqimage(600,600),sumsqeqimage(600,600)
	integer*4	eqimagecnt(600,600)
	real*4		plorbit(400,400)
	real*4		plimage(400,400),sumsqplimage(400,400)
	integer*4	plimagecnt(400,400)
	real*4		allsky(1800,900),radius,x,y,angle
	integer*4	allskycnt(1800,900)
	character	eqinfile*31, plinfile*29
	character	dum, outfile*33


	ic=2
	mode=4
	numorbstart=2
	numorbend = 6
	numorb=numorbend-numorbstart+1
	do i=1,1800
	 do j=1,900
	   allsky(i,j)=0.
	   allskycnt(i,j)=0
	 enddo
	enddo
	
	do s=1,5

	 if(s.eq.1) then
	   write (outfile,'(A12,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/EQ1_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	 elseif(s.eq.2) then
	   write (outfile,'(A12,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/EQ2_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	 elseif(s.eq.3) then
	   write (outfile,'(A12,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/EQ3_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	 elseif(s.eq.4) then
	   write (outfile,'(A12,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/SPL_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	 else
	   write (outfile,'(A12,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/NPL_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	 endif
	 if(s.lt.4)then
	   do i=1,600
	    do j=1,600
	      eqimage(i,j)=0.
	      sumsqeqimage(i,j)=0.
	      eqimagecnt(i,j)=0
	    enddo
	   enddo
	 else
	   do i=1,400
	    do j=1,400
	      plimage(i,j)=0.
	      sumsqplimage(i,j)=0.
	      plimagecnt(i,j)=0
  	    enddo
  	   enddo
	 endif


	 do k=numorbstart,numorbend
	 print*,'Camera ',ic,', Map #',s,' Orbit ',k,'Averaging ...'

	  if(k.eq.1) then
	    date=0529
	    hour=03
	  elseif(k.eq.2) then
	    date=0529
	    hour=04
	  elseif(k.eq.3) then
	    date=0529
	    hour=06
	  elseif(k.eq.4) then
	    date=0529
	    hour=08
	  elseif(k.eq.5) then
	    date=0529
	    hour=09
	  elseif(k.eq.6) then
	    date=0529
	    hour=11
	  elseif(k.eq.7) then
	    date=0529
	    hour=13
	  elseif(k.eq.8) then
	    date=0529
	    hour=14
	  elseif(k.eq.9) then
	    date=0529
	    hour=16
	  elseif(k.eq.10) then
	    date=0529
	    hour=18
	  elseif(k.eq.11) then
	    date=0529
	    hour=20
	  elseif(k.eq.12) then
	    date=0529
	    hour=21
	  elseif(k.eq.13) then
	    date=0529
	    hour=23
	  elseif(k.eq.14) then
	    date=0530
	    hour=01
	  elseif(k.eq.15) then
	    date=0530
	    hour=02
	  elseif(k.eq.16) then
	    date=0530
	    hour=04
	  elseif(k.eq.17) then
	    date=0530
	    hour=06
	  elseif(k.eq.18) then
	    date=0530
	    hour=07
	  endif

	  if(s.eq.1)then
	    write (eqinfile,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date,'_',hour,'_EqPlot_1_bin.grd'
	  elseif(s.eq.2)then
	    write (eqinfile,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date,'_',hour,'_EqPlot_2_bin.grd'
	  elseif(s.eq.3)then
	    write (eqinfile,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date,'_',hour,'_EqPlot_3_bin.grd'
	  elseif(s.eq.4)then
	    write (plinfile,'(A7,I4.4,A1,I2.2,A15)') 'Orbits/',date,'_',hour,'_spolar_bin.grd'
	  else
	    write (plinfile,'(A7,I4.4,A1,I2.2,A15)') 'Orbits/',date,'_',hour,'_npolar_bin.grd'
	  endif

	  if(s.lt.4)then  !!do eq plots

	    !! read in headers and orbit image arrays
	    open(12,file=eqinfile,form='formatted',readonly)
	    read(12,120)dum
	    read(12,120)dum
	    read(12,120)dum
	    read(12,120)dum
	    read(12,120)dum
120	    format(a20)
	    read(12,12)((eqorbit(i,j),i=1,600),j=1,600)
12	    format(20f9.1)
	    close(12)

	    !! Put orbit into total orbit picture
	   
	    do i=1,600
	     iarg=i+600*(3-s)
	     do j=1,600
	       if(eqorbit(i,j).gt.0.)then
	         sumsqeqimage(i,j)=sumsqeqimage(i,j)+eqorbit(i,j)**2
	         eqimage(i,j)=eqimage(i,j)+eqorbit(i,j)
		 eqimagecnt(i,j)=eqimagecnt(i,j)+1
	       endif
	     enddo
	    enddo

	  else  !!do polar plots

	    !! read in headers and orbit image arrays
	    open(17,file=plinfile,form='formatted',readonly)
	    read(17,122)dum
	    read(17,122)dum
	    read(17,122)dum
	    read(17,122)dum
	    read(17,122)dum
122	    format(a20)
	    read(17,17)((plorbit(i,j),i=1,400),j=1,400)
17	    format(20f9.1)
	    close(17)

	    do i=1,400
	     do j=1,400
	       if(plorbit(i,j).gt.0.)then
	         sumsqplimage(i,j)=sumsqplimage(i,j)+plorbit(i,j)**2
	         plimage(i,j)=plimage(i,j)+plorbit(i,j)
		 plimagecnt(i,j)=plimagecnt(i,j)+1
	       endif
	     enddo
	    enddo

	   endif
	  enddo !!end k loop to include many orbits

	do k=numorbstart,numorbend
c	 print*, 'Map #',s,' Orbit ',k,'Averaging ...'

	  if(k.eq.1) then
	    date=0529
	    hour=03
	  elseif(k.eq.2) then
	    date=0529
	    hour=04
	  elseif(k.eq.3) then
	    date=0529
	    hour=06
	  elseif(k.eq.4) then
	    date=0529
	    hour=08
	  elseif(k.eq.5) then
	    date=0529
	    hour=09
	  elseif(k.eq.6) then
	    date=0529
	    hour=11
	  elseif(k.eq.7) then
	    date=0529
	    hour=13
	  elseif(k.eq.8) then
	    date=0529
	    hour=14
	  elseif(k.eq.9) then
	    date=0529
	    hour=16
	  elseif(k.eq.10) then
	    date=0529
	    hour=18
	  elseif(k.eq.11) then
	    date=0529
	    hour=20
	  elseif(k.eq.12) then
	    date=0529
	    hour=21
	  elseif(k.eq.13) then
	    date=0529
	    hour=23
	  elseif(k.eq.14) then
	    date=0530
	    hour=01
	  elseif(k.eq.15) then
	    date=0530
	    hour=02
	  elseif(k.eq.16) then
	    date=0530
	    hour=04
	  elseif(k.eq.17) then
	    date=0530
	    hour=06
	  elseif(k.eq.18) then
	    date=0530
	    hour=07
	  endif

	  if(s.eq.1)then
	    write (eqinfile,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date,'_',hour,'_EqPlot_1_bin.grd'
	  elseif(s.eq.2)then
	    write (eqinfile,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date,'_',hour,'_EqPlot_2_bin.grd'
	  elseif(s.eq.3)then
	    write (eqinfile,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date,'_',hour,'_EqPlot_3_bin.grd'
	  elseif(s.eq.4)then
	    write (plinfile,'(A7,I4.4,A1,I2.2,A15)') 'Orbits/',date,'_',hour,'_spolar_bin.grd'
	  else
	    write (plinfile,'(A7,I4.4,A1,I2.2,A15)') 'Orbits/',date,'_',hour,'_npolar_bin.grd'
	  endif

	  if(s.lt.4)then  !!do eq plots

	    !! read in headers and orbit image arrays
	    open(12,file=eqinfile,form='formatted',readonly)
	    read(12,120)dum
	    read(12,120)dum
	    read(12,120)dum
	    read(12,120)dum
	    read(12,120)dum
	    read(12,12)((eqorbit(i,j),i=1,600),j=1,600)
	    close(12)
	    !! Put orbit into total orbit picture
	   
	    do i=1,600
	     iarg=i+600*(3-s)
	     do j=1,600
	       if(eqorbit(i,j).gt.0..and.eqimagecnt(i,j).gt.2)then
	         arg=float(eqimagecnt(i,j))
	         amean=eqimage(i,j)/arg
		 variance=(sumsqeqimage(i,j) - arg*(amean**2))/(arg-1.)
		 delta=eqorbit(i,j)-amean
		 if(delta*delta/(arg-1.).gt.0.5*variance)then
		  eqimagecnt(i,j)=eqimagecnt(i,j)-1
		  eqimage(i,j)=eqimage(i,j)-eqorbit(i,j)
		 endif
	       endif
	     enddo
	    enddo

	  else  !!do polar plots

	    !! read in headers and orbit image arrays
	    open(17,file=plinfile,form='formatted',readonly)
	    read(17,122)dum
	    read(17,122)dum
	    read(17,122)dum
	    read(17,122)dum
	    read(17,122)dum
	    read(17,17)((plorbit(i,j),i=1,400),j=1,400)
	    close(17)

	    do i=1,400
	     do j=1,400
	   	if(plorbit(i,j).gt.0..and.plimagecnt(i,j).gt.2)then
	         arg=float(plimagecnt(i,j))
	         amean=plimage(i,j)/arg
		 variance=(sumsqplimage(i,j) - arg*(amean**2))/(arg-1.)
		 delta=plorbit(i,j)-amean
		 if(delta*delta/(arg-1.).gt.0.5*variance)then
		  plimagecnt(i,j)=plimagecnt(i,j)-1
		  plimage(i,j)=plimage(i,j)-plorbit(i,j)
		 endif
	       	endif
	     enddo
	    enddo

	   endif
	  enddo		!!end 2nd k-loop

	    !! Put orbit into total orbit picture
	  
	  if(s.lt.4)then   
	    do i=1,600
	     iarg=i+600*(3-s)
	     do j=1,600
	       if(eqimage(i,j).gt.0.)then
		 allsky(iarg,j+150)=allsky(iarg,j+150)+eqimage(i,j)
		 allskycnt(iarg,j+150)=allskycnt(iarg,j+150)+eqimagecnt(i,j)
	       endif
	     enddo
	    enddo
	  else  !do polar regions
	    do i=1,400
	     do j=1,400
	       if(plimage(i,j).gt.0.)then
		 x=float(i)-200.5
		 y=float(j)-200.5
		 radius=sqrt(x*x+y*y)
		 if(radius.le.150.5)then
		 	jarg=nint(radius)
		 	mmax=nint(225./radius)
		 	if(s.eq.5)jarg=901-jarg
		 	angle=atan2d(y,x)
		 	if(s.eq.4)iarg=1+nint(5.*(180.-angle))
		 	if(s.eq.5)then
		 		if(angle.lt.0.)angle=angle+360.
		 		iarg=1+nint(5.*angle)
		 	endif
		 	do m=-mmax,mmax
		 	iiarg=iarg+m
		 	if(iiarg.le.0)iiarg=1800+iiarg
		 	if(iiarg.gt.1800)iiarg=iiarg-1800
		 	allsky(iiarg,jarg)=allsky(iiarg,jarg)+plimage(i,j)
		 	allskycnt(iiarg,jarg)=allskycnt(iiarg,jarg)+plimagecnt(i,j)
		 	enddo
		 endif
	       endif
	     enddo
	    enddo
	  endif
	    
	  if(s.lt.4)then

	    !! Put orbit into total orbit picture
	    do i=1,600
	     do j=1,600
	       if(eqimagecnt(i,j).gt.0) eqimage(i,j)=eqimage(i,j)/float(eqimagecnt(i,j))
	     enddo
	    enddo

	    !!Write image(orbit difference) to file
	    open (14,file = outfile, form = 'formatted')
	    if(s.eq.1)then
	      write(14,14)
14	      format('DSAA'/'600 600'/'0 120'/'-60 60'/'0 1400')
	    elseif(s.eq.2)then
	      write(14,15)
15	      format('DSAA'/'600 600'/'120 240'/'-60 60'/'0 1400')
	    else
	      write(14,16)
16	      format('DSAA'/'600 600'/'240 360'/'-60 60'/'0 1400')
	    endif
	    do j=1,600
	      write(14,'(20f9.1)')(eqimage(i,j),i=1,600)
	    enddo
	    close(14)

	  else
	  
	    do i=1,400
	     do j=1,400
		if(plimagecnt(i,j).gt.0) plimage(i,j)=plimage(i,j)/float(plimagecnt(i,j))
	     enddo
	    enddo

	    !!Write image(orbit difference) to file
	    open (19,file = outfile, form = 'formatted')
	    write(19,19)
19	    format('DSAA'/'400 400'/'-40 40'/'-40 40'/'0 1400')
	    do j=1,400
	      write(19,'(20f9.1)')(plimage(i,j),i=1,400)
	    enddo
	    close(19)

	  endif
	
	enddo  !!end s loop to include all 5 maps of the sky
	
	open (20,file='g:\sean\orbits\allsky.grd', form = 'formatted')
	write(20,20)
20	format('DSAA'/'1800 900'/'0 360'/'-90 90'/'0 1400')
	do j=1,900
	    do i=1,1800
	    if(allskycnt(i,j).gt.0) allsky(i,j)=allsky(i,j)/float(allskycnt(i,j))
	    enddo
	    write(20,'(20f9.1)')(allsky(i,j),i=1,1800)
	enddo
	close(20)	    

	print*, 'Fortran out!!!'
	end