C+
C NAME:	allskydiff
C PURPOSE:
C	This program will take difference images and combines them together
C CATEGORY:
C	SMEI data analysis
C OUTPUTS:
C	A grid file that represents the combination map all skyt in RA and Dec
C PROCEDURE:
C	Read in the orbits to be used, combine them together
C MODIFICATION HISTORY:
C	September 2003, Sean & Andy (UCSD/CASS)
C-
	PROGRAM allskydiff
	implicit none

	integer*4	i, j, k, s, hour, date, numorbstart,numorbend,numorb, ic, mode
	integer*4	m,mmax,iarg,iiarg,jarg,auroraswitch
	real*4		eqorbit(600,600)
	real*4		eqimage(600,600)
	integer*4	eqimagecnt(600,600)
	real*4		plorbit(400,400)
	real*4		plimage(400,400)
	integer*4	plimagecnt(400,400)
	real*4		allsky(1800,900),radius,x,y,angle
	integer*4	allskycnt(1800,900)
	character	dum, infile*33, outfile*35, outfile1*35
	character	aurinfile*40, auroutfile*42, auroutfile1*42


	ic=2
	mode=4
	numorbstart=1
	numorbend = 18
	numorb=8
	auroraswitch=1
	 do k=numorbstart,numorbend
	 if(auroraswitch.eq.0)print*, 'Camera 2, Orbit ',k,'Combining ...'
	 if(auroraswitch.eq.1)print*, 'Camera 2, Aurora, Orbit ',k,'Combining ...'

	  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

	do i=1,1800
	 do j=1,900
	   allsky(i,j)=0.
	   allskycnt(i,j)=0
	 enddo
	enddo
	if(auroraswitch.eq.0)then
	   write (outfile,'(A19,I4.4,A1,I2.2,A9)') 'Orbits/Allsky_Cam2_',date,'_',hour,'_diff.grd'
	   write (outfile1,'(A19,I4.4,A1,I2.2,A9)') 'Orbits/Allsky_Cam2_',date,'_',hour,'_diff.asd'
	else
	   write (auroutfile,'(A26,I4.4,A1,I2.2,A9)') 'Orbits/Aurora/Allsky_Cam2_',date,'_',hour,'_diff.grd'
	   write (auroutfile1,'(A26,I4.4,A1,I2.2,A9)') 'Orbits/Aurora/Allsky_Cam2_',date,'_',hour,'_diff.asd'
	endif   
	
	do s=1,5

	 if(s.lt.4)then
	   do i=1,600
	    do j=1,600
	      eqimage(i,j)=0.
	      eqimagecnt(i,j)=0
	    enddo
	   enddo
	 else
	   do i=1,400
	    do j=1,400
	      plimage(i,j)=0.
	      plimagecnt(i,j)=0
  	    enddo
  	   enddo
	 endif

	 if(auroraswitch.eq.0)then
	 if(s.eq.1) then
	   write (infile,'(A15,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/EQ1_Sub_',date,'_',hour,'_',numorb,'_avg.grd'
	 elseif(s.eq.2) then
	   write (infile,'(A15,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/EQ2_Sub_',date,'_',hour,'_',numorb,'_avg.grd'
	 elseif(s.eq.3) then
	   write (infile,'(A15,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/EQ3_Sub_',date,'_',hour,'_',numorb,'_avg.grd'
	 elseif(s.eq.4) then
	   write (infile,'(A15,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/SPL_Sub_',date,'_',hour,'_',numorb,'_avg.grd'
	 else
	   write (infile,'(A15,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/NPL_Sub_',date,'_',hour,'_',numorb,'_avg.grd'
	 endif
	 else
	 if(s.eq.1) then
	   write (aurinfile,'(A22,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/Aurora/EQ1_Sub_',date,'_',hour,'_',numorb,'_avg.grd'
	 elseif(s.eq.2) then
	   write (aurinfile,'(A22,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/aurora/EQ2_Sub_',date,'_',hour,'_',numorb,'_avg.grd'
	 elseif(s.eq.3) then
	   write (aurinfile,'(A22,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/Aurora/EQ3_Sub_',date,'_',hour,'_',numorb,'_avg.grd'
	 elseif(s.eq.4) then
	   write (aurinfile,'(A22,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/Aurora/SPL_Sub_',date,'_',hour,'_',numorb,'_avg.grd'
	 else
	   write (aurinfile,'(A22,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/Aurora/NPL_Sub_',date,'_',hour,'_',numorb,'_avg.grd'
	 endif
	 endif

	 if(s.lt.4)then  !!do eq plots: the only reason to convert to eqimage, plimage is to separate out zeroes/missing

	    !! read in headers and orbit image arrays
	    if(auroraswitch.eq.0)open(12,file=infile,form='formatted',readonly)
	    if(auroraswitch.eq.1)open(12,file=aurinfile,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(20f10.2)
	    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).ne.-99999.99)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
	    if(auroraswitch.eq.0)open(17,file=infile,form='formatted',readonly)
	    if(auroraswitch.eq.1)open(17,file=aurinfile,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(20f10.2)
	    close(17)

	    do i=1,400
	     do j=1,400
	       if(plorbit(i,j).ne.-99999.99)then
	         plimage(i,j)=plimage(i,j)+plorbit(i,j)
		 plimagecnt(i,j)=plimagecnt(i,j)+1
	       endif
	     enddo
	    enddo

	 endif

	    !! 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(eqimagecnt(i,j).ne.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(plimagecnt(i,j).ne.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)
c			mmax=0
		 	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
		 	if(m.eq.0)then
		 	 	allsky(iiarg,jarg)=plimage(i,j)
		 		allskycnt(iiarg,jarg)=plimagecnt(i,j)
		 	else
		 		if(allskycnt(iiarg,jarg).eq.0)allskycnt(iiarg,jarg)=1000
		 	endif
		 	enddo
		 endif
	       endif
	     enddo
	    enddo
	 endif
	    
  
	enddo  !!end s loop to include all 5 maps of the sky
	
	if(auroraswitch.eq.0)then
		open (20,file=outfile, form = 'formatted')
		open (21,file=outfile1, form = 'formatted')
	else
		open (20,file=auroutfile, form = 'formatted')
		open (21,file=auroutfile1, form = 'formatted')
	endif
	write(20,20)
20	format('DSAA'/'1800 900'/'0 360'/'-90 90'/'-10 10')
	do j=1,900
	    do i=1,1800
	    if(allskycnt(i,j).eq.0)allsky(i,j)=-99999.99
	    if(allskycnt(i,j).gt.0)allsky(i,j)=allsky(i,j)/float(allskycnt(i,j))
	    if(allskycnt(i,j).ge.1000)allsky(i,j)=+99999.99
	    enddo
	    write(20,'(20f10.2)')(allsky(i,j),i=1,1800)
c	    do i=1,1800
c	    if(allskycnt(i,j).eq.0)allsky(i,j)=-99999.99	! Removed leading 9 here
c	    enddo
	    write(21,'(20f10.2)')(allsky(i,j),i=1,1800)
	enddo
	close(20)
	close(21)
	enddo	!end k-loop	    

	print*, 'Fortran out!!!'
	end