C+
C NAME:
C	indexframe
C PURPOSE:
C	To index an entire frame of TMO data
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call indexframe(level,TOP_KMAX,node,nodeid,hits,night,iframe,affine,euler)
C INPUTS:
C
C OUTPUTS:
C
C CALLS:
C	rotate8(alfa, beta, gamma, phi,rlat)
C	airplane(night, iframe, i)
C
C PROCEDURE:
C
C MODIFICATION HISTORY:
C	Jan, 2003 Aaron Smith (UCSD)
C-
	subroutine indexframe(level,TOP_KMAX,node,nodeid,hits,night,iframe,affine,euler)
	implicit none
	
	integer*8 TOP_KMAX
	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
	real*8 x0, y0, r0, tempx, phi, rlat, rccd, arg, arg1, response, dusk
	integer*4 airplane, iairplane, night, iframe, i, j
	integer*2 k, level

	x0=630.6D0
	y0=1291.38D0
	r0=1194.5D0
	
	alfa=euler(1)
	beta=euler(2)
	gamma=euler(3)

	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 j=1,350
	do i=22,1263
	     iairplane=0
	     !!test pixel for airplanes
	     iairplane=airplane(night,iframe,i)

	     !!this rccd/dr are only preliminary quantities for determining
	     !!if the pixel is in the FOV
	     rccd = dsqrt((dfloat(i)-x0)*(dfloat(i)-x0) + (dfloat(j)-y0)*(dfloat(j)-y0))  !!watch out when i,j get switched
	     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((dfloat(i)-x0)/(dfloat(j)-y0))).lt.0.57735D0 )then !!less than tan 30

	       arg = rccd/dabs(dfloat(j)-y0)         !!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
	       response = (frame(i,j)-dusk*arg)*(rccd/r0)
	       !!correct for the 1/r pixel size in sky area, remove dawn/dusk

	       x(1) = dfloat(i) - 0.5D0
	       y(1) = dfloat(j) - 0.5D0
	       x(2) = dfloat(i) - 0.5D0
	       y(2) = dfloat(j) + 0.5D0
	       x(3) = dfloat(i) + 0.5D0
	       y(3) = dfloat(j) + 0.5D0
	       x(4) = dfloat(i) + 0.5D0
	       y(4) = dfloat(j) - 0.5D0


	       !!loop for each vertex of a pixel
	       do k=1,4
	         !!affine transformation to std. ccd 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)*(x(k) - x0) + (y(k) - y0)*(y(k) - y0));
	         thetax = dasind( (x(k) - x0) / 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
		 !!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))
		 if(k.eq.0)response=response/cz  !!Here correct surface brightness for aperture effective area
	         rlat = dasind(cz)

cc	         beta = camlat-90.D0
cc	         gamma = (-f*degperminute)-ra0 !!(frame#+1)*degperminute

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

	      enddo !!end vertex loop

	     CALL INDEXPIXEL(level,node,nodeid,hits,response,ra,dec)

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


	return
	end
