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, j, 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.,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.5
	if(ic.eq.2) r0=1188.5
	if(ic.eq.3) r0=1198.5

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

	ramax=cra+80.D0
	ramin=cra-80.D0
	ra=cra
	dec=cdec

C	calulate center of FOV
	phic=0
	rlatc=90
	call rotateq(qq, phic, rlatc)

C	print*,'Center RA/DEC: ',phic,rlatc

C	make sure star has a possibility of being in FOV
C	note: without this check, ra+180 will appear to be in FOV, making
C	      it seem as though the star transits twice an orbit -> WRONG!!!
	if(phic.gt.ramin.and.phic.lt.ramax) then

	  !!find ccd coordinates of star
	  call rotateq_sky_to_ccd(q,ra,dec,cx,cy,cz)

c	  cz=dsind(dec)
c	  crsq=1.D0-(cz*cz)
c	  cx=dsqrt( crsq/(1.D0+(dtand(ra)*dtand(ra))) )
c	  cy=dsqrt(crsq-(cx*cx))
c	  if(dtand(ra).lt.0.)cy=-cy

c	  print*, 'Cx : ',cx
c	  print*, 'Cy : ',cy
c	  print*, 'Cz : ',cz

	  thetax=dasind(cx)
	  thetay=dasind(cy)
c	  print*, 'Tx, Ty: ',thetax, thetay
	  rccd=-18.0226D0*thetay + r0

c	  print*, 'Rccd: ',rccd

	  xccd=rccd*dsind(thetax) + x0(ic+1)
	  yccd=y0(ic+1)-rccd*dcosd(thetax)

c	  print*, 'X,Y: ',xccd, yccd
c	  xccd=633
c	  yccd=47

	ii=idnint(xccd)
	jj=idnint(yccd)

	i=xccd
	j=yccd

	     !! 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.or.ic.eq.1).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

C	       print*,'		Centroid: ',xccd,yccd
	       ibad=1

		 rccd = dsqrt((ii - x0(ic+1))*(ii - x0(ic+1)) + (jj - y0(ic+1))*(jj - y0(ic+1)));

c		 print*, 'Rccd: ', rccd

	         rq = (rccd - r0)/1.013D0 !!y-plate scale adjustment
		 thetax = dasind( (ii - x0(ic+1)) / rccd );
		 thetay = rq / -18.0226D0

c		 print*, 'ThetaX, ThetaY: ',thetax, thetay

	         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)

c		 print*, 'Cx : ',cx
c		 print*, 'Cy : ',cy
c		 print*, 'Cz : ',cz
		 
		 call rotateq(qq, ra, dec) !!rotate to standard system using original qq quaternion

C		 print*, '		RA/DEC_: ',ra,dec

	     endif
	    endif

	endif

	return
	end