* SPECFITS * * Lists in ASCII format the object spectrum, variance array and sky spectrum * from a specified 2dFGRS database FITS file and extension number. * * 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 * all extensions are used; if XTN>0 then only the specified extension * is used. If NAME is a list of FITS files then if XTN=0 all extensions * are used; if XTN>0 it is just computed for the extensions specified * on each line of the file. * * Matthew Colless * * Original version: 25 June 2001 * 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,nargs,idot,nfiles,ioerr,xtn,extn(MAXFILES),quality, : abemma real spec_data(SPEC_SIZE),vari_data(SPEC_SIZE),sky_data(SPEC_SIZE), : sp_nullval,bjsel,z real*8 cdelt1,crpix1,crval1 character filename*20,ff*20,fff*50,fitsfile(MAXFILES)*20,comment*80, : name*10,ftype*10,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 use all extensions (0) or' print *,' 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 * 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. fff='/2dFGRS_Database/'//ff ! directory containing the FITS files call ftnopn(funit,fff,0,status) if (status.ne.0) then print *,'ERROR: cannot open FITS file named ',fff 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 * * DSS image descriptive keywords * X_BJ Y_BJ DX DY BJG RMAG PMAG FMAG * SMAG REDMAG IFIELD IGFIELD REGION OBJEQNX OBJRA OBJDE* * PLTSCALE XPIXELSZ YPIXELSZ OBJPLTX OBJPLTY DATAMAX DATAMIN * 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 * 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 end if * Check that the object spectrum (1), variance array (2) and * sky spectrum (3) components exist and have sizes 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 ftgkyj(funit,'ABEMMA',abemma,comment,status) if (status.ne.0) then print *,'ERROR: failed to get spectral keywords' goto 7777 end if * Extract the object spectrum, variance array and sky spectrum. * 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 call ftgpve(funit,0,2*sp_naxis1+1,sp_naxis1,sp_nullval, : sky_data,anynull,status) if (status.ne.0) then print *,'ERROR: failed to get sky spectrum' goto 7777 end if * Print these in ASCII format to a text file. call printspec(fff,ff,seqnum,i,name,crval1,crpix1,cdelt1, : spec_data,vari_data,sky_data,sp_naxis1,bjsel, : z,quality,abemma) * 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 * PRINTSPEC * * Print the object spectrum, variance array and sky spectrum to * an ASCII text file. subroutine printspec(fff,ff,seqnum,xtn,name,crval1,crpix1,cdelt1, : spec,var,sky,nelem,bjsel,z,quality,abemma) implicit none character fff*50,ff*20,off*6,name*10,outfile*80 integer seqnum,xtn,nelem,i,ioerr,in,quality,abemma,lnb real spec(nelem),var(nelem),sky(nelem),lambda,outspec,outvar,outsky, : bjsel,z real*8 crval1,crpix1,cdelt1 * Open the ASCII text file for output. in = index(ff,'.fits') off = ff(in-6:in-1) write(outfile,'(a6,''_'',i1,''.txt'')') off,xtn open(50,file=outfile,status='unknown',iostat=ioerr) if (ioerr.ne.0) then print *,'ERROR: failed to open output file',outfile return end if * Write the object header information. write (50,'(''# FILE: '',a)') fff write (50,'(''#'')') write (50,'(''# SEQNUM = '',i6.6,'' EXTNUM = '',i1, : '' NAME = '',a10)') seqnum,xtn,name write (50,'(''# BJSEL = '',f5.2,'' Z ='',f9.6,'' QUALITY = '',i1, : '' ABEMMA = '',i1)') bjsel,z,quality,abemma write (50,'(''#'')') write (50,'(''# Pixel Wavlength Object Variance Sky'')') write (50,'(''# Number Angstroms Counts Counts Counts'')') * Write out the pixel number, wavelength, object spectrum, variance array * and sky spectrum. do i = 1,nelem lambda = crval1+cdelt1*(real(i)-crpix1) outspec = spec(i) if (spec(i).lt.0) outspec = 0 outvar = var(i) if (var(i).lt.0) outvar = 0 outsky = sky(i) if (sky(i).lt.0) outsky = 0 * write(50,*) i,lambda,spec(i),var(i),sky(i) write(50,'(i8,x,4f10.2)') i,lambda,outspec,outvar,outsky end do * Close the outpute file close(50) * Notify progress to the user lnb = index(outfile,' ')-1 write (*,'(''SEQNUM='',i6.6,'' EXTNUM='',i1,'' NAME='',a, : '' --> '',a)') seqnum,xtn,name,outfile(1:lnb) end