C+
C NAME:
C	indexframe
C PURPOSE:
C	To index an entire frame of TMO data using an NxN binning scheme
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,night,iframe,frame,affine,euler)
C INPUTS:
C	mode			camera binning style; mode=N -> NxN binning
C	ic			camera number(0-4); 0 -> TMO camera
C	id			data frame max x-value (along rows)
C	jd			data frame max y-value (down columns)
C	level			spatial index level to be used
C	TOP_KMAX		dimension of the three triangle arrays(node,nodeid, hits)
C	node(TOP_KMAX)		array storing summed response for a given spatial triangle
C	nodeid(TOP_KMAX)	array storing triangle name in spatial index format(uint64 in C++)
C	hits(TOP_KMAX)		array storing amount of frames contibuting to a triangle's response
C	night			the night of the frame being indexed
C	iframe			the number of the frame being indexed
C	frame(1280,600)		2-D array storing the pixelated ccd picture
C	affine(6)		array storing parameters for affine transformation
C	euler(3)		array stroing three euler angles(alfa,beta,gamma)
C OUTPUTS:
C	node, nodeid, and hits are the only function arguments that this subroutine changes
C	and returns top the main program
C CALLS:
C	rotate8(alfa, beta, gamma, phi,rlat)
C	airplane(night, iframe, i)
C	INDEXPIXEL(level,node,nodeid,hits,response,ra,dec)
C PROCEDURE:
C	find the 4 vertices of a pixel, transform those ccd coordinates to sky coordinates,
C	and send the 4 ra/dec vertices to the C++ INDEXPIXEL routine.
C MODIFICATION HISTORY:
C	Jan, 2003 Aaron Smith (UCSD)
C-
	subroutine indexframe(mode,ic, id, jd,level,TOP_KMAX,node,nodeid,hits,night,iframe,frame,affine,euler)
	implicit none

	integer*8 TOP_KMAX, kk
	real*4 frame(id,jd), yoffset(5), xoffset(5)
	real*4 node(TOP_KMAX)
	integer*8 nodeid(TOP_KMAX)
	integer*2 hits(TOP_KMAX)
	real*8 x(4), y(4), ra(4), dec(4), euler(3), alfa, beta, gamma, affine(6)
	real*8 rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq,czsum
	real*8 x0(5), y0(5),x00(5), y00(5), r0, tempx, phi, rlat, rccd, arg, arg1, response, dusk, i, j
	integer*4 airplane, iairplane, night, iframe, ii, jj
	integer*4 itemp, jtemp, itempmin, itempmax, jtempmin, jtempmax, rcount
	integer*4 xmin, xmax, ymin, ymax, mode, ic, xoff, id, jd
	integer*2 k, level

	data yoffset /0.,65.,59.,62.,55./
	data xoffset /0.,1.,1.,1.,1./
	data x00 /630.D0,634.D0,628.7D0,637.7D0,633.D0/
	data y00 /1291.38D0,1306.1D0,1297.1D0,1304.9D0,1294.3D0/
	r0=1194.5D0
	x0(ic+1)=x00(ic+1)-xoffset(ic+1)
	y0(ic+1)=y00(ic+1)-yoffset(ic+1)
	alfa=euler(1)
	beta=euler(2)
	gamma=euler(3)

	if(mode.eq.2.and.ic.gt.0)then
	  ymin=1
	  ymax=jd
	  xmin=11
	  xmax=632
	  xoff=0
	else if(mode.eq.4.and.ic.gt.0)then
	  ymin=1
	  ymax=jd
	  xmin=6
	  xmax=316
	  xoff=0
	else if(mode.eq.1.and.ic.eq.0)then
	  ymin=1
	  ymax=350
	  xmin=1
	  xmax=1242
	  xoff=22
	else if(mode.eq.2.and.ic.eq.0)then
	  ymin=1
	  ymax=176
	  xmin=1
	  xmax=621
	  xoff=22
	else if(mode.eq.4.and.ic.eq.0)then
	  ymin=1
	  ymax=88
	  xmin=1
	  xmax=310
	  xoff=22
	else if(mode.eq.8.and.ic.eq.0)then
	  ymin=1
	  ymax=44
	  xmin=1
	  xmax=154
	  xoff=27
	else if(mode.eq.16.and.ic.eq.0)then
	  ymin=1
	  ymax=22
	  xmin=1
	  xmax=77
	  xoff=27
	endif

	arg =dexp(dfloat(1-iframe  )/4.9D0)		!calculate correction timeconstant for dusk
	arg1=dexp(dfloat(iframe-601)/4.9D0)		!and dawn
	if(night.eq.42)  dusk= 16100.D0*arg + 3300.D0*arg1
	if(night.eq.41)  dusk= 24000.D0*arg + 4000.D0*arg1
	if(night.eq.40)  dusk=235000.D0*arg + 1500.D0*arg1

	do jj=ymin,ymax
	do ii=xmin,xmax
	     iairplane=0
	     response = frame(ii,jj)

cc	     find unbinned ccd coordinates
cc	     -1.0 because ccd pixels should start 1 less to make compliant with old C version
	     i = dfloat( ii*mode ) - dfloat(mode-1)/2.D0 - 1.D0 + xoff
	     j = dfloat( jj*mode ) - dfloat(mode-1)/2.D0 - 1.D0

	     x(1) = i - dfloat(mode)/2.D0
	     if(mode.gt.1.and.ii.eq.xmin)x(1) = x(1) + 1.D0	!fix geometry of stamped pixel
	     y(1) = j - dfloat(mode)/2.D0
	     x(2) = x(1)
	     y(2) = j + dfloat(mode)/2.D0
	     y(3) = y(2)
	     x(3) = i + dfloat(mode)/2.D0
	     if(mode.gt.1.and.ii.eq.xmax)x(3) = x(3) - 1.D0	!fix geometry of stamped pixel
	     x(4) = x(3)
	     y(4) = y(1)
	     
	     if(ic.eq.0)then
cc	       calculate ccd pixel bounds within binned pixel:
cc	       be careful of rounding and chopping and precision(now assuming chopping only->adding 1.0 to a fraction to be chopped effectively rounds up)
               itempmin=x(1)+1.0
	       itempmax=x(3)
	       jtempmin=y(1)+1.0
	       jtempmax=y(3)

	       !!test pixel for airplanes
	       iairplane=0
	       do itemp=itempmin,itempmax
	         iairplane=iairplane+airplane(night,iframe,itemp)
	       enddo
	     endif

	     !!this rccd/dr are only preliminary quantities for determining
	     !!if the pixel is in the FOV
	     rccd = dsqrt( (i-x0(ic+1))*(i-x0(ic+1)) + (j-y0(ic+1))*(j-y0(ic+1)) )
	     dr = dabs(rccd - r0)

	     !!less than 1.5 deg x 18.0226 pix/deg
	     if( iairplane.eq.0.and.dr.lt.27.0339D0 )then
	      if( (dabs( (i-x0(ic+1))/(j-y0(ic+1)) ) ).lt.0.57735D0 )then !!less than tan 30

	       if(ic.eq.0)then
	         response=0.0
	         rcount=0
	         do jtemp=jtempmin,jtempmax
		   arg = rccd/dabs(dfloat(jtemp)-y0(ic+1))         !!x secant zenith angle
		   arg1 = (rccd-r0-27.D0)/54.D0
      		   arg=arg*(1.D0 - 0.2D0*arg1)         !!ad hoc E-W dawn/dusk gradient
		   do itemp=itempmin,itempmax
		     response = response + (frame(itemp+1,jtemp+1)-dusk*arg)*(rccd/r0)  !!correct for the 1/r pixel size in sky area, remove dawn/dusk
		     rcount = rcount + 1
		   enddo
	         enddo
	         response=response/dfloat(rcount)
	       endif

cc	       loop for each vertex of a pixel
	       czsum=0.D0
	       do k=1,4
	         !!affine transformation to std. sky coordinates
		 tempx = affine(1) + (1.D0+affine(2))*x(k) + affine(3)*y(k);
		 y(k) = affine(4) + affine(5)*x(k) + (1.D0+affine(6))*y(k);
		 x(k) = tempx;
	         rccd = dsqrt((x(k) - x0(ic+1))*(x(k) - x0(ic+1)) + (y(k) - y0(ic+1))*(y(k) - y0(ic+1)));
	         thetax = dasind( (x(k) - x0(ic+1)) / rccd );
	         thetax = thetax / 0.999744D0 !!atmosphere adjustment

	         !!quadratic coefficiants for the formula: rq=aqTY^2+bqTY+cq
	         !!from July 11, 2001 memo(ignoring 3rd and 4th order terms)

cc	         thetaxsq = thetax*thetax
cc	         aq = 0.1164D0 - (10.D0**-5)*thetaxsq )
cc	         bq = ( 3.6D0*(10.D0**-5)*thetaxsq ) - 18.0266D0
cc	         cq = ( 0.000936D0*thetaxsq ) - ( 1.044D0*(10.D0**-7)*thetaxsq*thetaxsq )

	         rq = (rccd - r0)/1.013D0 !!y-plate scale adjustment
	         thetay = rq / 18.0226D0  !!note: +18.0226 to flip thetay/cy sign relative to Andy's program
cc
cc		 thetay = ( -bq - sqrt((bq*bq)-(4.D0*aq*(cq-rq))) ) / ( 2.D0*aq )

	         cx = -dsind(thetax)
	         cy = dsind(thetay)
	         crsq = (cx*cx)+(cy*cy)
	         if(crsq.gt.0.D0) then
	           phi = datan2d(cy,cx)
	         else
	           phi = 0.D0
	         endif
	         cz = dsqrt(1.D0-(crsq))
		 czsum = czsum + cz
	         rlat = dasind(cz)

	         call rotate8(alfa, beta, gamma, phi, rlat)
		 ra(k)=phi
		 dec(k)=rlat

	      enddo !!end vertex loop
	      if(czsum.gt.0.0) response=4.D0*response/czsum  !!Here correct surface brightness for aperture effective area
	      
	      CALL INDEXPIXEL(level,node,nodeid,hits,response,ra,dec)

	    endif
	   endif
	  enddo
	enddo !!end pixel by pixel double loop

	return
	end