C+
C NAME:
C	Orbit_Subtract
C PURPOSE:
C	This program will take composite images of two orbits and subtract them
C	from each other
C CATEGORY:
C	SMEI data analysis
C OUTPUTS:
C	A grid file that represents the difference between two orbits
C PROCEDURE:
C	Read in the orbits to be used, subtract them from each other
C MODIFICATION HISTORY:
C	June 2003, Aaron Smith (UCSD/CASS)
C-
	PROGRAM Orbit_Subtract
	implicit none

	integer*4	i, j, k, z, date1, date2, hour1, hour2, file3, ic, mode, numorb
	integer*4	imin, imax, jmin, jmax, ii, jj
	integer*4	eqorbit1(600,600)	/360000*0/
	integer*4	eqorbit2(600,600)	/360000*0/
	integer*4	eqimage(600,600)	/360000*0/
	real*4		feqimage(600,600)	/360000*0/
	integer*4	feqimagecnt(600,600)	/360000*0/
	integer*4	plorbit1(400,400)	/160000*0/
	integer*4	plorbit2(400,400)	/160000*0/
	integer*4	plimage(400,400)	/160000*0/
	real*4		fplimage(400,400)	/160000*0/
	integer*4	fplimagecnt(400,400)	/160000*0/
	character	eqinfile1*24, eqinfile2*24, eqinfile3*26,plinfile1*22, plinfile2*22, plinfile3*26
	character	dum, outfile*27, outfile3*26


	do z=1,10  !!do many orbits

	if(z.eq.1) then
	    date1=0529
	    hour1=03
	elseif(z.eq.2) then
	    date1=0529
	    hour1=04
	elseif(z.eq.3) then
	    date1=0529
	    hour1=06
	elseif(z.eq.4) then
	    date1=0529
	    hour1=08
	elseif(z.eq.5) then
	    date1=0529
	    hour1=09
	elseif(z.eq.6) then
	    date1=0529
	    hour1=11
	elseif(z.eq.7) then
	    date1=0529
	    hour1=13
	elseif(z.eq.8) then
	    date1=0529
	    hour1=14
	elseif(z.eq.9) then
	    date1=0529
	    hour1=16
	elseif(z.eq.10) then
	    date1=0529
	    hour1=18
	elseif(z.eq.11) then
	    date1=0529
	    hour1=19
	elseif(z.eq.12) then
	    date1=0529
	    hour1=21
	endif

	date2=0529
	hour2=18

C      1 if infile3 and outfile3 to be used instead of infile2 and outfile (orbit average used, not consecutive orbit)
	file3=1
	ic=2
	mode=1
	numorb=10

	do k=1,5
	  print*, k,' Beginning Orbit Subtraction Process...'

	  if(k.eq.1)then
	    write (eqinfile1,'(I4.4,A1,I2.2,A17)') date1,'_',hour1,'_EqPlot_1_bin.grd'
	    write (eqinfile2,'(I4.4,A1,I2.2,A17)') date2,'_',hour2,'_EqPlot_1_bin.grd'
	    write (eqinfile3,'(A5,I1.1,A2,I1.1,A1,I2.2,A14)') 'EQ1_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (outfile,'(A8,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'EQ1_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (outfile3,'(A8,I4.4,A1,I2.2,A1,I2.2,A8)') 'EQ1_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  elseif(k.eq.2)then
	    write (eqinfile1,'(I4.4,A1,I2.2,A17)') date1,'_',hour1,'_EqPlot_2_bin.grd'
	    write (eqinfile2,'(I4.4,A1,I2.2,A17)') date2,'_',hour2,'_EqPlot_2_bin.grd'
	    write (eqinfile3,'(A5,I1.1,A2,I1.1,A1,I2.2,A14)') 'EQ2_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (outfile,'(A8,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'EQ2_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (outfile3,'(A8,I4.4,A1,I2.2,A1,I2.2,A8)') 'EQ2_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  elseif(k.eq.3)then
	    write (eqinfile1,'(I4.4,A1,I2.2,A17)') date1,'_',hour1,'_EqPlot_3_bin.grd'
	    write (eqinfile2,'(I4.4,A1,I2.2,A17)') date2,'_',hour2,'_EqPlot_3_bin.grd'
	    write (eqinfile3,'(A5,I1.1,A2,I1.1,A1,I2.2,A14)') 'EQ3_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (outfile,'(A8,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'EQ3_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (outfile3,'(A8,I4.4,A1,I2.2,A1,I2.2,A8)') 'EQ3_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  elseif(k.eq.4)then
	    write (plinfile1,'(I4.4,A1,I2.2,A15)') date1,'_',hour1,'_spolar_bin.grd'
	    write (plinfile2,'(I4.4,A1,I2.2,A15)') date2,'_',hour2,'_spolar_bin.grd'
	    write (plinfile3,'(A5,I1.1,A2,I1.1,A1,I2.2,A14)') 'SPL_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (outfile,'(A8,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'SPL_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (outfile3,'(A8,I4.4,A1,I2.2,A1,I2.2,A8)') 'SPL_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  else
	    write (plinfile1,'(I4.4,A1,I2.2,A15)') date1,'_',hour1,'_npolar_bin.grd'
	    write (plinfile2,'(I4.4,A1,I2.2,A15)') date2,'_',hour2,'_npolar_bin.grd'
	    write (plinfile3,'(A5,I1.1,A2,I1.1,A1,I2.2,A14)') 'NPL_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (outfile,'(A8,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'NPL_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (outfile3,'(A8,I4.4,A1,I2.2,A1,I2.2,A8)') 'NPL_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  endif


	  if(k.lt.4)then

	    !! read in headers and orbit image arrays
	    open(12,file=eqinfile1,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)((eqorbit1(i,j),i=1,600),j=1,600)
12	    format(20i6)
	    close(12)

	    if(file3.eq.0)then
	      open(13,file=eqinfile2,form='formatted',readonly)
	      read(13,121)dum
	      read(13,121)dum
	      read(13,121)dum
	      read(13,121)dum
	      read(13,121)dum
121	      format(a20)
	      read(13,13)((eqorbit2(i,j),i=1,600),j=1,600)
13	      format(20i6)
	      close(13)
	    else
	      open(25,file=eqinfile3,form='formatted',readonly)
	      read(25,125)dum
	      read(25,125)dum
	      read(25,125)dum
	      read(25,125)dum
	      read(25,125)dum
125	      format(a20)
	      read(25,25)((eqorbit2(i,j),i=1,600),j=1,600)
25	      format(20i6)
	      close(25)
	    endif

	    !! Subtract second orbit from the first
	    do i=1,600
	     do j=1,600
	       eqimage(i,j)=0
	       feqimage(i,j)=0
	       feqimagecnt(i,j)=0
	       if(eqorbit1(i,j).gt.0.and.eqorbit2(i,j).gt.0) eqimage(i,j)=eqorbit1(i,j)-eqorbit2(i,j)
	     enddo
	    enddo
	    
	    do i=1,600
	     do j=1,600
	       imin=i-2
	       if(imin.lt.1) imin=1
	       imax=i+2
	       if(imax.gt.600) imax=600
	       jmin=j-2
	       if(jmin.lt.1) jmin=1
	       jmax=j+2
	       if(jmax.gt.600) jmax=600
	       do ii=imin,imax
	        do jj=jmin,jmax
	          if(eqorbit1(ii,jj).gt.0.and.eqorbit2(ii,jj).gt.0)then
		    feqimage(i,j)=feqimage(i,j)+float(eqimage(ii,jj))
		    feqimagecnt(i,j)=feqimagecnt(i,j)+1
	          endif
	        enddo
	       enddo
	     enddo
	    enddo

	    do i=1,600
	     do j=1,600
	       if(feqimagecnt(i,j).gt.0) feqimage(i,j)=feqimage(i,j)/float(feqimagecnt(i,j))
	     enddo
	    enddo

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

	  else
	  
	    !! read in headers and orbit image arrays
	    open(20,file=plinfile1,form='formatted',readonly)
	    read(20,122)dum
	    read(20,122)dum
	    read(20,122)dum
	    read(20,122)dum
	    read(20,122)dum
122	    format(a20)
	    read(20,20)((plorbit1(i,j),i=1,400),j=1,400)
20	    format(20i6)
	    close(20)

	    if(file3.eq.0)then
	      open(21,file=plinfile2,form='formatted',readonly)
	      read(21,123)dum
	      read(21,123)dum
	      read(21,123)dum
	      read(21,123)dum
	      read(21,123)dum
123	      format(a20)
	      read(21,21)((plorbit2(i,j),i=1,400),j=1,400)
21	      format(20i6)
	      close(21)
	    else
	      open(22,file=plinfile3,form='formatted',readonly)
	      read(22,126)dum
	      read(22,126)dum
	      read(22,126)dum
	      read(22,126)dum
	      read(22,126)dum
126	      format(a20)
	      read(22,22)((plorbit2(i,j),i=1,400),j=1,400)
22	      format(20i6)
	      close(22)
	    endif

	    !! Subtract second orbit from the first
	    do i=1,400
	     do j=1,400
	       plimage(i,j)=0
	       fplimage(i,j)=0
	       fplimagecnt(i,j)=0
	       if(plorbit1(i,j).gt.0.and.plorbit2(i,j).gt.0) plimage(i,j)=plorbit1(i,j)-plorbit2(i,j)
	     enddo
	    enddo
	    
	    do i=1,400
	     do j=1,400
	       imin=i-2
	       if(imin.lt.1) imin=1
	       imax=i+2
	       if(imax.gt.400) imax=400
	       jmin=j-2
	       if(jmin.lt.1) jmin=1
	       jmax=j+2
	       if(jmax.gt.400) jmax=400
	       do ii=imin,imax
	        do jj=jmin,jmax
		  if(plorbit1(ii,jj).gt.0.and.plorbit2(ii,jj).gt.0) then
		    fplimage(i,j)=fplimage(i,j)+float(plimage(ii,jj))
		    fplimagecnt(i,j)=fplimagecnt(i,j)+1
		  endif
	        enddo
	       enddo
	     enddo
	    enddo
	    
	    do i=1,400
	     do j=1,400
	       if(fplimagecnt(i,j).gt.0) fplimage(i,j)=fplimage(i,j)/float(fplimagecnt(i,j))
	     enddo
	    enddo

	    if(file3.eq.0)then
	      !!Write image(orbit difference) to file
	      open (23,file = outfile, form = 'formatted')
	      write(23,23)
23	      format('DSAA'/'400 400'/'-40 40'/'-40 40'/'-10 10')
	      do j=1,400
	        write(23,'(20f9.2)')(fplimage(i,j),i=1,400)
	      enddo
	      close(23)
	    else
	      !!Write image(orbit difference) to file
	      open (24,file = outfile3, form = 'formatted')
	      write(24,24)
24	      format('DSAA'/'400 400'/'-40 40'/'-40 40'/'-10 10')
	      do j=1,400
	        write(24,'(20f9.2)')(fplimage(i,j),i=1,400)
	      enddo
	      close(24)
	    endif


	  endif
	enddo  !!end k loop to do all five plots for each orbit
	enddo  !!end z loop to do subtraction for many orbits

	print*, 'Fortran out!!!'
	end