C+
C NAME:
C	orbit_average
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
	implicit none

	integer*4	i, j, k, s, hour, date, numorb, ic, mode
	integer*4	eqorbit(1200,1200)	/1440000*0/
	integer*8	eqimage(1200,1200)	/1440000*0/
	integer*4	eqimagecnt(1200,1200)	/1440000*0/
	integer*4	plorbit(1600,1600)	/2560000*0/
	integer*4	plimage(1600,1600)	/2560000*0/
	integer*4	plimagecnt(1600,1600)	/2560000*0/
	character	eqinfile*27, plinfile*25
	character	dum, outfile*26


	ic=2
	mode=1
	
	do s=1,5
	 print*, s,' Beginning Orbit Averaging Process...'

	 do i=1,1600
	  do j=1,1600
	    plimage(i,j)=0
	    plimagecnt(i,j)=0
	    if(i.le.1200.and.j.le.1200)then
	      eqimage(i,j)=0
	      eqimagecnt(i,j)=0
	    endif
	  enddo
	 enddo
	 
	 numorb = 10

	 if(s.eq.1) then
	   write (outfile,'(A5,I1.1,A2,I1.1,A1,I2.2,A14)') 'EQ1_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	 elseif(s.eq.2) then
	   write (outfile,'(A5,I1.1,A2,I1.1,A1,I2.2,A14)') 'EQ2_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	 elseif(s.eq.3) then
	   write (outfile,'(A5,I1.1,A2,I1.1,A1,I2.2,A14)') 'EQ3_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	 elseif(s.eq.4) then
	   write (outfile,'(A5,I1.1,A2,I1.1,A1,I2.2,A14)') 'SPL_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	 else
	   write (outfile,'(A5,I1.1,A2,I1.1,A1,I2.2,A14)') 'NPL_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	 endif

	 do k=1,numorb

	  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=19
	  elseif(k.eq.12) then
	    date=0529
	    hour=21
	  endif

	  if(s.eq.1)then
	    write (eqinfile,'(A7,I4.4,A1,I2.2,A13)') 'Orbits/',date,'_',hour,'_EqPlot_1.grd'
	  elseif(s.eq.2)then
	    write (eqinfile,'(A7,I4.4,A1,I2.2,A13)') 'Orbits/',date,'_',hour,'_EqPlot_2.grd'
	  elseif(s.eq.3)then
	    write (eqinfile,'(A7,I4.4,A1,I2.2,A13)') 'Orbits/',date,'_',hour,'_EqPlot_3.grd'
	  elseif(s.eq.4)then
	    write (plinfile,'(A7,I4.4,A1,I2.2,A11)') 'Orbits/',date,'_',hour,'_spolar.grd'
	  else
	    write (plinfile,'(A7,I4.4,A1,I2.2,A11)') 'Orbits/',date,'_',hour,'_npolar.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,1200),j=1,1200)
12	    format(20i6)
	    close(12)

	    !! Put orbit into total orbit picture
	    do i=1,1200
	     do j=1,1200
	       if(eqorbit(i,j).gt.0)then
	         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,1600),j=1,1600)
17	    format(20i6)
	    close(17)

	    !! Subtract second orbit from the first
	    do i=1,1600
	     do j=1,1600
	       if(plorbit(i,j).gt.0)then
	         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

	  if(s.lt.4)then

	    !! Put orbit into total orbit picture
	    do i=1,1200
	     do j=1,1200
	       if(eqimagecnt(i,j).gt.0) eqimagecnt(i,j)=nint( float(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'/'1200 1200'/'0 120'/'-60 60'/'0 1400')
	    elseif(s.eq.2)then
	      write(14,15)
15	      format('DSAA'/'1200 1200'/'120 240'/'-60 60'/'0 1400')
	    else
	      write(14,16)
16	      format('DSAA'/'1200 1200'/'240 360'/'-60 60'/'0 1400')
	    endif
	    do j=1,1200
	      write(14,'(20i6)')(eqimagecnt(i,j),i=1,1200)
	    enddo
	    close(14)

	  else
	  
	    !! Subtract second orbit from the first
	    do i=1,1600
	     do j=1,1600
	       if(plimagecnt(i,j).gt.0) plimagecnt(i,j)=nint( float(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'/'1600 1600'/'-40 40'/'-40 40'/'0 1400')
	    do j=1,1600
	      write(19,'(20i6)')(plimagecnt(i,j),i=1,1600)
	    enddo
	    close(19)

	  endif
	
	enddo  !!end s loop to include all 5 maps of the sky

	print*, 'Fortran out!!!'
	end