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

	character	infile*17, outfile*16, dum*20
	integer*4	image(200,200), rotimage(200,200), rotimagecnt(200,200), rotimage2(100,100)
	integer*4       xnew, ynew, i,ii, jj, j,ss, x0, y0, maxr,minr,response, test
	integer*8	xsum, ysum, sum
	real*8		xcen, ycen, ra, dec,dra, ddec, cra, cdec, theta
	real*4		radius, x, y
	
	print*, 'Fortran in...'

	xsum=0
	ysum=0
	sum=0
	maxr=0
	minr=60000

	do ss=8,8

	do i=1,200
	do j=1,200
	  if(i.le.100.and.j.le.100) rotimage2(i,j)=0
	  image(i,j)=0
	  rotimage(i,j)=0
	  rotimagecnt(i,j)=0
	enddo
	enddo

	if(ss.eq.1)then
CCC	  Coordinates and Files for Aldebaran
	  cra=68.9802D0
	  cdec=16.5093D0
	  theta=-180.D0
	  theta=-6.75
	  write (infile,'(A17)') 'Aldebar_trans.grd'
	  write (outfile,'(A16)') 'Aldebar__rot.grd'
	elseif(ss.eq.2)then
CCC	  Coordinates and Files for HD 20902
	  cra=51.0807D0
	  cdec=49.8611D0
	  write (infile,'(A17)') 'HD20902_trans.grd'
	  write (outfile,'(A16)') 'HD20902__rot.grd'
	elseif(ss.eq.3)then
CCC	  Coordinates and Files for Sirius (HD 48915)
	  cra=101.28717D0
	  cdec=-16.71611D0
	  theta=180.D0-21.D0
	  write (infile,'(A17)') 'Sirius__trans.grd'
	  write (outfile,'(A16)') 'Sirius___rot.grd'
	elseif(ss.eq.4)then
CCC	  coordinates and files for Vega (alfa lyra)
	  cra=279.23475D0
	  cdec=38.78361D0
	  theta=21.9D0+6.75D0+8.D0
	  write (infile,'(A17)') 'Vega____trans.grd'
	  write (outfile,'(A16)') 'Vega_____rot.grd'
	elseif(ss.eq.5)then
CCC	  coordinates and files for Alcaid
	  cra=206.8852D0
	  cdec=49.3133D0
	  theta=-24.D0
	  write (infile,'(A17)') 'Alcaid__trans.grd'
	  write (outfile,'(A16)') 'Alcaid___rot.grd'
	elseif(ss.eq.6)then
CCC	  coordinates and files for Eta Draconis
	  cra=245.99786D0
	  cdec=61.51421D0
	  write (infile,'(A17)') 'EtaDrac_trans.grd'
	  write (outfile,'(A16)') 'EtaDrac__rot.grd'
	elseif(ss.eq.7)then
CCC	  coordinates of Polaris
	  cra=37.752D0
	  cdec=89.264D0
	  write (infile,'(A17)') 'Polaris_trans.grd'
	  write (outfile,'(A16)') 'Polaris_cent.grd'
	elseif(ss.eq.8)then
CCC	  Coordinates and Files for shape tester
	  theta=-6.75
	  cra=70.2662D0
	  cdec=0.07366D0
	  write (infile,'(A17)') 'rect_00_trans.grd'
	  write (outfile,'(A16)') 'rect_00__rot.grd'
	elseif(ss.eq.9)then
CCC	  Coordinates and Files for shape tester
	  theta=-6.75
	  cra=90.7083D0
	  cdec=1.4879D0
	  write (infile,'(A17)') 'rect_00_tran2.grd'
	  write (outfile,'(A16)') 'rect_00_rot2.grd'
	elseif(ss.eq.10)then
CCC	  Coordinates and Files for shape tester
	  theta=-6.75
	  cra=90.4381D0
	  cdec=38.60867D0
	  write (infile,'(A17)') 'rect_40_trans.grd'
	  write (outfile,'(A16)') 'rect_40__rot.grd'
	endif


	!! read in headers, then image array
	open(12,file=infile,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)((image(i,j),i=1,200),j=1,200)
12	format(20i6)
	close(12)

	!! find brightest pixel to draw box around
	do j=1,200
	 do i=1,200
	   if(image(i,j).gt.maxr)then
	     maxr=image(i,j)
	     x0=i
	     y0=j
	   endif
	 enddo
	enddo
	
	y0=y0-3
	
	!!find background within circle
	do j=1,200
	 do i=1,200
c	   radius = sqrt( float((i-x0)*(i-x0) + (j-y0)*(j-y0)) )
c	   if(radius.lt.18)then
	     if(image(i,j).lt.minr) minr=image(i,j)
c	   endif
	 enddo
	enddo

	!! find centroid of star
	do j=1,200
	 do i=1,200
c	   radius = sqrt( float((i-x0)*(i-x0) + (j-y0)*(j-y0)) )
c	   if(radius.lt.18)then
	      response=image(i,j)-minr
	      sum = sum + response
	      xsum = xsum + (i*response)
	      ysum = ysum + (j*response)
c	   endif
	 enddo
	enddo
	xcen=dfloat(xsum)/dfloat(sum)
	ycen=dfloat(ysum)/dfloat(sum)

	image(nint(xcen),nint(ycen))=-100

	ra = cra - (xcen - 100.5D0)*0.025D0
	dec = cdec + (ycen - 100.5D0)*0.025D0

	ddec=dec-cdec
	dra=ra-cra

	print*, 'X,Y:     ',xcen,ycen
	print*, 'RA/DEC:  ',ra,dec
	print*, 'Discrep: ',dra,ddec

	do i=1,200
	 do j=1,200
c	   rotimage(i,j)=image(i,j)

	   x=float(i)-100.5D0
	   y=float(j)-100.5D0

	   !!find rotated coordinates
	   xnew=nint( x*dcosd(theta)+y*dsind(theta)+100.5D0 )
	   ynew=nint( -x*dsind(theta)+y*dcosd(theta)+100.5D0 )

	   if(xnew.gt.0.and.xnew.lt.201.and.ynew.gt.0.and.ynew.lt.201)then
	     rotimage(xnew,ynew)=rotimage(xnew,ynew)+image(i,j)
	     rotimagecnt(xnew,ynew)=rotimagecnt(xnew,ynew)
	   endif
	 enddo
	enddo

	do i=1,200
	 do j=1,200
	   if(rotimagecnt(i,j).gt.0) rotimage(i,j)=rotimage(i,j)/rotimagecnt(i,j)
	   rotimagecnt(i,j)=0
	 enddo
	enddo
	
c	do i=1,100
c	 do j=1,100
c	  ii=i*2
c	  jj=j*2
c	  rotimage2(i,j)=rotimage(ii,jj)+rotimage(ii-1,jj)+rotimage(ii-1,jj-1)+rotimage(ii,jj-1)
c	 enddo
c	enddo
	
	do i=1,200
	 do j=1,200
	   if(rotimage(i,j).gt.0)then
	      rotimagecnt(i,j)=0
	      ra = cra - (i - 100.5D0)*0.025D0
	      dec = cdec + (j - 100.5D0)*0.025D0

	      xnew = nint( -(ra-cra)*40.D0*dcosd(dec) + 50.5D0 )
	      ynew = nint( (dec-cdec)*40.D0 + 50.5D0 )
	      if(xnew.ge.1.and.xnew.le.100.and.ynew.ge.1.and.ynew.le.100)then
	        rotimage2(xnew,ynew) = rotimage2(xnew,ynew) + rotimage(i,j)
	        rotimagecnt(xnew,ynew) = rotimagecnt(xnew,ynew) + 1
	      endif
	   endif
	 enddo
	enddo
	
	do i=1,100
	 do j=1,100
	   if(rotimagecnt(i,j).gt.0) rotimage(i,j)=rotimage(i,j)/rotimagecnt(i,j)
	 enddo
	enddo

	!!Write rotated image to file
	open (10,file = outfile, form = 'formatted')
	write(10,10)
10	format('DSAA'/'100 100'/'0 100'/'0 100'/'0 10000')
	do j=1,100
	  write(10,'(20i6)')(rotimage2(i,j),i=1,100)
	enddo
	close(10)
	
	enddo !!end ss do loop to include all the stars

	print*, 'Fortran out!!!'

	end
