*============================================================================== * This program calculates the completeness of observations over * the 2QZ survey strips and writes out a mask file in 1'x1' pixels. * This version uses the new sectors data in the catalogue. * * to compile: * f77 -o make_mask make_mask.f 2dfcatread.f check_raholes.f test_lims.f -fast * program compmask * * written by S.M. Croom 3/6/99, last updated: SMC 29/04/03 * * version 2.0: SMC 01/08/01. Now use sector data to make completeness mask. * version 2.1: SMC 30/10/01. Bug fix to read 3rd dummy line for limits file. * version 2.2: SMC 03/05/02. string length for nsecname increased after * addition of MGC data. * version 2.3: SMC 17/05/02. Add option of producing a mask based on the * number of good quality observations rather * than just the number of observed objects. * version 2.4: SMC 19/02/03 Prints total and effective area to screen. * version 2.5: SMC 29/04/03 good quality completeness is based on the * ratio n_good/n_observe, not n_good/n_cat. * Also include an estimate of N_obs/N_est for * each sector. * *############################################################################### * # * LICENCE # * # * This program is free software; you can redistribute it and/or # * modify it under the terms of the GNU General Public License # * as published by the Free Software Foundation; either version 2 # * of the License, or (at your option) any later version. # * # * This program is distributed in the hope that it will be useful, # * but WITHOUT ANY WARRANTY; without even the implied warranty of # * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # * GNU General Public License for more details. # * # * You should have received a copy of the GNU General Public License # * along with this program; if not, write to the Free Software # * Foundation, Inc., 59 Temple Place - Suite 330, # * Boston, MA 02111-1307, USA. # * # *############################################################################### * * Definition of variables: implicit none character*80 dumchar,catfile,secfile,holfile,limfile,mfile,mfileb character*80 mfilegood,mfilebgood,mfileratio,mfilebratio character*30 secname,nsecname*51,strip*1 integer maxsec,maxpos parameter(maxsec=3000,maxpos=12) character*30 secnames(maxsec) integer sec(maxsec),not(maxsec),nran,nsectors,nlims,ncat,qlim integer nobj(maxsec),nobserve(maxsec),nh,inhole,ngood(maxsec) integer i,mobj,j,j1,j2,ranum,decnum,nsec,nnot,iargc,narg integer nobstot,insector,secnum parameter(mobj=30000,ranum=4510,decnum=320) integer jsec(mobj),inobs(mobj) real bjmag(mobj),nnratio(maxsec),nest(maxsec) real maskratio(ranum,decnum),mratio,areagood real area(maxsec),secra(maxsec,maxpos),secdec(maxsec,maxpos) real seccdec(maxsec,maxpos),notcdec(maxsec,maxpos) real notra(maxsec,maxpos),notdec(maxsec,maxpos),frac(maxsec) real dtor,diff,maxdiff,ra,dec,pi,cdec,ratiog,areaobs,areatot real ratio,mask(ranum,decnum),fracgood(maxsec),maskg(ranum,decnum) real ramin,ramax,decmin,decmax,decstep,rastep,decav,version double precision raxhole,rayhole,dra,ddec double precision ramins(15),ramaxs(15),decmins(15),decmaxs(15) parameter(dtor=0.017453292,pi=3.141592654) parameter(version=2.5) logical withinfld * * catalogue parameters: include '2dfcatformat.h' * * Common and data blocks: common /raholes/ raxhole(4,400),rayhole(4,400),nh common /fldlims/ ramins,ramaxs,decmins,decmaxs,nlims * maximum radius (in degrees) data maxdiff /1.0/ * quality limit for good quality objects: data qlim /20/ ! quality 1 objects only * data qlim /40/ ! all observed objects *------------------------------------------------------------------------------ * write(*,1) version 1 format('make_mask (v',f3.1,')') * * set some parameters: maxdiff=(maxdiff*dtor)**2.0 * * iitialize some arrays: do i=1,maxsec nobj(i)=0 nobserve(i)=0 ngood(i)=0 nest(i)=0.0 nnratio(i)=0.0 end do * * get command line input: narg=iargc() if (narg.ne.4) then write(*,*) &'usage: make_mask ' call exit(2) end if call getarg(1,catfile) call getarg(2,secfile) call getarg(3,holfile) call getarg(4,limfile) * * check for SGP or NGP: do i=1,80-2 if (catfile(i:i+2).eq.'sgp') then strip='S' write(*,*) 'SGP strip' goto 50 else if (catfile(i:i+2).eq.'ngp') then strip='N' write(*,*) 'NGP strip' goto 50 else write(*,*) 'ERROR: dont know which strip this is. You need' write(*,*) ' to have sgp or ngp in the catalogue file' write(*,*) ' name' call exit(2) end if end do 50 write(*,*) * * read sectors file: write(*,*) 'reading sectors file...' nsectors=0 open(20,file=secfile,status='old') do i=1,maxsec read(20,'(a30,a51)',err=9999,end=100) secname,nsecname nsec=0 do j=1,30 if (secname(j:j).eq.'_') nsec=nsec+1 end do nnot=0 do j=1,50 if (nsecname(j:j).eq.'_') nnot=nnot+1 end do sec(i)=nsec not(i)=nnot backspace(20) secnames(i)=secname read(20,*,err=9999,end=100) secname,nsecname,area(i),nran & ,(secra(i,j),secdec(i,j),j=1,nsec),(notra(i,j),notdec(i,j) & ,j=1,nnot) do j=1,nsec if (secra(i,j).gt.3.0*pi/2.0) secra(i,j)=secra(i,j)-2.0*pi seccdec(i,j)=cos(secdec(i,j)) end do do j=1,nnot if (notra(i,j).gt.3.0*pi/2.0) notra(i,j)=notra(i,j)-2.0*pi notcdec(i,j)=cos(notdec(i,j)) end do nsectors=nsectors+1 end do 100 close(20) write(*,*) 'number of sectors: ',nsectors write(*,*) * * read in holes: write(*,*) 'reading holes file...' nh=0 open(20,file=holfile,status='old',err=9999) read(20,*) dumchar read(20,*) dumchar do i=1,1000 read(20,*,err=9999,end=200) raxhole(1,i),rayhole(1,i),raxhole(2 & ,i),rayhole(2,i),raxhole(3,i),rayhole(3,i),raxhole(4,i) & ,rayhole(4,i) do j=1,4 if (raxhole(j,i).gt.3.0*pi/2.0) raxhole(j,i)=raxhole(j,i)-2 & .0*pi end do nh=nh+1 end do 200 close(20) write(*,*) 'number of holes: ',nh write(*,*) * * read in field limits: write(*,*) 'reading limits file...' nlims=0 open(20,file=limfile,status='old',err=9999) read(20,*) dumchar read(20,*) dumchar read(20,*) dumchar do i=1,1000 read(20,*,err=9999,end=300) ramins(i),ramaxs(i),decmins(i) & ,decmaxs(i) if (ramins(i).gt.3.0*pi/2.0) ramins(i)=ramins(i)-2.0*pi if (ramaxs(i).gt.3.0*pi/2.0) ramaxs(i)=ramaxs(i)-2.0*pi nlims=nlims+1 end do 300 close(20) write(*,*) 'number of field limits: ',nlims write(*,*) * * read catalogue file and find fraction of objects observed: write(*,*) 'reading catalogue file...' ncat=0 nobstot=0 open(20,file=catfile,status='old',err=9999) do i=1,1000000 call catread(20,line,ier) if (ier.eq.1) goto 500 insector=0 bjmag(i)=b inobs(i)=nobs do j=1,nsectors if (sector.eq.secnames(j)(1:25)) then jsec(i)=j insector=insector+1 nobj(j)=nobj(j)+1 if (nobs.ge.1) then nobstot=nobstot+1 nobserve(j)=nobserve(j)+1 if (zq1.lt.qlim.and.zq1.gt.0) then ngood(j)=ngood(j)+1 end if end if end if end do if (insector.ne.1) then write(*,*) 'WARNING: error matching sector for ',iauname,' ' & ,insector end if ncat=ncat+1 end do 500 close(20) write(*,*) 'number of objects read: ',ncat write(*,*) 'number of objects observed: ',nobstot write(*,*) * * next go through all the sectors and find the fraction of observed * objects. This bit has been altered (SMC: 29/04/03) to give * fracgood(i)=ngood(i)/nobserve(i) rather than fracgood(i)=ngood(i)/nobj(i) * so that the estimates of spectroscopic completeness are separate from * the coverage estimates: areaobs=0.0 areatot=0.0 areagood=0.0 do i=1,nsectors if (nobj(i).gt.0) then frac(i)=float(nobserve(i))/float(nobj(i)) else frac(i)=1.0 end if if (nobserve(i).gt.0) then fracgood(i)=float(ngood(i))/float(nobserve(i)) else fracgood(i)=1.0 end if areaobs=areaobs+area(i)*frac(i) areatot=areatot+area(i) areagood=areagood+area(i)*frac(i)*fracgood(i) end do write(*,20) areatot write(*,30) areaobs write(*,40) areagood 20 format('total area=',f8.3,' sq. deg.') 30 format('observed area=',f8.3,' sq. deg.') 40 format('area with good spectra=',f8.3,' sq. deg.') * * next sum over the f_s(b_J) values for each object to define the * values for N_obs/N_est: do i=1,ncat if (inobs(i).gt.0) then secnum=jsec(i) nest(secnum)=nest(secnum)+(1-exp(bjmag(i)-20.388)*(1 & -fracgood(secnum))**0.919) end if end do * now calculate N_obs/N_est: write(*,*) 'calculating N_obs/N_est...' do i=1,nsectors if (nest(i).gt.0.0) then nnratio(i)=float(nobserve(i))/nest(i) else nnratio(i)=1.0 end if end do * * next set up the RA, Dec grid: if (strip.eq.'N') then decmin=-0.04783743 decmax=0.04366513 ramin=2.57286209 ramax=3.88225610 decav=0.0 else if (strip.eq.'S') then decmin=-0.567232 decmax=-0.479966 ramin=-0.613252 ramax=0.853235 decav=-30.0 end if * 1' steps: decstep=1.0/60.0*dtor rastep=1.0/60.0*dtor/cos(decav*dtor) decmin=decmin+decstep/2.0 ramin=ramin+rastep/2.0 * * open outout files: mfile='new_cov.mask' mfileb='new_cov.mask_bin' mfilegood='new_comp.mask' mfilebgood='new_comp.mask_bin' mfileratio='new_nnratio.mask' mfilebratio='new_nnratio.mask_bin' open(22,file=mfile) open(23,file=mfileb,form='unformatted') open(24,file=mfilegood) open(25,file=mfilebgood,form='unformatted') open(26,file=mfileratio) open(27,file=mfilebratio,form='unformatted') write(22,*) ramin/dtor write(22,*) decmin/dtor write(22,*) ranum write(22,*) decnum write(22,*) rastep/dtor write(22,*) decstep/dtor write(23) ramin/dtor write(23) decmin/dtor write(23) ranum write(23) decnum write(23) rastep/dtor write(23) decstep/dtor write(24,*) ramin/dtor write(24,*) decmin/dtor write(24,*) ranum write(24,*) decnum write(24,*) rastep/dtor write(24,*) decstep/dtor write(25) ramin/dtor write(25) decmin/dtor write(25) ranum write(25) decnum write(25) rastep/dtor write(25) decstep/dtor write(26,*) ramin/dtor write(26,*) decmin/dtor write(26,*) ranum write(26,*) decnum write(26,*) rastep/dtor write(26,*) decstep/dtor write(27) ramin/dtor write(27) decmin/dtor write(27) ranum write(27) decnum write(27) rastep/dtor write(27) decstep/dtor * * loop over grid: write(*,*) 'making mask...' do j1=1,ranum ra=ramin+rastep*float(j1-1) if (mod(j1,500).eq.0) write(*,10) j1,ranum 10 format(i5,' of ',i5,' RA steps done') do j2=1,decnum dec=decmin+decstep*float(j2-1) cdec=cos(dec) ratio=0.0 ratiog=0.0 dra=dble(ra) ddec=dble(dec) call test_lims(dra,ddec,withinfld) if (.not.withinfld) goto 790 call check_raholes(dra,ddec,inhole) if (inhole.eq.0) goto 790 * find out which sector the pixel is in: do i=1,nsectors do j=1,sec(i) diff=(ra-secra(i,j))**2.0*cdec*seccdec(i,j)+(dec & -secdec(i,j))**2.0 if (diff.gt.maxdiff) goto 600 end do do j=1,not(i) diff=(ra-notra(i,j))**2.0*cdec*notcdec(i,j)+(dec & -notdec(i,j))**2.0 if (diff.le.maxdiff) goto 600 end do ratio = frac(i) ratiog= fracgood(i) mratio = nnratio(i) goto 790 600 end do 790 mask(j1,j2)=ratio maskg(j1,j2)=ratiog maskratio(j1,j2)=mratio end do end do write(22,*) mask write(23) mask write(24,*) maskg write(25) maskg write(26,*) maskratio write(27) maskratio close(22) close(23) close(24) close(25) close(26) close(27) write(*,*) write(*,*) 'mask written to: ',mfile(1:30) write(*,*) 'binary version written to: ',mfileb(1:30) write(*,*) 'good obs mask written to: ',mfilegood(1:30) write(*,*) 'good obs binary version written to: ',mfilebgood(1:30) write(*,*) 'N_obs/N_est mask written to: ',mfileratio(1:30) write(*,*) 'N_obs/N_est binary version written to: ' & ,mfilebratio(1:30) stop 9999 stop 'file read error' end *==============================================================================