C+
C NAME:
C	indexframe_q_test2
C PURPOSE:
C	To index an entire frame of SMEI data using a generalized NxN binning scheme
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,frame,affine,q,badcnt,totcnt,BadPix,
C     &	cutout,cutcnt,cutcounter,pattotal,patdark,coarsesky,coarseskyhits,skytime,orbitfraction)
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	level			spatial index level to be used
C	TOP_KMAX		dimension of the three triangle arrays(node,nodeid, hits)
C	node(TOP_KMAX)		array storing summed response for a given spatial triangle
C	nodeid(TOP_KMAX)	array storing triangle name in spatial index format(uint64 in C++)
C	hits(TOP_KMAX)		array storing amount of frames contibuting to a triangle's response
C	frame(id,jd)		2-D array storing the pixelated ccd picture
C	affine(6)		array storing parameters for affine transformation
C	q(4)			quaternion for coordinate rotation
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	rotate8(alfa, beta, gamma, phi,rlat)
C	rotateq(q, phi, rlat)
C	INDEXPIXEL(level,node,nodeid,hits,response,ra,dec,tcount,ttotal)
C PROCEDURE:
C	find the 4 vertices of a pixel, transform those ccd coordinates to sky coordinates,
C	and send the 4 ra/dec vertices to the C++ INDEXPIXEL routine.
C MODIFICATION HISTORY:
C	Jan-May, 2003 Aaron Smith (UCSD)
C	  note: not compliant with TMO data anymore
C-

	subroutine indexframe(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,image1,image2,image3,
     &	imagecnt1,imagecnt2,imagecnt3,npolarimage,spolarimage,npolarimagecnt,spolarimagecnt,eqmask,npmask,spmask)
	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
	integer*2 eqmask(3600,1200), npmask(800,800), spmask(800,800)
	real*4 node(TOP_KMAX),patdark,patarg,z(3,3),z0,zcut,sumzsqr,skytime(720,360),orbitfraction,a,b,abslope
	real*4 frame(id,jd), yoffset(5), xoffset(5), minr,cutcnt(720,360),cutcounter(720,360),cutout(720,360),pattotal(636,128)
	real*8 x(5), y(5), ra(4), dec(4), euler(4), q(4), affine(6), aq, bq, cq
	real*8 rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq, czsum, drmin, drmax, thetsum, thetarot
	real*8 x0(5), y0(5),x00(5), y00(5), yref(5), r0, tempx, phi, rlat, rccd, arg, arg1, response, dusk, i, j
	real*8 dval, px, py, pz
	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*4 ii, iii,jjj,jj, BadPix(id,jd),coarsesky(720,360), badanglecnt
	integer*4 itemp, jtemp, itempmin, itempmax, jtempmin, jtempmax, rcount
	integer*4 xmin, xmax, ymin, ymax, mode, ic, xoff, id, jd, ira, jdec, xcen, ycen
	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
c	badanglecnt=0

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


	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,a,b)
	     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.
		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
		endif
c1000	format('response = ',f8.2,'  sumzsqr = ',f10.2,'  z(2,2) = ',f7.2,'  z0 = ',f7.2,'  z(2,2)**2/sumzsqr = ',f10.2,
c     &		'  (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)
	     
	     !!Middle of pixel
	     x(5) = (x(3)+x(1))/2.D0
	     y(5) = (y(2)+y(1))/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
ccc	      if( (dabs( (i-x0(ic+1))/(j-y0(ic+1)) ) ).lt.0.1051D0 )then !!less than tan 6

	       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 );
		 
		 !!Ad hoc correction of thetax, fixed along with quaternions to 
		 !!make observed star coordinates match the catalogues
		 if(ic.eq.1)then
		   thetax = thetax*0.9995
		 elseif(ic.eq.2)then
		   thetax = thetax*0.997
		 elseif(ic.eq.3)then
		   thetax = thetax*1.0048 - thetax*thetax*0.00007
		 endif

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

		 !!See thetax correction note a few lines above
		 if(ic.eq.1) then
		   thetay = thetay - thetax*thetax*0.0001D0
		 elseif(ic.eq.2) then
		   thetay = thetay - thetax*thetax*0.00005D0
		 elseif(ic.eq.3) then
		   thetay = thetay - thetax*thetax*0.00016D0
		 endif

	         cx = dsind(thetax)
	         cy = dsind(thetay)
	         crsq = (cx*cx)+(cy*cy)
		 cz = dsqrt(1.D0-(crsq))

		 if(crsq.gt.0.D0) then
	             phi = datan2d(cy,cx)
	  	 else
	             phi = 0.D0
	    	 endif
	  	 rlat = dasind(cz)
		 call rotateq(q, phi, rlat) !!rotate to standard system

		 if(k.eq.5)then
		 	!!Paul's analytic way of getting thetarot
		 	dval = dsqrt(1.D0-cx*cx)

                        px =  0.D0                               ! Unit vector for PSF orientation
                        py =  cz/dval                           ! .. in camera coordinate frame
                        pz = -cy/dval

                        call rotateq_xyz(q,cx,cy,cz)	! Camera --> equatorial frame
                        call rotateq_xyz(q,px,py,pz)

                        thetarot = -datan2d( cy*px-cx*py,pz)         ! Angle with equatorial north
			if(thetarot.lt.0.D0)then
c	 		  print*, thetarot,cx,cy,cz
			  thetarot=thetarot+360.D0
c			  badanglecnt=badanglecnt+1
			elseif(thetarot.gt.360.D0)then
c			  print*, thetarot,cx,cy,cz
			  thetarot=thetarot-360.D0
c			  badanglecnt=badanglecnt+1
			endif

		   CALL index_theta_rot(thetarot,phi,rlat,image1,image2,image3,imagecnt1,imagecnt2,imagecnt3,
     &		     npolarimage,spolarimage,npolarimagecnt,spolarimagecnt)

		   ira=720-int(2.*ra(1))
		   jdec=1+int(180.+2.*dec(1))
		   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
		 else
		   czsum = czsum + cz
		   ra(k)=phi
		   dec(k)=rlat
		 endif

	       enddo !!end vertex loop
	       
	       if(ra(1).gt.340.6.and.ra(1).lt.340.7.and.dec(1).gt.-46.9.and.dec(1).lt.-46.8)then
	         print*, ii,jj, ra(1), dec(1)
	       endif

	       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
c	       if(ii.le.250)then
		 !!See if location is close to a star (0=nostar, 1=star) and get a slope if it is,
		 !!if not set slope=-1 sending the pixel to the reguler crx algorithm
	         if(dec(1).gt.59.5)then
		   xcen = nint( 10.0*(90.0-dec(1))*cosd(-ra(1)) + 400.5 )
		   ycen = nint( 10.0*(90.0-dec(1))*sind(-ra(1)) + 400.5 )
		   if(npmask(xcen,ycen).eq.1) then
		     do iii=1,3
	     	      do jjj=1,3
	     		z(iii,jjj)=frame(ii+iii-2,jj+jjj-2)
	     	      enddo
	     	     enddo
		     call plane(z,z0,a,b)
		     abslope=sqrt(a*a+b*b)
		   else
		     abslope=-1.0
		   endif
		 elseif(dec(1).lt.-59.5)then
		   xcen = nint( -10.0*(90.0+dec(1))*cosd(ra(1)) + 400.5 )
		   ycen = nint( -10.0*(90.0+dec(1))*sind(ra(1)) + 400.5 )
		   if(spmask(xcen,ycen).eq.1) then
		     do iii=1,3
	     	      do jjj=1,3
	     		z(iii,jjj)=frame(ii+iii-2,jj+jjj-2)
	     	      enddo
	     	     enddo
		     call plane(z,z0,a,b)
		     abslope=sqrt(a*a+b*b)
		   else
		     abslope=-1.0
		   endif
		 else
		   xcen = nint( -(ra(1)-180.0)*10.0 + 1800.5 )
		   ycen = nint( dec(1)*10.0 + 600.5 )
		   if(eqmask(xcen,ycen).eq.1) then
		     do iii=1,3
	     	      do jjj=1,3
	     		z(iii,jjj)=frame(ii+iii-2,jj+jjj-2)
	     	      enddo
	     	     enddo
		     call plane(z,z0,a,b)
		     abslope=sqrt(a*a+b*b)
		   else
		     abslope=-1.0
		   endif
		 endif
		 !!Use this if-block to only index stars within the starmask, and within ra/dec bounds
c		 if(abslope.gt.0.0.and.dec(1).lt.59.0.and.dec(1).gt.-59.0) then
		   CALL INDEXPIXEL(level,node,nodeid,hits,response,ra,dec,tcount,ttotal,iflag,mflag,a,b,abslope)
c		 endif
ccc		 CALL INDEXPIXELNOALG(level,node,nodeid,hits,response,ra,dec,tcount)
c	       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
	       coarsesky(ira,jdec)=coarsesky(ira,jdec)+int(response)
	       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
	       coarseskyhits(ira,jdec)=coarseskyhits(ira,jdec)+1
	     endif
	     endif !!end FOV if-block
	   endif  !!end non-zero response if-block
	enddo
	enddo !!end pixel by pixel double loop
	
c	if(badanglecnt.gt.0)  print*, '		******* Bad Angles: ',badanglecnt
	
	return
	end