C+
C NAME:
C	avgstar_subtract_brx
C MODIFICATION HISTORY:
C	July, 2003 A. Smith (UCSD/CASS)
C-
	program avgstar_subtract_brx
	implicit none

	character	infile1*50, infile2*50, infile3*16, outfile*17, outfile2*21, outfile3*25, dum*20
	integer*4	star(200,200), avgstar(200,200), favgstarcnt(200,200), fsubstarcnt(200,200)
	integer*4       i, j,ss, date, orbitnum, maxr, minr, x0, y0, ii, jj
	integer*4	imin, imax, jmin, jmax, iarg, jarg
	real*8		startotal, avgstartotal, thetax
	real*8		radius, mult, substar(200,200), fsubstar(200,200), favgstar(200,200)
	real*8		total, tempavgstar(800,800), ftempavgstar(800,800)

	print*, 'Fortran in...'

	!!Read in average star to be subtracted from star transits
	write (infile3,'(A16)') 'avg_star_str.grd'
	open(11,file=infile3,form='formatted',readonly)
	read(11,119)dum
	read(11,119)dum
	read(11,119)dum
	read(11,119)dum
	read(11,119)dum
119	format(a20)
	read(11,11)((avgstar(i,j),i=1,200),j=1,200)
11	format(20i6)
	close(11)

	do ss=1,7

	do i=1,800
	 do j=1,800
	   tempavgstar(i,j)=0.D0
	   ftempavgstar(i,j)=0.D0
	 enddo
	enddo

	do i=1,200
	 do j=1,200
	   star(i,j)=0
	   substar(i,j)=0.D0
	   fsubstar(i,j)=0.D0
	   fsubstarcnt(i,j)=0
	   favgstar(i,j)=0.D0
	   favgstarcnt(i,j)=0
	 enddo
	enddo

	date=0227
	orbitnum=1
	if(ss.eq.1)then
	  write (infile1,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/brx_transits/Aldebar_',date,'_',orbitnum,'_brx.grd'
	  write (infile2,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/str_transits/Aldebar_',date,'_',orbitnum,'_str.grd'
	  thetax=1.23D0
	  write (outfile,'(A17)') 'Aldebar___sub.grd'
	  write (outfile2,'(A21)') 'Aldebar___bin_sub.grd'
	elseif(ss.eq.2)then
	  write (infile1,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/brx_transits/Betel___',date,'_',orbitnum,'_brx.grd'
	  write (infile2,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/str_transits/Betel___',date,'_',orbitnum,'_str.grd'
	  thetax=19.37D0
	  write (outfile,'(A17)') 'Betel_____sub.grd'
	  write (outfile2,'(A21)') 'Betel_____bin_sub.grd'
	elseif(ss.eq.3)then
	  write (infile1,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/brx_transits/Sirius__',date,'_',orbitnum,'_brx.grd'
	  write (infile2,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/str_transits/Sirius__',date,'_',orbitnum,'_str.grd'
	  thetax=26.19D0
	  write (outfile,'(A17)') 'Sirius____sub.grd'
	  write (outfile2,'(A21)') 'Sirius____bin_sub.grd'
	elseif(ss.eq.4)then
	  write (infile1,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/brx_transits/Vega____',date,'_',orbitnum,'_brx.grd'
	  write (infile2,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/str_transits/Vega____',date,'_',orbitnum,'_str.grd'
	  thetax=-21.89D0
	  write (outfile,'(A17)') 'Vega______sub.grd'
	  write (outfile2,'(A21)') 'Vega______bin_sub.grd'
	elseif(ss.eq.5)then
	  write (infile1,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/brx_transits/Alcaid__',date,'_',orbitnum,'_brx.grd'
	  write (infile2,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/str_transits/Alcaid__',date,'_',orbitnum,'_str.grd'
	  thetax=28.74D0
	  write (outfile,'(A17)') 'Alcaid____sub.grd'
	  write (outfile2,'(A21)') 'Alcaid____bin_sub.grd'
	elseif(ss.eq.6)then
	  write (infile1,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/brx_transits/EtaDrac_',date,'_',orbitnum,'_brx.grd'
	  write (infile2,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/str_transits/EtaDrac_',date,'_',orbitnum,'_str.grd'
	  thetax=4.77D0
	  write (outfile,'(A17)') 'EtaDrac___sub.grd'
	  write (outfile2,'(A21)') 'EtaDrac___bin_sub.grd'
	elseif(ss.eq.7)then
	  write (infile1,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/brx_transits/Almach__',date,'_',orbitnum,'_brx.grd'
	  write (infile2,'(A35,I4.4,A1,I2.2,A8)') 'Star_Transits/str_transits/Almach__',date,'_',orbitnum,'_str.grd'
	  thetax=-22.47D0
	  write (outfile,'(A17)') 'Almach____sub.grd'
	  write (outfile2,'(A21)') 'Almach____bin_sub.grd'
	endif

	!! read in headers, then image array
	open(12,file=infile1,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)((star(i,j),i=1,200),j=1,200)
12	format(20i6)
	close(12)
	
	!! find brightest pixel to draw circle around
	maxr=0
	do j=1,200
	 do i=1,200
	   if(star(i,j).gt.maxr)then
	     maxr=float(star(i,j))
	     x0=i
	     y0=j-2
	   endif
	 enddo
	enddo

	!! put avgstar into much finer grid
	do i=1,200
	 do j=1,200
	   imax=i*4
	   jmax=j*4
	   imin=imax-3
	   jmin=jmax-3
	   do ii=imin,imax
	    do jj=jmin,jmax
	      tempavgstar(ii,jj)=float(avgstar(i,j))/16.D0
	    enddo
	   enddo
	 enddo
	enddo

	!! squeeze the fine grid avg star
	do i=1,800
	 do j=1,800
	   ii = nint( (float(i)-400.5)*dcosd(thetax) + 400.5 )
	   jj = j
	   ftempavgstar(ii,jj)=tempavgstar(i,j)
	 enddo
	enddo

	!! rebin the fine grid to the normal size grid
	do i=1,200
	 do j=1,200
	   imax=i*4
	   jmax=j*4
	   imin=imax-3
	   jmin=jmax-3
	   do ii=imin,imax
	    do jj=jmin,jmax
	      favgstar(i,j)=favgstar(i,j)+ftempavgstar(ii,jj)
	    enddo
	   enddo
	 enddo
	enddo

	startotal=0.D0
	avgstartotal=0.D0
	do i=1,200
	 do j=1,200
	   if(favgstarcnt(i,j).gt.0) favgstar(i,j)=favgstar(i,j)/float(favgstarcnt(i,j))
	   radius = sqrt( float((i-x0)*(i-x0))*1.55 + float((j-y0)*(j-y0)) )
	   if(radius.lt.35.and.star(i,j).gt.0)then
	     startotal=startotal+float(star(i,j))
	     avgstartotal=avgstartotal+favgstar(i,j)
	   endif
	 enddo
	enddo

	mult=startotal/avgstartotal

	do i=1,200
	 do j=1,200
	   if(star(i,j).gt.0)  substar(i,j)=float(star(i,j)) - ( favgstar(i,j)*mult )
	 enddo
	enddo

	total=0.D0
	do i=1,200
	 do j=1,200
	   radius = sqrt( float((i-x0)*(i-x0))*1.55 + float((j-y0)*(j-y0)) )
	   if(radius.lt.35.and.star(i,j).gt.0) total=total+substar(i,j)
	 enddo
	enddo

	print*, '(1) Total Residue: ',total
	
cccccccccccccccccccccccccccccc bin-filter the subtracted picture Cccccccccccccccccccccccccccccccc
	do i=1,200
	  do j=1,200
	       imin=i-1
	       if(imin.lt.1) imin=1
	       imax=i+1
	       if(imax.gt.200) imax=200
	       jmin=j-1
	       if(jmin.lt.1) jmin=1
	       jmax=j+1
	       if(jmax.gt.200) jmax=200
	       do ii=imin,imax
	        do jj=jmin,jmax
		  if(star(ii,jj).gt.0) then
		    fsubstar(i,j)=fsubstar(i,j)+substar(ii,jj)
		    fsubstarcnt(i,j)=fsubstarcnt(i,j)+1
		  endif
		enddo
	       enddo
	  enddo
	enddo

	do i=1,200
	 do j=1,200
	   if(fsubstarcnt(i,j).gt.0) fsubstar(i,j)=fsubstar(i,j)/float(fsubstarcnt(i,j))
	 enddo
	enddo
	
	total=0.D0
	do i=1,200
	 do j=1,200
	   radius = sqrt( float((i-x0)*(i-x0))*1.55 + float((j-y0)*(j-y0)) )
	   if(radius.lt.35) total=total+fsubstar(i,j)
	 enddo
	enddo

	print*, '(2) Total Residue: ',total


	!!Write subtracted image to file
	open (10,file = outfile, form = 'formatted')
	write(10,10)
10	format('DSAA'/'200 200'/'0 200'/'0 200'/'-5000 5000')
	do j=1,200
	  write(10,'(20i6)')(nint(substar(i,j)),i=1,200)
	enddo
	close(10)

	!!Write subtracted image to file
	write (outfile3,'(a4,a21)') '999_',outfile2
	open (13,file = outfile3, form = 'formatted')
	write(13,13)
13	format('DSAA'/'200 200'/'0 200'/'0 200'/'0 65000')
	do j=1,200
	  write(13,'(20i6)')(nint(favgstar(i,j)),i=1,200)
	enddo
	close(13)

	enddo !!end ss do loop to include all the stars

	print*, 'Fortran out!!!'

	end
