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,xmin,xmax,tan1, tan2)
	implicit none

	integer*8 TOP_KMAX, kk, counter
	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(4), q(4), qq(4), affine(6)
	real*8 aq, bq, cq, xccd, yccd, tan1, tan2
	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
	integer*4 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, 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.2) r0=1188.5
	if(ic.eq.3) r0=1198.5

C	numbers for camera 1,2,3 from SMEI Point v.1.0
C	data yoffset /0.,0.,0.,0.,55./
C	data yref /0.,61.612891D0,61.579259D0,59.0449220,0./
C	data x00 /630.D0,631.96653D0,625.78564D0,627.96738D0,633.D0/
C	data y00 /1291.38D0,1241.6218D0,1241.4524D0,1237.6987D0,1294.3D0/
C	r0=y00(ic+1)-yref(ic+1)

	x0(ic+1)=x00(ic+1)-xoffset(ic+1)
	y0(ic+1)=y00(ic+1)-yoffset(ic+1)

	counter=0

	do jj=1,jd
	do ii=22,1263

	     response = 1000.D0

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 2 and 3
	     if( (ic.eq.2.and.dr.lt.32.4D0.and.dr.gt.-27.D0).or.(ic.eq.3.and.dr.lt.30.D0.and.dr.gt.-26.D0) )then
	      if( (dabs( (i-x0(ic+1))/(j-y0(ic+1)) ) ).lt.0.57735D0 )then !!less than tan 30

	     if(dr.lt.9.D0.and.dr0.gt.-9.D0)then
	      if( ((dabs( (i-x0(ic+1))/(j-y0(ic+1)) ) ).lt.tan1).and.((dabs( (i-x0(ic+1))/(j-y0(ic+1)) ) ).gt.tan2) ) then

	       counter=counter+1

	       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 );

		 !! quadratic coeffitients 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

	         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
c		 if(k.eq.1) print*, ra(k),dec(k),thetax

	      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
	   endif
	   endif
	  enddo
	enddo !!end pixel by pixel double loop
	
	print*, 'Count: ',counter

	return
	end