C This is a template program for reading the 2dFgg database object fits files 
C For the filename provided on the command line, the DSS image is read into   
C an array, and APM and DSS data are read from primary header keywords        
C The number of spectral extensions ( observations ) are determined by        
C stepping though the FITS file. The user is prompted to select one of the    
C spectral extensions to read the data from. The object spectra, variance and 
C sky are loaded into arrays and simple operations performed on the data.     
C Data related to the spectra are read from keywords in the header of the     
C chosen spectral extension ( redshift, quality, observed ra & dec. The file  
C is closed and allocated memory released. Written by Ian Price, June 1999.    
C                                                                             
C To compile include files "fitsio.h" and "longnam.h", and the   
C library libcfitsio.a, so that the compile command will look something like
C
C f77 -o rdbfits 2dF_rdbfits.f -L/usr/local/lib -L/usr/local -lsla -lcfitsio -lsocket -lnsl
C
C Translated from the C version by  Karl Glazebrook, June 1999  
C Minor changes: more long names,main stuff in subroutine, mean stats, etc.   
C   Karl Glazebrook, June 1999.                                               
C Minor updates for new fields: Carole Jackson, Apr 2001

      program tdfgrs_example
 
      implicit none
      
      integer BUFFSIZE

      parameter( BUFFSIZE=1024 )    ! buffer size when reading from fits file
      
C    *** Critical hard-wired parameters *** (though won't change for 2dF)      
      
      integer DSS_SIZE              ! Size of a DSS image (pixels)
      parameter (DSS_SIZE=49) 
      
      integer SPEC_SIZE            ! Size of a spectrum (pixels)
      parameter (SPEC_SIZE=1024)


      integer status,header_unit
      real  spec_data(SPEC_SIZE),vari_data(SPEC_SIZE)
      real  sky_data(SPEC_SIZE)
      real  dss_image(DSS_SIZE,DSS_SIZE)
      real  sp_nullval,bjsel, mean
      integer funit
      character*256 filename
      character*80 comment
      integer*2 dss_buff(BUFFSIZE),im_nullval
      integer im_naxis1,im_naxis2,sp_naxis1,sp_naxis2,nbuffer,fpixel
      integer npixel,serial
      integer igal,num_spec,quality
      real*8 ra,dec,cdelt1,crpix1,crval1,obsra,obsdec
      real plate_scale,z,z_helio, snr, eta
      logical anynull
      
      integer i,ii,jj
      
      ! Internal functions
      
      real*8 calc_mean
      
      integer nargs
      integer iargc  ! Unix function

      nargs = iargc()
      if (nargs.ne.1) then
   	  print *,' Usage : template 2dF_database_fits_file'
   	  call exit(1) ! Unix function
      endif
      
      call getarg(1,filename)

C    FITS status flag is 0 until an error occurs. If not 0, the fits
C    routine will  be skipped over. 				

      status = 0

C    open the database object FITS file for reading

      call ftnopn(funit,filename,0, status)

      if (status.ne.0) then
   	  print *,' File NOT Found - ',filename
      endif

C     Read APM and DSS image related keywords from primary FITS header      
C
C     Available keywords are :- 					    
C      	     								    
C     FITS descriptive keywords 					    
C     SIMPLE   BITPIX	NAXIS	 NAXIS1   NAXIS2   EXTEND   BSCALE   BZERO  
C      	     								    
C     2dFGRS/APM descriptive keywords					    
C     SEQNUM   NAME	IMAGE	 RA	  DEC	   EQUINOX  BJSEL    PROB   
C     PARK     PARMU	IGAL	 JON	  ORIENT   ECCENT   AREA     X_BJ     
C     Y_BJ     DX       DY	 BJG	  RMAG     PMAG     FMAG     SMAG
C     IFIELD   IGFIELD  REGION   OBJEQNX  OBJRA    OBJDEC   PLTSCALE XPIXELSZ
C     YPIXELSZ OBJPLTX  OBJPLTY  DATAMAX  DATAMIN  BJSELOLD BJG_OLD
 
C    get the DSS image size keywords

      call ftgkyj(funit, 'NAXIS1', im_naxis1, comment, status)
      call ftgkyj(funit, 'NAXIS2', im_naxis2, comment, status)
      
C read the DSS(x,y) image via a buffer - probably overkill for the small 2dFgg files 
C Note unlike the C example we can not dynamically allocate memory for the array with
C malloc so we use a fixed size array dss_image(DSS_SIZE,DSS_SIZE). (The alternative
C would be to call an OS-dependent subroutine).

      im_nullval = 0.0 ! do not perform checking for null data values

      npixel = im_naxis1*im_naxis2
      fpixel = 1
      
      do while (npixel.gt.0)
      
          nbuffer = npixel
    	  if (nbuffer.gt.BUFFSIZE) nbuffer = BUFFSIZE

C   	 read in a buffer full of data from the image array
 
	  call ftgpvi(funit,0,fpixel,nbuffer,im_nullval, dss_buff,
     :	      anynull,status)

C   	 copy the buffer - useful for type casts ( ie. read float image, and store
C  	 as doubles 							      

          do i=1,nbuffer
	      jj = int( (fpixel+i-2)/DSS_SIZE )
	      ii = fpixel+i-2 - jj * DSS_SIZE
	      dss_image(ii+1,jj+1) = dss_buff(i)
   	  end do
   	  npixel = npixel - nbuffer
   	  fpixel = fpixel + nbuffer
      
      end do
      
      mean = calc_mean(dss_image,   im_naxis1*im_naxis2)
      print *,'Mean of DSS image = ',mean,' counts'
      
C    read APM keywords - BJSEL RA DEC SEQNUM IGAL 

      call ftgkyd(funit, 'RA', ra, comment, status)
      call ftgkyd(funit, 'DEC', dec, comment, status)
      call ftgkye(funit, 'BJSEL', bjsel, comment, status)
      call ftgkyj(funit, 'SEQNUM', serial, comment, status)
      call ftgkyj(funit, 'IGAL', igal, comment, status)
      call ftgkye(funit, 'PLTSCALE', plate_scale, comment, status)
      
      if (status.eq.0) then
         print *,' '
	 print *,'APM parameters:'
         print *,' '
	 print *,'BJSEL = ',bjsel
   	 print *,'RA = ',ra
	 print *,'DEC = ',dec
	 print *,'SEQNUM = ',serial
	 print *,'IGAL = ',igal
	 print *,'PLTSCL = ',plate_scale
	 print *,' '
      else
   	  print *,
     :	    'Failed to read data and keywords from the primary image'
      endif

C    Finished with the primary header and DSS image - move on to the spectra

      num_spec = 0
      do while (status.eq.0)
         call ftmrhd(funit,1,header_unit,status)
	 if (status.eq.0) num_spec = num_spec +1
      end do
      
      print *,num_spec,' spectral extensions found'

C    reset status to zero, and select spectral extension interegate
     
      if (num_spec.ne.0) then
   	  i = 0
   	  do while (i.lt.1 .or. i.gt.num_spec)
   	      print *,
     :	        'Which spectrum do you want to read from file (1..',
     :		num_spec,')'
   	      read(5,*) i
   	  end do

C        resetting the status flag 

   	  status = 0

C   	 move to the selected extension

	  call ftmahd(funit,i+1, header_unit,status)
 
   	  if (status.ne.0) print *,
     :	          'Failed to move to the chosen extension'
	  
C   	 read the spectra/variance/sky image dimension keywords 

          call ftgkyj(funit, 'NAXIS1', sp_naxis1, comment, status)
          call ftgkyj(funit, 'NAXIS2', sp_naxis2, comment, status)
      
C   	check that all components exist - could also check sp_naxis1==1024 

   	  if (sp_naxis2 .ne. 3) then
   	      status = 1
   	      print *,'fits file has invalid structure'
   	  endif
 
   	  print *,'Dimensions of spectrum array = ', sp_naxis1
	  
C   	read the wavelength scale keywords

C  	wavelength defined as CRVAL1 + CDELT1*(pixel - CRPIX1), where pixel
C       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)
	  
C   	read the spectral keywords Z Z_HELIO OBSRA OBSDEC QUALITY 
	  
          call ftgkye(funit, 'Z', z, comment, status)
          call ftgkye(funit, 'Z_HELIO', z_helio, comment, status)
          call ftgkyd(funit, 'OBSRA', obsra, comment, status)
          call ftgkyd(funit, 'OBSDEC', obsdec, comment, status)
          call ftgkyj(funit, 'QUALITY', quality, comment, status)

      	  if (status.eq.0) then
      	     print *,' '
	     print *,'Spectrum parameters:'
      	     print *,' '
	     print *,'Z = ',z
	     print *,'Z_HELIO = ',z_helio
   	     print *,'OBSRA = ',obsra
	     print *,'OBSDEC = ',obsdec
	     print *,'QUALITY = ',quality
	     print *,' '
      	  else
   	      print *,'Failed to read spectral keywords'
      	  endif

c       read the optional keywords S/N and eta

          call ftgkye(funit, 'SNR', snr, comment, status)
      	  if (status.eq.0) then
      	     print *,'SNR = ', snr
          endif
          call ftgkye(funit, 'ETA', eta, comment, status)
      	  if (status.eq.0) then
      	     print *,'ETA = ', eta
          endif


C   	load the data without buffering, and replace null data values with
C   	-1.0e20								

   	  sp_nullval = -1.0e20

	  call ftgpve(funit,0,1,sp_naxis1,sp_nullval, 
     :	       spec_data,anynull,status)
     	  call ftgpve(funit,0,sp_naxis1+1,sp_naxis1,sp_nullval, 
     :	       vari_data,anynull,status)
	  call ftgpve(funit,0,2*sp_naxis1+1,sp_naxis1,sp_nullval, 
     :	       sky_data,anynull,status)
	  
C   	 Do some stuff with the data 
 
   	  if (status.eq.0) then
   	     call do_stuff_to_data (spec_data, vari_data, 
     :	          sky_data, sp_naxis1, crval1, crpix1, cdelt1)
   	  endif
 
       end if

C    reset the status flag before closing to ensure it closes 

      status = 0
 
      call ftclos(funit,status)

      end
      
C    Trvial subroutine to calculate the mean of an array   
   
      real*8 function calc_mean (a, n)
      
      implicit none
      
      integer n
      real a(n)
      
      integer i
      real sum
      
      sum = 0.0
      do i=1,n
         sum = sum + a(i)
      end do
      
      calc_mean = sum/dble(n)
      
      end
      	
      subroutine do_stuff_to_data (mydata, var, sky, nelem, 
     :               crval1, crpix1, cdelt1)

C     Do some more stuff with the data e.g. 
C     perform operations on the data - 
C     take square root of variance			    
C     turn sky counts into flux in arbitrary units	    
C     Work out mean of good data 
 
      implicit none
      integer nelem
      real mydata(nelem), var(nelem), sky(nelem)
      real*8  crval1, crpix1, cdelt1

      integer i,n
      real lambda
      real*8 sum

      n=0
      sum=0.0
      
      do i=1,nelem
      
C       Wavelength of each pixel 
 
          lambda = crval1 + cdelt1*( real(i) - crpix1 )
          if (i.eq.1) print *,'Wavelength min = ', lambda,' A'
          if (i.eq.nelem) print *,'Wavelength max = ', lambda,' A'
 
          if ( var(i).gt.0.0 ) then  ! Good data has positive variance
              var(i) = sqrt( var(i) )
              sum = sum + mydata(i)
              n = n+1
          else 
              var(i) = -1.0
          end if
        end do
 
        print *, 'Mean of good data = ', sum/dble(n),' units'

      end