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 night the night of the frame being indexed C iframe the number of the frame being indexed 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 airplane(night, iframe, i) 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,night,iframe,frame,affine,q) implicit none integer*8 TOP_KMAX, kk real*4 frame(id,jd), yoffset(5), xoffset(5) real*4 node(TOP_KMAX) integer*8 nodeid(TOP_KMAX) integer*2 hits(TOP_KMAX) real*8 x(4), y(4), ra(4), dec(4), euler(4), q(4), affine(6) real*8 rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq, czsum 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 airplane, iairplane, night, iframe, 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, zz 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/ C r0=1194.5D0 !!old numbers... r0=1198.5D0 !!goose it up to eliminate bright band at edge... C numbers for camera 1,2,3 from SMEI Point v.1.0 C data yoffset /0.,0.,0.,0.,55./ C data yref /0.,61.612891D0,61.579259D0,59.0449220,0./ C data x00 /630.D0,631.96653D0,625.78564D0,627.96738D0,633.D0/ C data y00 /1291.38D0,1241.6218D0,1241.4524D0,1237.6987D0,1294.3D0/ C r0=y00(ic+1)-yref(ic+1) x0(ic+1)=x00(ic+1)-xoffset(ic+1) y0(ic+1)=y00(ic+1)-yoffset(ic+1) if(mode.eq.1.and.ic.gt.0)then ymin=1 ymax=jd xmin=22 xmax=1263 !!(1264-1) to remove the last half covered pixel xoff=0 else if(mode.eq.2.and.ic.gt.0)then ymin=1 ymax=jd xmin=11 xmax=632 xoff=0 else if(mode.eq.4.and.ic.gt.0)then ymin=1 ymax=jd xmin=6 xmax=316 xoff=0 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 !!less than 1.5 deg x 18.0226 pix/deg if( dr.lt.30.D0.and.dr.gt.-26.D0 )then if( (dabs( (i-x0(ic+1))/(j-y0(ic+1)) ) ).lt.0.57735D0 )then !!less than tan 30 response=response*(rccd/r0) cc loop for each vertex of a pixel czsum=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))); thetax = dasind( (x(k) - x0(ic+1)) / rccd ); c thetax = thetax !!/ 0.999744D0 !!atmosphere adjustment !!quadratic coefficiants for the formula: rq=aqTY^2+bqTY+cq !!from July 11, 2001 memo(ignoring 3rd and 4th order terms) cc thetaxsq = thetax*thetax cc aq = 0.1164D0 - (10.D0**-5)*thetaxsq ) cc bq = ( 3.6D0*(10.D0**-5)*thetaxsq ) - 18.0266D0 cc cq = ( 0.000936D0*thetaxsq ) - ( 1.044D0*(10.D0**-7)*thetaxsq*thetaxsq ) rq = (rccd - r0) !!/1.013D0 !!y-plate scale adjustment thetay = rq / -18.0226D0 cc thetay = ( -bq - sqrt((bq*bq)-(4.D0*aq*(cq-rq))) ) / ( 2.D0*aq ) 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 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