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 find_theta_rot(mode,ic,q,ra,dec,dr,thetax,thetarot,badangle)
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	ra/dec			coordinates of star
C	thetarot		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 two 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	Final thetarot should be clockwise angle of rotation for fishtail to be pointed north
C MODIFICATION HISTORY:
C	Aug, 2003 Aaron Smith (CASS)
C	Aug, 2004 Aaron Smith (CASS)
C-
	subroutine index_theta_rot(thetarot,ra,dec,image1,image2,image3,imagecnt1,imagecnt2,imagecnt3,
     &    npolarimage,spolarimage,npolarimagecnt,spolarimagecnt)
	implicit none

	real*8 thetarot, crsq, ra, dec, arg, cra1, cra2, cra3
	integer*4 xarg, yarg
	real*4 image1(120,120), image2(120,120), image3(120,120)
	real*4 npolarimage(160,160), spolarimage(160,160)
	integer*4 imagecnt1(120,120), imagecnt2(120,120), imagecnt3(120,120)
	integer*4 spolarimagecnt(160,160), npolarimagecnt(160,160)
	integer*2 ibad

	if(thetarot.gt.360.0.or.thetarot.lt.0.0) print*, '************** Bad Angle! ***********'

cccccccccccccccccccccccccccc  Fill Coarse Rotator Maps  ccccccccccccccccccccccccccccccccccccccc
	cra1=60.D0
	cra2=180.D0
	cra3=300.D0

c
c	   Make first RA/dec plot (0-120 in ra)
c
	      ibad=0
	      xarg = nint( -(ra-cra1) + 60.5D0 )
	      yarg = nint( dec + 60.5D0 )
	      if(xarg.lt.1.or.xarg.gt.120.or.yarg.lt.1.or.yarg.gt.120)ibad=1
	      if(ibad.eq.0)then
	        !!At the border between 0 and 360 degrees averaging gives bad values,
		!!so only use one of the values
		if(imagecnt1(xarg,yarg).ge.1)then
		  arg=image1(xarg,yarg)/float(imagecnt1(xarg,yarg))
		  if(arg.gt.320.00.and.thetarot.lt.180.00) go to 976
		  if(arg.lt.40.00.and.thetarot.gt.180.00) go to 976
		endif
	        image1(xarg,yarg) = image1(xarg,yarg) + thetarot
	        imagecnt1(xarg,yarg) = imagecnt1(xarg,yarg) + 1
	      endif
c
c	   Make second RA/dec plot (120-240 in ra)
c
	      ibad=0
	      xarg = nint( -(ra-cra2) + 60.5D0 )
	      yarg = nint( dec + 60.5D0 )
	      if(xarg.lt.1.or.xarg.gt.120.or.yarg.lt.1.or.yarg.gt.120)ibad=1
	      if(ibad.eq.0)then
	        !!At the border between 0 and 360 degrees averaging gives bad values,
		!!so only use one of the values
		if(imagecnt2(xarg,yarg).ge.1)then
		  arg=image2(xarg,yarg)/float(imagecnt2(xarg,yarg))
		  if(arg.gt.320.00.and.thetarot.lt.180.00) go to 976
		  if(arg.lt.40.00.and.thetarot.gt.180.00) go to 976
		endif
	        image2(xarg,yarg) = image2(xarg,yarg) + thetarot
	        imagecnt2(xarg,yarg) = imagecnt2(xarg,yarg) + 1
	      endif
c
c	   Make third RA/dec plot (240-360 in ra)
c
	      ibad=0
	      xarg = nint( -(ra-cra3) + 60.5D0 )
	      yarg = nint( dec + 60.5D0 )
	      if(xarg.lt.1.or.xarg.gt.120.or.yarg.lt.1.or.yarg.gt.120)ibad=1
	      if(ibad.eq.0)then
	        !!At the border between 0 and 360 degrees averaging gives bad values,
		!!so only use one of the values
		if(imagecnt3(xarg,yarg).ge.1)then
		  arg=image3(xarg,yarg)/float(imagecnt3(xarg,yarg))
		  if(arg.gt.320.00.and.thetarot.lt.180.00) go to 976
		  if(arg.lt.40.00.and.thetarot.gt.180.00) go to 976
		endif
	        image3(xarg,yarg) = image3(xarg,yarg) + thetarot
	        imagecnt3(xarg,yarg) = imagecnt3(xarg,yarg) + 1
	      endif
c
c	   Make North Polar plot
c
	    if(dec.gt.33.D0)then
	      ibad=0
	      xarg = nint( 2.D0*(90.D0-dec)*cosd(-ra) + 80.5D0)
	      yarg = nint( 2.D0*(90.D0-dec)*sind(-ra) + 80.5D0)
	      if(xarg.lt.1.or.xarg.gt.160.or.yarg.lt.1.or.yarg.gt.160)ibad=1
	      if(ibad.eq.0)then
	        !! Convert eq angle to north polar angle
	        thetarot=360.D0-(thetarot+90.D0-ra)
		if(thetarot.gt.360.0) then
		  thetarot=thetarot-360.D0
		elseif(thetarot.lt.0.0)then
		  thetarot=thetarot+360.D0
		endif
		!!At the border between 0 and 360 degrees averaging gives bad values,
		!!so only use one of the values
		if(npolarimagecnt(xarg,yarg).ge.1)then
		  arg=npolarimage(xarg,yarg)/float(npolarimagecnt(xarg,yarg))
		  if(arg.gt.320.00.and.thetarot.lt.180.00) go to 976
		  if(arg.lt.40.00.and.thetarot.gt.180.00) go to 976
		endif
	        npolarimage(xarg,yarg) = npolarimage(xarg,yarg) + thetarot
	        npolarimagecnt(xarg,yarg) = npolarimagecnt(xarg,yarg) + 1
	      endif
	    endif
c
c	   Make South Polar plot
c
	    if(dec.lt.-33.D0)then
	      ibad=0
	      xarg = nint( -2.D0*(90.D0+dec)*cosd(ra) + 80.5D0)
	      yarg = nint( -2.D0*(90.D0+dec)*sind(ra) + 80.5D0)
	      if(xarg.lt.1.or.xarg.gt.160.or.yarg.lt.1.or.yarg.gt.160)ibad=1
	      if(ibad.eq.0)then
	        !! Convert eq angle to south polar angle
	        thetarot=360.D0-(thetarot+ra-270.D0)
		if(thetarot.gt.360.0) then
		  thetarot=thetarot-360.D0
		elseif(thetarot.lt.0.0)then
		  thetarot=thetarot+360.D0
		endif
		!!At the border between 0 and 360 degrees averaging gives bad values,
		!!so only use one of the values
		if(spolarimagecnt(xarg,yarg).ge.1)then
		  arg=spolarimage(xarg,yarg)/float(spolarimagecnt(xarg,yarg))
		  if(arg.gt.320.00.and.thetarot.lt.180.00) go to 976
		  if(arg.lt.40.00.and.thetarot.gt.180.00) go to 976
		endif
	        spolarimage(xarg,yarg) = spolarimage(xarg,yarg) + thetarot
	        spolarimagecnt(xarg,yarg) = spolarimagecnt(xarg,yarg) + 1
	      endif
	    endif

976	return
	end