C+
C NAME:
C	indexframe_centroid
C PURPOSE:
C	index a small box around a star
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,frame,affine,q, cra,cdec)
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	frame(1280,600)		2-D array storing the pixelated ccd picture
C	affine(6)		array storing parameters for affine transformation
C	q(4)			quaternion for coordinate rotation
C	cra/cdec		coordinates of star to be indexed
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	rotateq(q, phi, rlat)
C	INDEXPIXEL(level,node,nodeid,hits,response,ra,dec)
C PROCEDURE:
C	find ccd coordinates for approximate centroid of a bright star and index a box around those coordinates
C MODIFICATION HISTORY:
C	May, 2003 Aaron Smith (UCSD)
C-
	subroutine indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,frame,affine,qq,cra,cdec,
     &		avgthetarot,fthetax,badangle,xra,ydec,goodframe)
	implicit none

	integer*8 TOP_KMAX, kk, frametest
	real*4 frame(id,jd), yoffset(5), xoffset(5)
	real*4 node(TOP_KMAX)
	integer*8 nodeid(TOP_KMAX)
	integer*2 hits(TOP_KMAX), goodframe
	real*8 x(4), y(4), ra(4), dec(4), euler(4), q(4), qq(4), affine(6)
	real*8 aq, bq, cq, xccd, yccd,thetarot,avgthetarot, xra(2), ydec(2), drmin, drmax
	real*8 rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq, czsum, cra, cdec
	real*8 x0(5), y0(5),x00(5), y00(5), yref(5), r0, tempx, phi, rlat, rccd, arg, arg1, response, dusk, i, j
	real*8 alfa, beta, gamma, xtri, ytri, hyptri, slope, fthetax
	integer*4 ii, jj
	integer*4 itemp, jtemp, itempmin, itempmax, jtempmin, jtempmax, rcount
	integer*4 xmin, xmax, ymin, ymax, mode, ic, xoff, id, jd, badangle
	integer*2 k, level, zz

C	numbers from Andy's memo
	data yoffset /0.,65.,59.,61.,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/
	if(ic.eq.1) r0=1193.5D0
	if(ic.eq.2) r0=1188.5D0
	if(ic.eq.3) r0=1198.5D0
	
	goodframe=0

	!!Ad hoc correction to x0 for camera 2
	if(ic.eq.2)then
	  x0(ic+1)=x00(ic+1)-xoffset(ic+1)-2.D0
	  y0(ic+1)=y00(ic+1)-yoffset(ic+1)
	else
	  x0(ic+1)=x00(ic+1)-xoffset(ic+1)
	  y0(ic+1)=y00(ic+1)-yoffset(ic+1)
	endif

	if(ic.eq.1.or.ic.eq.2)then
	    drmin=-27.D0
	    drmax=22.D0		!was 32.D0: should be 22.D0 according to 11 July 2001 memo, page 5
	elseif(ic.eq.3)then
	    drmin=-26.D0
	    drmax=20.D0
c	    drmin=-34.D0	!!-35.0
c	    drmax=24.D0		!! 24.0
	endif

CC	Find ccd coordinates of star, don't change cra/cdec!!!
CC	Take inverse quaternion to go from sky to ccd
	  q(1)=qq(1)
	  q(2)=-qq(2)
	  q(3)=-qq(3)
	  q(4)=-qq(4)
	  phi=cra
	  rlat=cdec
	  call rotateq_sky_to_ccd(q,phi,rlat, cx, cy, cz)

	  thetax=dasind(cx)
	  fthetax=thetax
	  thetay=dasind(cy)
	  rccd=-18.0226D0*thetay + r0
	  dr=rccd-r0
	  xccd=rccd*dsind(thetax) + x0(ic+1)
	  yccd=y0(ic+1)-rccd*dcosd(thetax)
	  
c	  print*, 'Centroid: ',xccd,yccd

CC	find pixel bounds for box to be drawn around star, don't go outside array bounds
	  ymin=idnint(yccd)-60
	  if(ymin.lt.1)ymin=1
	  ymax=idnint(yccd)+60
	  if(ymax.gt.256)ymax=256
	  xmin=idnint(xccd)-60
	  if(xmin.lt.22)xmin=22
	  xmax=idnint(xccd)+60
	  if(xmax.gt.1263)xmax=1263

CC      find angle of intersection between a line of constant thetax and a line of longitude
	  badangle=0
	  call find_theta_rot(mode,ic,qq,cra,cdec,dr,thetax,avgthetarot,badangle,xra,ydec)

c	ii=xccd
c	jj=yccd
	do jj=ymin,ymax
	do ii=xmin,xmax

	     response = frame(ii,jj)

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

	     i = dfloat(ii) - 1.D0
	     j = dfloat(jj) - 1.D0

C	     find 4 vertex coordinates for pixel
	     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)

	     !! these 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 = rccd - r0

	     !!this dr test is for cameras 1, 2, and 3
	     if( dr.lt.drmax.and.dr.gt.drmin)then
	      if( (dabs( (i-x0(ic+1))/(j-y0(ic+1)) ) ).lt.0.57735D0 )then !!less than tan 30

	       response = response*(rccd/r0)  !!correct for the 1/r pixel size in sky area

	       !! 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)));
	         rq = (rccd - r0) !!/1.013D0 !!y-plate scale adjustment
		 thetax = dasind( (x(k) - x0(ic+1)) / rccd );
		 
		 !!Ad hoc correction of thetax, fixed along with quaternions to
		 !!make observed star coordinates match the catalogues
		 if(ic.eq.1)then
		   thetax = thetax*0.9995
		 elseif(ic.eq.2)then
		   thetax = thetax*0.997
		 elseif(ic.eq.3)then
		   thetax = thetax*1.0048 - thetax*thetax*0.00007
		 endif

		 !! quadratic coefficients from Andy's July 2001 memo
		 aq = 0.1664D0 - (0.00001*thetax*thetax)
		 bq = (0.000036D0*thetax*thetax) - 18.0226D0
		 cq = (0.000936D0*thetax*thetax) - rq
		 thetay = ( -bq-sqrt((bq*bq)-(4.D0*aq*cq)) ) / ( 2.D0*aq )
C	         thetay = rq / -18.0226D0

		 !!See thetax correction note a few lines above
		 if(ic.eq.1) then
		   thetay = thetay - thetax*thetax*0.0001D0
		 elseif(ic.eq.2) then
		   thetay = thetay - thetax*thetax*0.00005D0
		 elseif(ic.eq.3) then
		   thetay = thetay - thetax*thetax*0.00016D0
		 endif

	         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 rotateq(qq, phi, rlat) !!rotate to standard system using original qq quaternion
		 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

	      if(nint(i).eq.nint(xccd).and.nint(j).eq.nint(yccd)) goodframe=1
	      CALL INDEXPIXEL(level,node,nodeid,hits,response,ra,dec)

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

	return
	end