Program Eq_Subtract
	Implicit None
	
	character	dum, outgrid*48, frame*28, infile*52
	integer*4	year, day, hour, min, sec, year2, day2, hour2, min2, sec2
	integer*4	minorb, maxorb, totalorbnum, i, j, k, ic
	integer*4	instar(200,200), orbloop, orbloopmin, orbloopmax
	integer*4	histogramcnt(1000,100)
	real*4		histogram(1000,100), maxmag

	print*, 'Fortran In...'
	
	ic=2
	
70	format(f5.2)
71	format(i5)	
	print*, 'please enter the maximum magnitude to be checked (usually 4.5)'
	read 70, maxmag	
	print*, 'please enter the minimum local orbit'
	read 71, orbloopmin
	print*, 'please enter the maximum local orbit'
	read 71, orbloopmax
	
	do i=1,1000
	do j=1,100
	  histogram(i,j)=0.0
	  histogramcnt(i,j)=0
	enddo
	enddo
	
	do i=1,200
	do j=1,200
	  instar(i,j)=0
	enddo
	enddo
	
	if(ic.eq.1)then
	  open(10,file='stdstar_c1_fixd.grd',readonly)
	elseif(ic.eq.2)then
	  open(10,file='stdstar_c2_fixd.grd',readonly)
	else
	  open(10,file='stdstar_c3_fixd.grd',readonly)
	endif
	read(10,9)dum
	read(10,9)dum
	read(10,9)dum
	read(10,9)dum
	read(10,9)dum
	do j=1,200
	  read(10,11)(instar(i,j),i=1,200)
	enddo
9	format(a10)
11	format(20i6)
	close(10)
	
	open(31,file='start_times2.txt')
	open(32,file='map_list.txt')
	
	print*, 'Beginning Orbit Loop...'
	
	minorb=55555
	do orbloop=1,orbloopmax
	  read(32,320) year, day, hour, min, sec
320	  format(6x,i4,1x,i3.3,1x,i2.2,i2.2,i2.2)
310       format(i7,4x,i4,4x,i3.3,4x,i2.2,1x,i2.2,1x,i2.2)

	  !!Search through start_times2.txt to find the totalorbit number (in all of smei)
	  if(orbloop.ge.orbloopmin)then	
819	    read(31,310) totalorbnum, year2, day2, hour2, min2, sec2
	    if(year.eq.year2.and.day.eq.day2.and.hour.eq.hour2)then
	      if(orbloop.eq.orbloopmin) minorb=totalorbnum
	      write(frame,820) 'c',ic,'sky_',year,'_',day,'_',hour,min,sec,'.fts.gz'
820	      format(A1,i1.1,A4,i4,A1,i3.3,A1,i2.2,i2.2,i2.2,A7)
	      write(infile,821) '/smeidb/sky1/2003may/c',ic,'/',frame
821	      format(A22,i1.1,A1,A28)
	    else
	      !!Haven't found the corresponding orbit, go back and read the next orbit 
	      goto 819
	    endif
	
	    print*, 'Checking orbit: ',totalorbnum
	    call subtracteqmaps(ic,instar,histogram,histogramcnt,infile,minorb,totalorbnum,maxmag)
	    call subtractpolarmaps(ic,instar,histogram,histogramcnt,infile,minorb,totalorbnum,maxmag)
	  endif
	  	
	enddo !!End loop over all the separate orbits
	
	close(31)
	close(32)
	
	do i=1,1000
	do j=1,100
	  if(histogramcnt(i,j).gt.0)  then
	    histogram(i,j)=histogram(i,j)/histogramcnt(i,j)
	  else
	    histogram(i,j)=9.9
	  endif
	enddo
	enddo
	
	!!Write histogram grid to a file
3001	format(A24,i4.4,A1,i3.3,A5,f4.2,A7)
	write(outgrid,3001) 'Chi_Output/histogram_c2_',minorb,'_',orbloopmax-orbloopmin+1,'orbs_',maxmag,'mag.grd'
	open (33,file = outgrid, form = 'formatted')
	write(33,300)
300	format('DSAA')
	write(33,301) orbloopmax-orbloopmin+1,'100'
301	format(i4.4,1x,A3)
	write(33,302) '0',orbloopmax-orbloopmin+1
302	format(A1,1x,i4.4)
	write(33,303)
303	format('0 100'/'-1 1')
	do j=1,100
	  write(33,'(20f6.3)')(histogram(i,j),i=1,orbloopmax-orbloopmin+1)
	enddo
	close(33)
	
999	print*, 'Fortran Out!!!'

	End