C+ C NAME: C indexframe C PURPOSE: C To index an entire frame of TMO data using an NxN binning scheme C CATEGORY: C Data processing C CALLING SEQUENCE: C call indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,night,iframe,frame,affine,q) 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(1280,600) 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) 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) implicit none integer*8 TOP_KMAX, kk integer*8 nodeid(TOP_KMAX),tcount integer*2 hits(TOP_KMAX) real*4 node(TOP_KMAX) real*4 frame(id,jd), yoffset(5), xoffset(5), minr real*8 x(4), y(4), 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 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, jj integer*4 itemp, jtemp, itempmin, itempmax, jtempmin, jtempmax, rcount integer*4 xmin, xmax, ymin, ymax, mode, ic, xoff, id, jd 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 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 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 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 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 response = frame(ii,jj) if(response.gt.0.)then 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) 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,4 !!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 ra(k)=phi dec(k)=rlat enddo !!end 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 CALL INDEXPIXEL(level,node,nodeid,hits,response,ra,dec) endif endif !!end FOV if-block endif !!end non-zero response if-block enddo enddo !!end pixel by pixel double loop return end