C+ C NAME: C averagenight_allsky.f C PURPOSE: C This subroutine takes the three triangle arrays(node, nodeid, and hits) C and averages them and puts them into a two-dimensional data array(image) C which can be printed out and viewed as a picture using Surfer C CATEGORY: C Data processing C CALLING SEQUENCE: C call averagenight_allsky(node,nodeid,hits,TOP_KMAX,level,image1,image2,image3,npolar,spolar) C INPUTS: C node real*4 array of size TOP_KMAX storing the response value C corresponding to a triangle of any index C nodeid integer*8 array of size TOP_KMAX storing the name of a triangle C in format(uint64 in C++) useful for JHU Spatial Index routines C hits integer*2 array of size TOP_KMAX storing the number of frames contributing C to a triangles response(node array) C TOP_KMAX integer*8 the size of the three triangle arrays C level integer*2 the level of the spatial index to be used C image(1,2,3) integer*4 array of size (1200,1200), storing the final 3 equatorial images to be viewed C (s,n)polar integer*4 array of size(1600,1600) storing final polar(dec>50 or dec<50) images C OUTPUTS: C outputs the final image to printed out and viewed by the Surfer imaging program C CALLS: C TRIANGLECENTER C++ function returning the ra/dec of the center of a triangle C MODIFICATION HISTORY: C April 2003, Aaron Smith (UCSD/CASS) C Updated to include capability for all sky imaging in 3 equatorial plots, and two polar plots C- subroutine averagenight_allsky(node, nodeid, hits, index, level,imagecnt1, imagecnt2, imagecnt3, npolarimagecnt, spolarimagecnt) implicit none integer*8 index, i, j integer*2 level, hits(index), ibad integer*4 imagecnt1(1200,1200), imagecnt2(1200,1200), imagecnt3(1200,1200) integer*4 npolarimagecnt(1600,1600), spolarimagecnt(1600,1600), xarg, yarg integer*8 nodeid(index), imagenode real*4 node(index), image1(1200,1200), image2(1200,1200), image3(1200,1200) real*4 npolarimage(1600,1600), spolarimage(1600,1600), noderesponse, arg real*8 cra1, cra2, cra3, ra, dec, sdec cra1=60.D0 cra2=180.D0 cra3=300.D0 do i=1,1200 do j=1,1200 image1(i,j)=0.0 imagecnt1(i,j)=0 image2(i,j)=0.0 imagecnt2(i,j)=0 image3(i,j)=0.0 imagecnt3(i,j)=0 enddo enddo c do i=1,1600 do j=1,1600 npolarimage(i,j)=0.0 npolarimagecnt(i,j)=0 spolarimage(i,j)=0.0 spolarimagecnt(i,j)=0 enddo enddo do i=1,index c if(node(i).gt.0.0.and.hits(i).ge.1)then if(hits(i).ge.1)then !Hey, Aaron, why check BOTH hits and node? c print*,node(i),hits(i) imagenode = nodeid(i) noderesponse = node(i)/float(hits(i)) CALL TRIANGLECENTER(level, imagenode, ra, dec) c print*,ra,dec c c Make first RA/dec plot (0-120 in ra) c ibad=0 xarg = nint( -(ra-cra1)*10.D0 + 600.5D0 ) yarg = nint( dec*10.D0 + 600.5D0 ) if(xarg.lt.1.or.xarg.gt.1200.or.yarg.lt.1.or.yarg.gt.1200)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.355.00.and.noderesponse.lt.180.00) go to 976 if(arg.lt.5.00.and.noderesponse.gt.180.00) go to 976 endif image1(xarg,yarg) = image1(xarg,yarg) + noderesponse 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)*10.D0 + 600.5D0 ) yarg = nint( dec*10.D0 + 600.5D0 ) if(xarg.lt.1.or.xarg.gt.1200.or.yarg.lt.1.or.yarg.gt.1200)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.355.00.and.noderesponse.lt.180.00) go to 976 if(arg.lt.5.00.and.noderesponse.gt.180.00) go to 976 endif image2(xarg,yarg) = image2(xarg,yarg) + noderesponse 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)*10.D0 + 600.5D0 ) yarg = nint( dec*10.D0 + 600.5D0 ) if(xarg.lt.1.or.xarg.gt.1200.or.yarg.lt.1.or.yarg.gt.1200)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.355.00.and.noderesponse.lt.180.00) go to 976 if(arg.lt.5.00.and.noderesponse.gt.180.00) go to 976 endif image3(xarg,yarg) = image3(xarg,yarg) + noderesponse 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( 20.*(90.-dec)*cosd(-ra) + 800.5) yarg = nint( 20.*(90.-dec)*sind(-ra) + 800.5) if(xarg.lt.1.or.xarg.gt.1600.or.yarg.lt.1.or.yarg.gt.1600)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(npolarimagecnt(xarg,yarg).ge.1)then arg=npolarimage(xarg,yarg)/float(npolarimagecnt(xarg,yarg)) if(arg.gt.355.00.and.noderesponse.lt.180.00) go to 976 if(arg.lt.5.00.and.noderesponse.gt.180.00) go to 976 endif npolarimage(xarg,yarg) = npolarimage(xarg,yarg) + noderesponse 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( -20.*(90.+dec)*cosd(ra) + 800.5) yarg = nint( -20.*(90.+dec)*sind(ra) + 800.5) if(xarg.lt.1.or.xarg.gt.1600.or.yarg.lt.1.or.yarg.gt.1600)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(spolarimagecnt(xarg,yarg).ge.1)then arg=spolarimage(xarg,yarg)/float(spolarimagecnt(xarg,yarg)) if(arg.gt.355.00.and.noderesponse.lt.180.00) go to 976 if(arg.lt.5.00.and.noderesponse.gt.180.00) go to 976 endif spolarimage(xarg,yarg) = spolarimage(xarg,yarg) + noderesponse spolarimagecnt(xarg,yarg) = spolarimagecnt(xarg,yarg) + 1 endif endif 976 endif enddo do i=1,1200 do j=1,1200 if(imagecnt1(i,j).gt.0) then image1(i,j)=image1(i,j)/float(imagecnt1(i,j)) imagecnt1(i,j)=nint(image1(i,j)) endif if(imagecnt2(i,j).gt.0) then image2(i,j)=image2(i,j)/float(imagecnt2(i,j)) imagecnt2(i,j)=nint(image2(i,j)) endif if(imagecnt3(i,j).gt.0) then image3(i,j)=image3(i,j)/float(imagecnt3(i,j)) imagecnt3(i,j)=nint(image3(i,j)) endif enddo enddo do i=1,1600 do j=1,1600 if(npolarimagecnt(i,j).gt.0) then npolarimage(i,j)=npolarimage(i,j)/float(npolarimagecnt(i,j)) npolarimagecnt(i,j)=nint(npolarimage(i,j)) endif if(spolarimagecnt(i,j).gt.0) then spolarimage(i,j)=spolarimage(i,j)/float(spolarimagecnt(i,j)) spolarimagecnt(i,j)=nint(spolarimage(i,j)) endif enddo enddo c The "cross" at the poles do i=800,801 do j=796,805 npolarimagecnt(i,j)=10000 spolarimagecnt(i,j)=10000 enddo enddo do i=796,805 do j=800,801 npolarimagecnt(i,j)=10000 spolarimagecnt(i,j)=10000 enddo enddo return end