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, ibin
	integer*4	imin, imax, jmin, jmax, ii, jj,auroraswitch
	real*4		eqorbit1(600,600)
	real*4		eqorbit2(600,600)
	real*4		eqimage(600,600)
	real*4		feqimage(600,600)
	integer*4	feqimagecnt(600,600)
	real*4		plorbit1(400,400)
	real*4		plorbit2(400,400)
	real*4		plimage(400,400)
	real*4		fplimage(400,400)
	integer*4	fplimagecnt(400,400)
	character	eqinfile1*31, eqinfile2*31, eqinfile3*33,plinfile1*29, plinfile2*29, plinfile3*33
	character	dum, outfile*34, outfile3*33
	character	aureqinfile1*38, aureqinfile2*38, aureqinfile3*40,aurplinfile1*36
	character	auroutfile*41, auroutfile3*40, aurplinfile2*36, aurplinfile3*40
c
C      when file3 = 0, the given orbit has the above "date2, hour2" orbit subtracted
c      when file3 = 1, the given orbit has the "orbit_avg.grd" subtracted
c	smooths by +-ibin
c
	file3=1
	ic=2
	mode=4
	numorb=8
	ibin=0
	auroraswitch = 1
	Print *, 'Camera ',ic,', ',numorb,' orbit-average removed, smoothed by +-',ibin

	do z=1,18  !!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=20
	elseif(z.eq.12) then
	    date1=0529
	    hour1=21
	elseif(z.eq.13) then
	    date1=0529
	    hour1=23
	elseif(z.eq.14) then
	    date1=0530
	    hour1=01
	elseif(z.eq.15) then
	    date1=0530
	    hour1=02
	elseif(z.eq.16) then
	    date1=0530
	    hour1=04
	elseif(z.eq.17) then
	    date1=0530
	    hour1=06
	elseif(z.eq.18) then
	    date1=0530
	    hour1=07
	endif

	date2=0529
	hour2=16

	do k=1,5
	  print*, z,k,' Beginning Orbit Subtraction Process...'
	  if(auroraswitch.eq.0)then
	  if(k.eq.1)then
	    write (eqinfile1,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date1,'_',hour1,'_EqPlot_1_bin.grd'
	    write (eqinfile2,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date2,'_',hour2,'_EqPlot_1_bin.grd'
	    write (eqinfile3,'(A12,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/EQ1_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (outfile,'(A15,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'Orbits/EQ1_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (outfile3,'(A15,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/EQ1_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  elseif(k.eq.2)then
	    write (eqinfile1,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date1,'_',hour1,'_EqPlot_2_bin.grd'
	    write (eqinfile2,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date2,'_',hour2,'_EqPlot_2_bin.grd'
	    write (eqinfile3,'(A12,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/EQ2_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (outfile,'(A15,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'Orbits/EQ2_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (outfile3,'(A15,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/EQ2_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  elseif(k.eq.3)then
	    write (eqinfile1,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date1,'_',hour1,'_EqPlot_3_bin.grd'
	    write (eqinfile2,'(A7,I4.4,A1,I2.2,A17)') 'Orbits/',date2,'_',hour2,'_EqPlot_3_bin.grd'
	    write (eqinfile3,'(A12,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/EQ3_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (outfile,'(A15,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'Orbits/EQ3_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (outfile3,'(A15,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/EQ3_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  elseif(k.eq.4)then
	    write (plinfile1,'(A7,I4.4,A1,I2.2,A15)') 'Orbits/',date1,'_',hour1,'_spolar_bin.grd'
	    write (plinfile2,'(A7,I4.4,A1,I2.2,A15)') 'Orbits/',date2,'_',hour2,'_spolar_bin.grd'
	    write (plinfile3,'(A12,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/SPL_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (outfile,'(A15,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'Orbits/SPL_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (outfile3,'(A15,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/SPL_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  else
	    write (plinfile1,'(A7,I4.4,A1,I2.2,A15)') 'Orbits/',date1,'_',hour1,'_npolar_bin.grd'
	    write (plinfile2,'(A7,I4.4,A1,I2.2,A15)') 'Orbits/',date2,'_',hour2,'_npolar_bin.grd'
	    write (plinfile3,'(A12,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/NPL_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (outfile,'(A15,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'Orbits/NPL_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (outfile3,'(A15,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/NPL_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  endif
	  else
	  if(k.eq.1)then
	    write (aureqinfile1,'(A14,I4.4,A1,I2.2,A17)') 'Orbits/Aurora/',date1,'_',hour1,'_EqPlot_1_bin.grd'
	    write (aureqinfile2,'(A14,I4.4,A1,I2.2,A17)') 'Orbits/Aurora/',date2,'_',hour2,'_EqPlot_1_bin.grd'
	    write (aureqinfile3,'(A19,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/Aurora/EQ1_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (auroutfile,'(A22,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'Orbits/Aurora/EQ1_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (auroutfile3,'(A22,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/Aurora/EQ1_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  elseif(k.eq.2)then
	    write (aureqinfile1,'(A14,I4.4,A1,I2.2,A17)') 'Orbits/Aurora/',date1,'_',hour1,'_EqPlot_2_bin.grd'
	    write (aureqinfile2,'(A14,I4.4,A1,I2.2,A17)') 'Orbits/Aurora/',date2,'_',hour2,'_EqPlot_2_bin.grd'
	    write (aureqinfile3,'(A19,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/Aurora/EQ2_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (auroutfile,'(A22,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'Orbits/Aurora/EQ2_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (auroutfile3,'(A22,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/Aurora/EQ2_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  elseif(k.eq.3)then
	    write (aureqinfile1,'(A14,I4.4,A1,I2.2,A17)') 'Orbits/Aurora/',date1,'_',hour1,'_EqPlot_3_bin.grd'
	    write (aureqinfile2,'(A14,I4.4,A1,I2.2,A17)') 'Orbits/Aurora/',date2,'_',hour2,'_EqPlot_3_bin.grd'
	    write (aureqinfile3,'(A19,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/Aurora/EQ3_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (auroutfile,'(A22,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'Orbits/Aurora/EQ3_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (auroutfile3,'(A22,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/Aurora/EQ3_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  elseif(k.eq.4)then
	    write (aurplinfile1,'(A14,I4.4,A1,I2.2,A15)') 'Orbits/Aurora/',date1,'_',hour1,'_spolar_bin.grd'
	    write (aurplinfile2,'(A14,I4.4,A1,I2.2,A15)') 'Orbits/Aurora/',date2,'_',hour2,'_spolar_bin.grd'
	    write (aurplinfile3,'(A19,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/Aurora/SPL_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (auroutfile,'(A22,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'Orbits/Aurora/SPL_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (auroutfile3,'(A22,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/Aurora/SPL_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  else
	    write (aurplinfile1,'(A14,I4.4,A1,I2.2,A15)') 'Orbits/Aurora/',date1,'_',hour1,'_npolar_bin.grd'
	    write (aurplinfile2,'(A14,I4.4,A1,I2.2,A15)') 'Orbits/Aurora/',date2,'_',hour2,'_npolar_bin.grd'
	    write (aurplinfile3,'(A19,I1.1,A2,I1.1,A1,I2.2,A14)') 'Orbits/Aurora/NPL_c',ic,'_m',mode,'_',numorb,'_orbit_avg.grd'
	    write (auroutfile,'(A22,I4.4,A1,I2.2,A1,I4.4,A1,I2.2,A4)') 'Orbits/Aurora/NPL_Sub_',date1,'_',hour1,'_',date2,'_',hour2,'.grd'
	    write (auroutfile3,'(A22,I4.4,A1,I2.2,A1,I2.2,A8)') 'Orbits/Aurora/NPL_Sub_',date1,'_',hour1,'_',numorb,'_avg.grd'
	  endif
	  endif

	  if(k.lt.4)then

	    !! read in headers and orbit image arrays
	    if(auroraswitch.eq.0)open(12,file=eqinfile1,form='formatted',readonly)
	    if(auroraswitch.eq.1)open(12,file=aureqinfile1,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(20f9.1)
	    close(12)

	    if(file3.eq.0)then
	      if(auroraswitch.eq.0)open(13,file=eqinfile2,form='formatted',readonly)
	      if(auroraswitch.eq.1)open(13,file=aureqinfile2,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(20f9.1)
	      close(13)
	    else
	      if(auroraswitch.eq.0)open(25,file=eqinfile3,form='formatted',readonly)
	      if(auroraswitch.eq.1)open(25,file=aureqinfile3,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(20f9.1)
	      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-ibin
	       if(imin.lt.1) imin=1
	       imax=i+ibin
	       if(imax.gt.600) imax=600
	       jmin=j-ibin
	       if(jmin.lt.1) jmin=1
	       jmax=j+ibin
	       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)+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).eq.0) feqimage(i,j)=-99999.99
	       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
	      if(auroraswitch.eq.0)then
	        open (14,file = outfile, form = 'formatted')
	      else
	        open (14,file = auroutfile, form = 'formatted')
	      endif
	      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,'(20f10.2)')(feqimage(i,j),i=1,600)
	      enddo
	      close(14)
	    else
	      !!Write image(orbit difference) to file
	      if(auroraswitch.eq.0)then
	      open (17,file = outfile3, form = 'formatted')
	      else
	      open (17,file = auroutfile3, form = 'formatted')
	      endif
	      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,'(20f10.2)')(feqimage(i,j),i=1,600)
	      enddo
	      close(17)
	    endif

	  else
	  
	    !! read in headers and orbit image arrays
	    if(auroraswitch.eq.0)open(20,file=plinfile1,form='formatted',readonly)
	    if(auroraswitch.eq.1)open(20,file=aurplinfile1,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(20f9.1)
	    close(20)

	    if(file3.eq.0)then
	      if(auroraswitch.eq.0)open(21,file=plinfile2,form='formatted',readonly)
	      if(auroraswitch.eq.1)open(21,file=aurplinfile2,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(20f9.1)
	      close(21)
	    else
	      if(auroraswitch.eq.0)open(22,file=plinfile3,form='formatted',readonly)
	      if(auroraswitch.eq.1)open(22,file=aurplinfile3,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(20f9.1)
	      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-ibin
	       if(imin.lt.1) imin=1
	       imax=i+ibin
	       if(imax.gt.400) imax=400
	       jmin=j-ibin
	       if(jmin.lt.1) jmin=1
	       jmax=j+ibin
	       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)+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).eq.0) fplimage(i,j)=-99999.99
	       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
	      if(auroraswitch.eq.0)then
	      open (23,file = outfile, form = 'formatted')
	      else
	      open (23,file = auroutfile, form = 'formatted')
	      endif
	      write(23,23)
23	      format('DSAA'/'400 400'/'-40 40'/'-40 40'/'-10 10')
	      do j=1,400
	        write(23,'(20f10.2)')(fplimage(i,j),i=1,400)
	      enddo
	      close(23)
	    else
	      !!Write image(orbit difference) to file
	      if(auroraswitch.eq.0)then
	      open (24,file = outfile3, form = 'formatted')
	      else
	      open (24,file = auroutfile3, form = 'formatted')
	      endif
	      write(24,24)
24	      format('DSAA'/'400 400'/'-40 40'/'-40 40'/'-10 10')
	      do j=1,400
	        write(24,'(20f10.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