* SNRFITS * * Computes the S/N ratio of the 2dFGRS spectra. The S/N ratio is * computed as both the mean and the median S/N per pixel over the * specified wavelength range (4000-7500A). * * Parameters are: * 1. NAME: Either the name of a 2dFGRS FITS file or the name of a * file with list of FITS filenames and extensions (if XTN>0). * The filenames have format TTT/D/SSSSSS.fits, where TTT is the * tile number, D is the last digit of the serial number and * SSSSSS is the serial number. The file format is A17 if XTN=0 * and A17XI1 if XTN>0. * 2. XTN: If NAME is the name of a single FITS file, then if XTN=0 * the S/N is computed for all extensions and if XTN>0 it is just * computed for the specified extension. If NAME is a list of FITS * files then if XTN=0 the S/N is computed for all extensions and * if XTN>0 it is just computed for the extensions specified on * each line of the file. * * Compile with: f77 -e -o snrfits snrfits.f -lcfitsio -lsocket * * Author: Matthew Colless * * Original version: 26 October 2000 - S/N of the first spectral extension * Updated: 26 March 2001 - S/N for all observed spectra * Updated: 26 April 2001 - S/N for specified extension or all extensions * Parameters. implicit none integer MAXFILES ! maximum number of FITS files to process parameter (MAXFILES=500000) integer BUFFSIZE ! buffer size for reading FITS file parameter (BUFFSIZE=1024) integer SPEC_SIZE ! size of a spectrum (pixels) parameter (SPEC_SIZE=1024) * Variables. integer i,ii,status,header_unit,funit,seqnum,sp_naxis1,sp_naxis2, : num_spec,quality,nargs,idot,nfiles,ioerr,xtn, : extn(MAXFILES) real spec_data(SPEC_SIZE),vari_data(SPEC_SIZE),sp_nullval,z, : snmean,bjsel real*8 cdelt1,crpix1,crval1 character filename*20,ff*50,fitsfile(MAXFILES)*20,comment*80, : name*10,ftype*10,spfile*15,option,line*80 logical anynull * Internal functions. integer iargc ! Unix function * Parse command line. nargs = iargc() ! Unix function if (nargs.ne.2) then print *,'ERROR: requires two arguments:' print *,'1. Either the name of a 2dFGRS FITS file or the name' print *,' of a file with list of FITS filenames.' print *,'2. Whether to compute the S/N for all extensions (0)' print *,' or just one extension (>0).' goto 9999 end if call getarg(1,filename) ! Unix subroutine call getarg(2,option) ! Unix subroutine read(option,'(i1)') xtn * Decide if file on command line is a single FITS file or a list * of FITS files. idot = index(filename,'.') ftype = filename(idot+1:) if (ftype.eq.'fit' .or. ftype.eq.'fits' .or. : ftype.eq.'FIT' .or. ftype.eq.'FITS') then nfiles = 1 fitsfile(nfiles) = filename extn(nfiles) = xtn else nfiles = 0 open(10,file=filename,status='old',iostat=ioerr) if (ioerr.ne.0) then print *,'ERROR: cannot open list file named ',filename goto 9999 end if do while (ioerr.eq.0) read(10,'(a)',iostat=ioerr) line if (ioerr.eq.0) then nfiles = nfiles+1 if (xtn.eq.0) then read(line,'(a)') fitsfile(nfiles) extn(nfiles) = 0 else read(line,'(a17,x,i1)') fitsfile(nfiles),extn(nfiles) end if end if end do end if write (*,'(''# SEQNUM NAME b_J NS NX Z'', : '' Q SPfile SNRave SNRmed'')') * Loop over all the FITS files to be processed. do ii=1,nfiles ff = fitsfile(ii) xtn = extn(ii) * FITS status flag is 0 until an error occurs; if not 0, the FITS * routine will be skipped over. status = 0 * Open the database object FITS file for reading. ff='/priv/magus11/twodfgg/Database/'//ff call ftnopn(funit,ff,0,status) if (status.ne.0) then print *,'ERROR: cannot open FITS file named ',ff goto 8888 end if * Read APM and DSS image related keywords from primary FITS header. * Available keywords are: * * FITS descriptive keywords * SIMPLE BITPIX NAXIS NAXIS1 NAXIS2 EXTEND BSCALE BZERO * * 2dFGRS/APM descriptive keywords * SEQNUM NAME IMAGE RA DEC EQUINOX BJSEL PROB * PARK PARMU IGAL JON ORIENT ECCENT AREA X_BJ * Y_BJ DX DY BJG RMAG PMAG FMAG SMAG * IFIELD IGFIELD REGION OBJEQNX OBJRA OBJDEC PLTSCALE XPIXELSZ * YPIXELSZ OBJPLTX OBJPLTY DATAMAX DATAMIN BJSELOLD BJG_OLD * Read APM keywords. call ftgkyj(funit,'SEQNUM',seqnum,comment,status) call ftgkys(funit,'NAME',name,comment,status) call ftgkye(funit,'BJSEL',bjsel,comment,status) if (status.ne.0) then print *,'ERROR: failed to read APM keywords from primary image' goto 8888 end if * Get the number of spectral extensions. call ftthdu(funit,num_spec,status) num_spec = num_spec-1 if (status.ne.0 .or. num_spec.eq.0) then * print *,'ERROR: no (readable) spectral extensions' goto 8888 end if * Loop over all spectral extensions. do i=1,num_spec * Check if this extension is to be processed if (xtn.eq.0 .or. xtn.eq.i) then * Write the object information. write (*,'(2x,i6,2x,a10,2x,f5.2,2x,i2$)') seqnum,name,bjsel,num_spec * Select the next spectral extension. call ftmahd(funit,i+1,header_unit,status) if (status.ne.0) then print *,'ERROR: failed to select spectral extension',i goto 7777 end if * Read the spectrum and variance image dimension keywords. call ftgkyj(funit,'NAXIS1',sp_naxis1,comment,status) call ftgkyj(funit,'NAXIS2',sp_naxis2,comment,status) if (status.ne.0) then print *,'ERROR: failed to read spectrum dimension keywords' goto 7777 else write (*,'(2x,i2,$)') i end if * Check that spectrum (1) and variance (2) components exist and * size of spectrum is as expected. if (sp_naxis2 .ne. 3) then status = 1 print *,'ERROR: FITS file has invalid structure' goto 7777 end if if (sp_naxis1.ne.1024) then status = 1 print *,'ERROR: spectral dimension is not 1024' goto 7777 end if * Read the wavelength scale keywords. Wavelength defined as * CRVAL1 + CDELT1*(pixel - CRPIX1), where pixel starts at 1, * and finishes at NAXIS1. call ftgkyd(funit,'CDELT1',cdelt1,comment,status) call ftgkyd(funit,'CRVAL1',crval1,comment,status) call ftgkyd(funit,'CRPIX1',crpix1,comment,status) if (status.ne.0) then print *,'ERROR: failed to get wavelength scale keywords' goto 7777 end if * Read the spectral keywords. call ftgkye(funit,'Z',z,comment,status) call ftgkyj(funit,'QUALITY',quality,comment,status) call ftgkys(funit,'SPFILE',spfile,comment,status) if (status.ne.0) then print *,'ERROR: failed to get spectral keywords' goto 7777 else write (*,'(2x,f6.4,2x,i1,2x,a15$)') z,quality,spfile end if * Extract the spectrum and variance arrays. Replace null data values with 0. sp_nullval = 0 call ftgpve(funit,0,1,sp_naxis1,sp_nullval, : spec_data,anynull,status) if (status.ne.0) then print *,'ERROR: failed to get object spectrum' goto 7777 end if call ftgpve(funit,0,sp_naxis1+1,sp_naxis1,sp_nullval, : vari_data,anynull,status) if (status.ne.0) then print *,'ERROR: failed to get variance array' goto 7777 end if * Compute the S/N ratio. call setsnr(spec_data,vari_data,sp_naxis1,crval1,crpix1,cdelt1,snmean) write (*,'(2x,f6.1)') snmean * Reset the status flag before continuing. 7777 status = 0 * End loop over extensions to be processed. end if * End loop over spectral extensions. end do * Close FITS file 8888 status = 0 call ftclos(funit,status) end do 9999 status = 0 end * SETSNR * * Use the spectrum and the variance array to calculate the mean * S/N ratio over a standard wavelength range. subroutine setsnr(spec,var,nelem,crval1,crpix1,cdelt1,snmean) implicit none integer nelem,i,n real spec(nelem),var(nelem),lambda,minlam,maxlam, : sn(2048),snmean real*8 crval1,crpix1,cdelt1,sum minlam = 4000 maxlam = 7500 n = 0 sum = 0 do i = 1,nelem lambda = crval1+cdelt1*(real(i)-crpix1) if (lambda.ge.minlam .and. lambda.le.maxlam .and. : spec(i).gt.0 .and. var(i).gt.0) then n = n+1 sn(n) = spec(i)/sqrt(var(i)) sum = sum+sn(n) * write(7,*) n,lambda,sn(n) end if end do if (n.gt.0 .and. sum.gt.0) then snmean = sum/n else snmean = 0 end if end