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

	character	infile1*50, infile2*50, outfile*17, outfile2*43, dum*20
	integer*2	maxrad
	integer*4	star(200,200)
	integer*4       i, j,ss, date, orbitnum, maxr, minr, x0, y0, ii, jj
	integer*4	imin, imax, jmin, jmax, iarg, jarg, response, test
	integer*8	xsum, ysum, sum
	real*8		thetax, radius, mult
	real*8		xcen, ycen,dra, ddec, xoff, yoff

	print*, 'Fortran in...'

	do ss=1,7
	
	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__cent.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____cent.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___cent.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_____cent.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___cent.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__cent.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___cent.grd'
	endif
	print*, outfile

	do maxrad=21,45,3
C	  maxrad=36
	  write (outfile2,'(A23,A13,A1,I2.2,A4)') 'Star_Transits/centroid/',outfile(1:13),'_',maxrad,'.grd'


	do i=1,200
	 do j=1,200
	   star(i,j)=0
	 enddo
	enddo

	!! 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

	sum=0
	xsum=0
	ysum=0

	!! find centroid of star
	do j=1,200
	 do i=1,200
	   !!multiply xsqr to make that dimension smaller like an elipse
	   radius = sqrt( float((i-x0)*(i-x0))*1.55 + float((j-y0)*(j-y0)) )
	   if(radius.lt.maxrad.and.star(i,j).gt.0)then
	      response=star(i,j)
	      sum = sum + response
	      xsum = xsum + (i*response)
	      ysum = ysum + (j*response)
	   else
	     star(i,j)=0
	   endif
	 enddo
	enddo
	xcen=dfloat(xsum)/sum
	ycen=dfloat(ysum)/sum
	
	star(nint(xcen),nint(ycen))=-100

	print*, maxrad,'X,Y:     ',xcen,ycen
	print*, '          ',nint(xcen),'           ',nint(ycen)

	!!Write subtracted image to file
	open (10,file = outfile2, form = 'formatted')
	write(10,10)
10	format('DSAA'/'200 200'/'0 200'/'0 200'/'-100 65000')
	do j=1,200
	  write(10,'(20i6)')(star(i,j),i=1,200)
	enddo
	close(10)

	enddo !!end maxrad loop to try different radius'
	enddo !!end ss do loop to include all the stars

	print*, 'Fortran out!!!'

	end
