C+
C NAME:
C	find_theta_rot
C PURPOSE:
C	find the angle of intersection between a line on the ccd of constant thetax, and a line of longitude on the globe
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	q(4)			quaternion for coordinate rotation
C	cra/cdec		coordinates of star
C	avgthetarot		sum of intersection angle, to be averaged
C	xra/ydec		ra/dec to be used as cartesian coordinates in finding angle of intersection
C	dr			number of pixels cra/cdec is from r0
C PROCEDURE:
C	find to points on ccd of constant thetax, rotate those coordinates so that cra, cdec become [0,0]
C	then create a line from those two points and find the angle between it and y-axis (dec-axis)
C MODIFICATION HISTORY:
C	August, 2003 Aaron Smith (UCSD)
C-
	subroutine find_theta_rot(mode,ic,q,cra,cdec,dr,thetax,avgthetarot,badangle,xra,ydec)
	implicit none

	real*8 aq, bq, cq, xccd, yccd, thetarot, avgthetarot, xra(2), ydec(2)
	real*8 rq, dr, thetax, thetay, cx, cy, cz, crsq, cra, cdec
	real*8 r0, tempx, phi, rlat, rccd, q(4)
	real*8 alfa, beta, gamma, xtri, ytri, hyptri, slope
	integer*4 ii, jj, mode, ic, badangle
	integer*2 k, level, zz

	if(ic.eq.1) r0=1193.5D0
	if(ic.eq.2) r0=1188.5D0
	if(ic.eq.3) r0=1198.5D0

	badangle=0
	do ii=1,2
	     	if(ii.eq.1) rq=dr+9.D0
	     	if(ii.eq.2) rq=dr-9.D0
	     	if(rq.gt.32.4D0.or.rq.lt.-27.D0) badangle=1   !!dont use a line that is outside the ROI

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

	         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))
	         rlat = dasind(cz)

		 call rotateq(q, phi, rlat) !!rotate to standard system using original qq quaternion
		 xra(ii)=phi
		 ydec(ii)=rlat
	enddo

	alfa=cra !!line up with x-axis
	beta=-cdec  !!rotate about new y-axis so centroid now at [0,0]
	gamma=0.D0

	do ii=1,2
		phi=xra(ii)
	 	rlat=ydec(ii)
	  	call rotate8(alfa,beta,gamma,phi,rlat)
	  	if( phi.gt.180.D0 ) phi=phi-360.D0
	  	if( phi.lt.-180.D0 ) phi=phi+360.D0
	  	xra(ii)=phi
	  	ydec(ii)=rlat
	enddo
	slope=(ydec(1)-ydec(2))/(xra(1)-xra(2))

	xtri=xra(1)-xra(2)
	ytri=dabs(ydec(1))+dabs(ydec(2))
	hyptri=dsqrt( xtri*xtri + ytri*ytri )
	thetarot=dacosd(xtri/hyptri)

	if(slope.gt.0.D0)then
	  if(xra(1).lt.xra(2)) thetarot=90.D0-thetarot
	  if(xra(1).gt.xra(2)) thetarot=thetarot-270.D0
	else
	  if(xra(1).lt.xra(2)) thetarot=thetarot-270.D0
	  if(xra(1).gt.xra(2)) thetarot=90.D0-thetarot
	endif
	if(badangle.eq.0) avgthetarot=avgthetarot+thetarot

c	print*, '	      Theta: ',thetarot,avgthetarot

	return
	end