*     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
