C+
C NAME:
C	averagenight_transit.f
C PURPOSE:
C	This subroutine takes the three triangle arrays(node, nodeid, and hits)
C	and averages them and puts them into a two-dimensional data array(image)
C	which can be printed out and viewed as a picture using Surfer
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call averagenight(node,nodeid,hits,TOP_KMAX,level,cra,cdec,image)
C INPUTS:
C	node		real*4		array of size TOP_KMAX storing the response value
C					corresponding to a triangle of any index
C	nodeid		integer*8	array of size TOP_KMAX storing the name of a triangle
C					in format(uint64 in C++) useful for JHU Spatial Index routines
C	hits		integer*2	array of size TOP_KMAX storing the number of frames contributing
C					to a triangles response(node array)
C	TOP_KMAX	integer*8	the size of the three triangle arrays
C	level		integer*2	the level of the spatial index to be used
C	cra,cdec	real*8		the ra/dec coordinates of the center of the field of view for
C					the image array
C	imagecnt	integer*4	array of size (3200,1600), storing the final image to be viewed
C	polarimagecnt	integer*4	array of size(1600,1600) storing final polar(dec>50) image
C OUTPUTS:
C	outputs the final image to printed out and viewed by the Surfer imaging program
C CALLS:
C	TRIANGLECENTER			C++ function returning the ra/dec of the center of a triangle
C MODIFICATION HISTORY:
C	May 2003 Aaron Smith (UCSD/CASS)
C
C-
	subroutine averagenight_rotate_2(node,nodeid,hits,index,level,cra,cdec,thetarot,thetax,finalimage,finalimage2)
	implicit none

	integer*8 index, badcnt,cnt
	integer*2 level, hits(index), ibad
	integer*4 imagecnt(200,200), finalimage(200,200), imagecnt2(200,200), finalimage2(200,200)
	integer*4 xarg, yarg, xnew, ynew, i, j, xarg2, yarg2
	integer*8 nodeid(index), imagenode
	real*4 node(index), image(200,200), image2(200,200), noderesponse, fxarg, fyarg
	real*8 cra, cdec, rac, ra, dec, sdec, sum,sum2, rsum, dsum, rac, decc, xtri, ytri, hyptri
	real*8 alfa, beta, gamma, theta, xra(2), ydec(2), slope, thetarot, thetax

	integer*4	x0, y0,minr, minr2,response,response2, test
	integer*8	xsum, ysum, xsum2, ysum2
	real*8		xcen,xcen2, ycen2, ycen,dra, ddec, xoff,xoff2, yoff2, yoff, ii, jj
	real*8		rdiscrep, maxr, maxr2
	real*4		radius

	sum=0.D0
	rsum=0.D0
	dsum=0.D0

	do i=1,200
	 do j=1,200
	   image(i,j)=0.0
	   image2(i,j)=0.0
	   imagecnt(i,j)=0
	   imagecnt2(i,j)=0
	   finalimage(i,j)=0.0
	   finalimage(i,j)=0.0
	 enddo
	enddo

	!!Euler angles for rotating centroid to [0,0]
	alfa=cra !!line up with x-axis
	beta=-cdec  !!rotate about new y-axis so centroid now at [0,0]
	gamma=0.D0

	badcnt=0
	cnt=0
	do i=1,index
	  if(node(i).gt.0.0.and.hits(i).gt.0)then	!Here test on both response and hits
	    imagenode = nodeid(i)
	    noderesponse = node(i)/float(hits(i))

	    CALL TRIANGLECENTER(level, imagenode, ra, dec)

	    call rotate8(alfa,beta,gamma,ra,dec)

	    !! rotate making y-axis line up with centroid, then rotate around centroid (y-axis)
	    call rotate8(-90.D0,thetarot,0.D0,ra,dec)
	    rac=90.D0  !!center ra is now centered on 90 degrees (lined up with y-axis)

	    sum=sum+noderesponse
	    rsum=rsum+(ra*noderesponse)
	    dsum=dsum+(dec*noderesponse)

	    if( ra.gt.180.D0 ) ra=ra-360.D0
	    if( ra.lt.-180.D0 ) ra=ra+360.D0

	    !! Make first 200x200 image
	      ibad=0
	      xarg = nint( (-(ra-rac)*40.D0)/dcosd(thetax) + 100.5D0 )
	      yarg = nint( dec*40.D0 + 100.5D0 )
	      if(xarg.lt.1.or.xarg.gt.200.or.yarg.lt.1.or.yarg.gt.200)ibad=1
	      if(ibad.eq.0)then
	        image(xarg,yarg) = image(xarg,yarg) + noderesponse
	        imagecnt(xarg,yarg) = imagecnt(xarg,yarg) + 1
	      endif
	    !! make second image
	      ibad=0
	      xarg2 = nint( -(ra-rac)*40.D0 + 100.5D0 )
	      yarg2 = nint( dec*40.D0 + 100.5D0 )
	      if(xarg2.lt.1.or.xarg2.gt.200.or.yarg2.lt.1.or.yarg2.gt.200)ibad=1
	      if(ibad.eq.0)then
	        image2(xarg2,yarg2) = image2(xarg2,yarg2) + noderesponse
	        imagecnt2(xarg2,yarg2) = imagecnt2(xarg2,yarg2) + 1
	      endif

	  endif
	enddo

	do i=1,200
	 do j=1,200
	   if(imagecnt(i,j).gt.0) then
	     image(i,j)=image(i,j)/float(imagecnt(i,j))
	     imagecnt(i,j)=nint(image(i,j))
	   endif
	   if(imagecnt2(i,j).gt.0) then
	     image2(i,j)=image2(i,j)/float(imagecnt2(i,j))
	     imagecnt2(i,j)=nint(image2(i,j))
	   endif
	 enddo
	enddo

ccccccccccccccccccccccc CENTROIDING/CALIBRATING STAR ccccccccccccccccccccccccccccccccccccccc

	!! find brightest pixel to draw circle around
	!! only do a small rectangle around the center, anything bigger might pick up a cosmic ray as
	!! the brightest pixel and shift the centroiding circle -> centroid gets shifted along with transit image
	maxr=0.D0
	do j=94,106
	 do i=85,115
c	   radius = sqrt( float((i-100)*(i-100))*1.55 + float((j-100)*(j-100)) )
c	   if(radius.lt.)then
	     if(imagecnt(i,j).gt.nint(maxr))then
	       maxr=float(imagecnt(i,j))
	       x0=i
	       y0=j-2
	     endif
c	   endif
	 enddo
	enddo

	sum=0
	xsum=0
	ysum=0
	sum2=0
	xsum2=0
	ysum2=0

	!!find background within region surrounding star, find max within star
	minr=60000
	maxr=0.D0
	minr2=60000
	maxr2=0.D0
	do j=1,200
	 do i=1,200
	   radius = sqrt( float((i-x0)*(i-x0) + (j-y0)*(j-y0)) )
	   if(radius.lt.40)then
	     if(imagecnt(i,j).gt.0.and.imagecnt(i,j).lt.minr) minr=imagecnt(i,j)
	     if(imagecnt(i,j).gt.nint(maxr)) maxr=float(imagecnt(i,j))
	     if(imagecnt2(i,j).gt.0.and.imagecnt2(i,j).lt.minr2) minr2=imagecnt2(i,j)
	     if(imagecnt2(i,j).gt.nint(maxr2)) maxr2=float(imagecnt2(i,j))
	   endif
	 enddo
	enddo

	!!subtract
	do j=1,200
	 do i=1,200
	     if(imagecnt(i,j).gt.0)  imagecnt(i,j)=imagecnt(i,j)-minr
	 enddo
	enddo

	maxr=maxr-minr
	!!calibrate the star so it has the same brightness as the others
	do j=1,200
	 do i=1,200
	     imagecnt(i,j)=nint( float(imagecnt(i,j))*(65000.D0/maxr) )
	 enddo
	enddo

	!! 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.36)then
	      response=imagecnt(i,j)
	      sum = sum + response
	      xsum = xsum + (i*response)
	      ysum = ysum + (j*response)
	      response2=imagecnt2(i,j)
	      sum2 = sum2 + response2
	      xsum2 = xsum2 + (i*response2)
	      ysum2 = ysum2 + (j*response2)
	   endif
	 enddo
	enddo
	xcen=dfloat(xsum)/sum
	ycen=dfloat(ysum)/sum
	xcen2=dfloat(xsum2)/sum2
	ycen2=dfloat(ysum2)/sum2

	dec = (ycen - 100.5D0)*0.025D0
	ra = (xcen - 100.5D0)*0.025D0
	rdiscrep = dsqrt(ra*ra + dec*dec)

	print*, '(1)X,Y:     ',xcen,ycen
	print*, '(1)Discrep: ',ra,dec
	print*, '(1)Radial Discrep: ',rdiscrep
	dec = (ycen2 - 100.5D0)*0.025D0
	ra = (xcen2 - 100.5D0)*0.025D0
C	print*, '(2)X,Y:     ',xcen2,ycen2
C	print*, '(2)Discrep: ',ra,dec

	xoff=100.5D0-xcen
	yoff=100.5D0-ycen
	xoff2=100.5D0-xcen
	yoff2=100.5D0-ycen
	do j=1,200
	 do i=1,200
	   ii=float(i)
	   jj=float(j)
	   xarg=nint(ii+xoff)
	   yarg=nint(jj+yoff)
	   xarg2=nint(ii+xoff2)
	   yarg2=nint(jj+yoff2)
	   if(xarg.ge.1.and.xarg.le.200.and.yarg.ge.1.and.yarg.le.200)  finalimage(xarg,yarg)=imagecnt(i,j)
	   if(xarg2.ge.1.and.xarg2.le.200.and.yarg2.ge.1.and.yarg2.le.200)  finalimage2(xarg2,yarg2)=imagecnt2(i,j)
	 enddo
	enddo

	return
	end