C+
C NAME:
C	indexframe_rotator.f
C MODIFICATION HISTORY:
C	Aug, 2004 Aaron Smith (UCSD)
C-

	subroutine indexframe_rotator(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,frame,affine,q,badcnt,totcnt,BadPix,
     &	cutout,cutcnt,cutcounter,pattotal,patdark,coarsesky,coarseskyhits,skytime,orbitfraction)
	implicit none

	integer*8 TOP_KMAX, kk
	integer*8 nodeid(TOP_KMAX),tcount, ttotal, totcnt, badcnt
	integer*2 hits(TOP_KMAX), coarseskyhits(720,360),iflag,mflag
	real*4 node(TOP_KMAX),patdark,patarg,z(3,3),z0,zcut,sumzsqr,skytime(720,360),orbitfraction
	real*4 frame(id,jd), yoffset(5), xoffset(5), minr,cutcnt(720,360),cutcounter(720,360),cutout(720,360),pattotal(636,128)
	real*8 thetarot, xra(2), ydec(2), cra, cdec, x(5), y(5), ra(5), dec(5),ra4(4), dec4(4)
	real*8 euler(4), q(4), affine(6), aq, bq, cq
	real*8 rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq, czsum, drmin, drmax, thetsum
	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, iii,jjj,jj, BadPix(id,jd),coarsesky(720,360), badangle
	integer*4 itemp, jtemp, itempmin, itempmax, jtempmin, jtempmax, rcount
	integer*4 xmin, xmax, ymin, ymax, mode, ic, xoff, id, jd, ira, jdec
	integer*2 k, level

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.5D0
	if(ic.eq.2) r0=1188.5D0
	if(ic.eq.3) r0=1198.5D0
	totcnt=0
	badcnt=0
	iflag=0
	mflag=0

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

	if(mode.eq.1)then
	  ymin=1
	  ymax=jd
	  xmin=22
	  xmax=1263 !!(1264-1) to remove the last half covered pixel
	  xoff=0
	  zcut=111.	!500 electron cut
	  if(ic.eq.1.or.ic.eq.2)then
	    drmin=-27.D0
	    drmax=22.D0		!was 32.D0: should be 22.D0 according to 11 July 2001 memo, page 5
	  elseif(ic.eq.3)then
	    drmin=-26.D0
	    drmax=20.D0
	  endif
	else if(mode.eq.2)then
	  ymin=1
	  ymax=jd
	  xmin=11
	  xmax=632
	  xoff=0
	  zcut=27.76	!500 electron cut
	  if(ic.eq.1.or.ic.eq.2)then
	    drmin=-27.D0
	    drmax=22.D0
	  elseif(ic.eq.3)then
	    drmin=-26.D0
	    drmax=20.D0
	  endif
	else if(mode.eq.4)then
	  ymin=1
	  ymax=jd
	  xmin=6
	  xmax=316
	  xoff=0
	  zcut=6.94	!500 electron cut
	  if(ic.eq.1.or.ic.eq.2)then
	    drmin=-26.D0
	    drmax=21.D0		!was 31.D0: see note above
	  elseif(ic.eq.3)then
	    drmin=-25.D0
	    drmax=19.D0
	  endif
	endif

	do jj=ymin,ymax
	do ii=xmin,xmax
	  iflag=0
	  if(ic.eq.3)then
	     patarg=pattotal(ii,jj)-patdark*0.4*(1.+0.01*float(jj))	!remove dark current slope from hot/flipper pixel flag
	     if(patarg.gt.50.) iflag=1
	  endif
	  mflag=0
	  response = frame(ii,jj)
	  if(response.gt.0.)then
	     if(ic.eq.3.and.ii.gt.1.and.ii.lt.xmax.and.jj.gt.1.and.jj.lt.ymax)then	!? apply to cam 1 & 2?
	     do iii=1,3
	     do jjj=1,3
	     	z(iii,jjj)=frame(ii+iii-2,jj+jjj-2)
	     enddo
	     enddo
	     call plane(z,z0)
	     if(abs(z(2,2)).gt.zcut)then
	     	sumzsqr=0.
	     	do iii=1,3
	     	do jjj=1,3
	     	if(iii.ne.2.or.jjj.ne.2)sumzsqr=sumzsqr+z(iii,jjj)**2
	     	enddo
	     	enddo
	     	sumzsqr=sumzsqr/5.
c		if((z(2,2)*z(2,2))/sumzsqr.gt.5.)then
c		print 1000,response,sumzsqr,z(2,2),z0,z(2,2)**2/sumzsqr,ii,jj
c		endif
1000	format('response = ',f8.2,'  sumzsqr = ',f10.2,'  z(2,2) = ',f7.2,'  z0 = ',f7.2,'  z(2,2)**2/sumzsqr = ',f10.2,
     &		'  (i,j) = (',i3,',',i3,')')
		if((z(2,2)*z(2,2))/sumzsqr.gt.25.)mflag=1
		if(iflag.eq.1.and.z(2,2).lt.0.)mflag=1
	     endif
	     endif
CC	     find unbinned ccd coordinates
CC	     -1.0 because ccd pixels should start 1 less to make compliant with old C version
	     i = dfloat( ii*mode ) - dfloat(mode-1)/2.D0 - 1.D0 + xoff
	     j = dfloat( jj*mode ) - dfloat(mode-1)/2.D0 - 1.D0

CC	     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)
	     
	     x(5)=( x(1) + x(3) ) / 2.D0
	     y(5)=( y(1) + y(3) ) / 2.D0

CC	     this rccd/dr are only preliminary quantities for determining
CC	     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 1, 2, and 3
	     if( dr.lt.drmax.and.dr.gt.drmin)then
	      if( (dabs( (i-x0(ic+1))/(j-y0(ic+1)) ) ).lt.0.57735D0 )then !!less than tan 30


	       response=response*(rccd/r0)	!first order radial binsize correction

cc	       loop for each vertex of a pixel
	       czsum=0.D0
	       thetsum=0.D0
	       do k=1,5

	         !!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)
		 thetax = dasind( (x(k) - x0(ic+1)) / rccd );

		 !! quadratic coefficients 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 )
		 thetsum=thetsum+thetay
C	         thetay = rq / -18.0226D0

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

		 call rotateq(q, phi, rlat) !!rotate to standard system

		 if(k.lt.5)then
		   ra4(k)=phi
		   dec4(k)=rlat
		 else
		   cra=phi
		   cdec=rlat
		   thetarot=0.D0
		   call find_theta_rot(mode,ic,q,cra,cdec,rq,thetax,thetarot,badangle)
		   ira=720-int(2.*cra)
		   jdec=1+int(180.+2.*cdec)
		   if(ira.gt.720) then
		     ira=720
		   elseif(ira.lt.1)then
		     ira=1
		   endif
		   if(jdec.gt.360)then
		     jdec=360
		   elseif(jdec.lt.1)then
		     jdec=1
		   endif
		 endif

	       enddo !!end k vertex loop

	       response=response*(1.D0-0.0092D0*thetsum/4.D0)	!radial binsize, nonlinear term
	       if(czsum.gt.0.0) response=4.D0*response/czsum  !!Here correct surface brightness for aperture effective area

	       tcount=0
	       ttotal=0
	       if(badangle.ne.1)then
		 CALL INDEXPIXELNOALG(level,node,nodeid,hits,thetarot,ra4,dec4,tcount)
	       endif
	       BadPix(ii,jj)=BadPix(ii,jj)+tcount
	       totcnt=totcnt+ttotal
	       badcnt=badcnt+tcount
	       if(ttotal.ne.0)then
	         cutout(ira,jdec)=cutout(ira,jdec)+response*float(tcount)/float(ttotal)
	         cutcnt(ira,jdec)=cutcnt(ira,jdec)+float(tcount)/float(ttotal)
	         cutcounter(ira,jdec)=cutcounter(ira,jdec)+1.
	       endif
	       if(badangle.ne.1)then
	         !!At the border between 0 and 360 degrees averaging gives bad values,
		 !!so only use one of the values
		 if(coarseskyhits(ira,jdec).ge.1)then
		   arg=float(coarsesky(ira,jdec))/float(coarseskyhits(ira,jdec))
		   if(arg.gt.355.00.and.thetarot.lt.180.00) go to 976
		   if(arg.lt.5.00.and.thetarot.gt.180.00) go to 976
		 endif
	         coarsesky(ira,jdec)=coarsesky(ira,jdec)+thetarot
	         coarseskyhits(ira,jdec)=coarseskyhits(ira,jdec)+1
	       endif
976	       arg=skytime(ira,jdec)/float(coarseskyhits(ira,jdec))
	       if(orbitfraction/arg.gt.30.)skytime(ira,jdec)=-skytime(ira,jdec)
	       if(skytime(ira,jdec).lt.0.)skytime(ira,jdec)=skytime(ira,jdec)-orbitfraction
	       if(skytime(ira,jdec).ge.0.)skytime(ira,jdec)=skytime(ira,jdec)+orbitfraction
	     endif
	     endif !!end FOV if-block
	   endif  !!end non-zero response if-block
	enddo
	enddo !!end pixel by pixel double loop

	return
	end
