C+
C NAME:
C	averagenight.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	Jan 2003 Aaron Smith (UCSD/CASS)
C		C++ function converted to present fortran subroutine
C-
	subroutine averagenight(node, nodeid, hits, index, level, cra, cdec, imagecnt, npolarimagecnt, spolarimagecnt)

	integer*8 index, hitcnt
	integer*2 level, hits(index), ibad
	integer*4 imagecnt(1200,1200), npolarimagecnt(1600,1600), spolarimagecnt(1600,1600), xarg, yarg, i, j
	integer*8 nodeid(index), imagenode
	real*4 node(index), image(1200,1200), npolarimage(1600,1600), spolarimage(1600,1600), noderesponse
	real*8 cra, cdec, ra, dec, sdec, rmin, rmax
c
	hitcnt=0
	rmin=65000
	rmax=0

	do i=1,1200
	 do j=1,1200
	   image(i,j)=0.0
	   imagecnt(i,j)=0
	 enddo
	enddo
c
	do i=1,1600
	 do j=1,1600
	   npolarimage(i,j)=0.0
	   npolarimagecnt(i,j)=0
	   spolarimage(i,j)=0.0
	   spolarimagecnt(i,j)=0
	 enddo
	enddo
c
	do i=1,index
	  if(node(i).gt.0.0.and.hits(i).gt.0)then	!Here test on both response and hits
	    if(rmax.lt.node(i))rmax=node(i)
	    if(rmin.gt.node(i))rmin=node(i)
	    hitcnt=hitcnt+1
	    imagenode = nodeid(i)
	    noderesponse = node(i)/float(hits(i))

	    CALL TRIANGLECENTER(level, imagenode, ra, dec)

	    if( (ra-cra).gt.180.D0 ) ra=ra-360.D0
	    if( (ra-cra).lt.-180.D0 ) ra=ra+360.D0
c
c	   Make RA/dec plot
c
	      ibad=0
	      xarg = nint( -(ra-cra)*10.D0 + 600.D0 )
	      yarg = nint( (dec-cdec)*10.D0 + 600.D0 )
	      if(xarg.lt.1.or.xarg.gt.1200.or.yarg.lt.1.or.yarg.gt.1200)ibad=1
	      if(ibad.eq.0)then
	        image(xarg,yarg) = image(xarg,yarg) + noderesponse
	        imagecnt(xarg,yarg) = imagecnt(xarg,yarg) + 1
	      endif
c
c	   Make North Polar plot
c
	    if(dec.gt.33.D0)then
	      ibad=0
	      xarg = nint( 20.*(90.-dec)*cosd(-ra) + 800.)
	      yarg = nint( 20.*(90.-dec)*sind(-ra) + 800.)

	      if(xarg.lt.1.or.xarg.gt.1600.or.yarg.lt.1.or.yarg.gt.1600)ibad=1
	      if(ibad.eq.0)then
	        npolarimage(xarg,yarg) = npolarimage(xarg,yarg) + noderesponse
	        npolarimagecnt(xarg,yarg) = npolarimagecnt(xarg,yarg) + 1
	      endif
	    endif
c
c	   Make South Polar plot
c
	    if(dec.lt.-33.D0)then
	      ibad=0
	      sdec=-dec
	      xarg = nint( 20.*(90.-sdec)*cosd(-ra) + 800.)
	      yarg = nint( 20.*(90.-sdec)*sind(-ra) + 800.)

	      if(xarg.lt.1.or.xarg.gt.1600.or.yarg.lt.1.or.yarg.gt.1600)ibad=1
	      if(ibad.eq.0)then
	        spolarimage(xarg,yarg) = spolarimage(xarg,yarg) + noderesponse
	        spolarimagecnt(xarg,yarg) = spolarimagecnt(xarg,yarg) + 1
	      endif
	    endif

	  endif
	enddo

	do i=1,1200
	 do j=1,1200
	   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
	 enddo
	enddo
	do i=1,1600
	 do j=1,1600
	   if(npolarimagecnt(i,j).gt.0) then
	     npolarimage(i,j)=npolarimage(i,j)/float(npolarimagecnt(i,j))
	     npolarimagecnt(i,j)=nint(npolarimage(i,j))
	   endif
	   if(spolarimagecnt(i,j).gt.0) then
	     spolarimage(i,j)=spolarimage(i,j)/float(spolarimagecnt(i,j))
	     spolarimagecnt(i,j)=nint(spolarimage(i,j))
	   endif
	 enddo
	enddo
	
	print*, rmin,'  < RESPONSE2 < ',rmax
	print*, 'Hit Count: ',hitcnt

	return
	end