C+
C NAME:
C	indexframe_transit_test.f
C PURPOSE:
C	find the frames containing a star
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call indexframe(mode,ic, id, jd,affine,qq, ibad, 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	affine(6)		array storing parameters for affine transformation
C	q(4)			quaternion for coordinate rotation
C	cra/cdec		star coordinates for
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,affine,qq, ibad, cra, cdec)
	implicit none

	real*4 yoffset(5), xoffset(5)
	real*8 euler(4), q(4), qq(4), affine(6), xccd, yccd, cra, cdec
	real*8 rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq
	real*8 x0(5), y0(5),x00(5), y00(5), yref(5)
	real*8 r0, rccd, i(3), j(3), phic, rlatc, ra, dec, ramax, ramin
	integer*4  k, ii, jj, mode, ic, xoff, id, jd, ibad

C	numbers from Andy's memo
	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

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)

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

	i(1) = 50.D0
	i(2) = 640.D0
	i(3) = 1200.D0
	j(1) = 200.D0
	j(2) = 25.D0
	j(3) = j(1)

	   print*, 'Entering Loop...'

	   do ii=1,3

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

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

	       ibad=0
	       
	         rccd = dsqrt((i(ii) - x0(ic+1))*(i(ii) - x0(ic+1)) + (j(ii) - y0(ic+1))*(j(ii) - y0(ic+1)));
	         rq = (rccd - r0)/1.013D0 !!y-plate scale adjustment
		 thetax = dasind( (i(ii) - x0(ic+1)) / rccd );

		 print*, ii,'	ThetaX: ',thetax

		 thetay = rq / -18.0226D0

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

		 call rotateq(qq, ra, dec) !!rotate to standard system using original qq quaternion
		 
		 print*, ii,'	RA/DEC_: ',ra,dec

	     endif
	    endif
	    
	  enddo

	return
	end