
      SUBROUTINE rddtg1(unit, ifile, filenm, ftype, ispec, itype, 
     &                  phabuf, ndata, nungc, qsysdt, llchan, ldchan, 
     &                  hlchan, hdchan, ierr)


      INCLUDE '../../inc/xspec.inc'

      INTEGER nungc
      REAL    phabuf(nungc,5)
      INTEGER unit, ndata
      INTEGER llchan, ldchan, hlchan, hdchan
      INTEGER ierr, itype, ftype, ispec, ifile
      CHARACTER filenm*(*)
      LOGICAL qsysdt

c kaa 1/24/95
c	XSPEC subroutine called by RDDTOG to read from a spectral data file 
c	in OGIP FITS format. The file is assumed to be opened and stationed 
c       at the correct extension.
c
c	unit		i        i: The I/O unit for the PHA file
c       ifile           i        i: The dataset number
c       filenm          c        i: The dataset filename
c       ftype           i        i: The file type (I or II)
c       ispec           i        i: If type II then set to spectrum to read
c       itype           i        i: The type of data,
c					0 = counts,
c					1 = counts/sec,
c                                       2 = counts/sec/area
c       phabuf		r        r: The data, the first row of which contains
c                                   the data, the second contains the variance,
c                                   the third the systematics values, the
c                                   fourth the backscal values and the fifth
c                                   the areascal values.
c       ndata           i        r: The actual number of data points read
c       nungc           i        i: Max number of ungrouped channels
c       qsysdt          l        r: there was a systematics packages.
c       llchan          i        r: Lowest legal data channel
c       ldchan          i        r: Lowest actual data channel
c       hlchan          i        r: Highest legal data channel
c       hdchan          i        r: Highest actual data channel
c       ierr            i        r: Error flag (0= no error)

      REAL syserr, junk, backscal, areascal
      INTEGER index, filestat, ibin, ijunk, nsize
      CHARACTER comment*80, context*80, wrtstr*255
      LOGICAL qanyf, qpoiss

      INTEGER lenact
      EXTERNAL lenact

      filestat = 0

      context = 'RDDTG1 : FITSIO error'

c  Find the number of stored data points. For the type 2 files we need to
c  use different methods depending on whether we have variable or constant
c  length vectors.

      IF ( ftype .EQ. 1 ) THEN

         CALL ftgkyj(unit, 'NAXIS2', ndata, comment, filestat)
         context = 'RDDTG1: Failed to read NAXIS2'
         IF (filestat .NE. 0) GOTO 100

      ELSEIF ( ftype .EQ. 2 ) THEN

         IF ( itype .EQ. 0 ) THEN
            CALL ftgcno(unit, .FALSE., 'COUNTS', index, filestat)
            context = 'RDDTG1: Failed to find COUNTS column'
            IF (filestat .NE. 0) GOTO 100
         ELSEIF ( itype .EQ. 1 ) THEN
            CALL ftgcno(unit, .FALSE., 'RATE', index, filestat)
            context = 'RDDTG1: Failed to find RATE column'
            IF (filestat .NE. 0) GOTO 100
         ENDIF
         CALL ftgdes(unit, index, ispec, ndata, ijunk, filestat)
         IF ( filestat .NE. 0 ) THEN
            filestat = 0
            CALL ftgtdm(unit, index, 1, ijunk, ndata, filestat)
            context = 
     &        'RDDTG1: Failed to find size of COUNTS or RATE column'
            IF (filestat .NE. 0) GOTO 100
         ENDIF

      ENDIF

c  Check whether the TLMIN1 and TLMAX1 keywords are set. If so then
c  read the start and end channels in the data.

      CALL ftgkyj(unit, 'TLMIN1', llchan, comment, filestat)
      IF (filestat .NE. 0) THEN
         filestat = 0
         llchan = -1
      ENDIF
      CALL ftgkyj(unit, 'TLMAX1', hlchan, comment, filestat)
      IF (filestat .NE. 0) THEN
         filestat = 0
         hlchan = -1
      ENDIF

      IF (llchan .NE. -1 .AND. hlchan .NE. -1) THEN
         IF ( ftype .EQ. 1 ) THEN
            CALL ftgcvj(unit, 1, 1, 1, 1, 0, ldchan, qanyf, filestat)
            CALL ftgcvj(unit, 1, ndata, 1, 1, 0, hdchan, qanyf, 
     &                  filestat)
         ELSEIF ( ftype .EQ. 2 ) THEN
            CALL ftgcno(unit, .FALSE., 'CHANNEL', index, filestat)
            IF (filestat .EQ. 0 .AND. index .NE. 0 ) THEN
               CALL ftgcvj(unit, index, ispec, 1, 1, 0, ldchan, qanyf, 
     &                     filestat)
               CALL ftgcvj(unit, index, ispec, ndata, 1, 0, hdchan, 
     &                     qanyf, filestat)
               context = 'RDDTG1: Failed to read CHANNEL data'
               IF (filestat .NE. 0) GOTO 100
            ENDIF
            filestat = 0
         ENDIF
      ENDIF

c  If counts data then extract into phabuf

      IF (itype .EQ. 0) THEN

         CALL ftgcno(unit, .FALSE., 'COUNTS', index, filestat)
         IF (filestat .NE. 0) THEN
            context = 'RDDTG1: Failed to find COUNTS column'
            GOTO 100
         ENDIF
         IF ( ftype .EQ. 1 ) THEN
            CALL ftgcve(unit, index, 1, 1, ndata, 1.e-32, phabuf, qanyf,
     &                  filestat)
         ELSEIF ( ftype .EQ. 2 ) THEN
            CALL ftgcve(unit, index, ispec, 1, ndata, 1.e-32, phabuf, 
     &                  qanyf, filestat)
         ENDIF
         context = 'RDDTG1: Failed to read COUNTS data'
         IF (filestat .NE. 0) GOTO 100

c  otherwise read the `pha per sec' data

      ELSEIF (itype .EQ. 1) THEN

         CALL ftgcno(unit, .FALSE., 'RATE', index, filestat)
         IF (filestat .NE. 0) THEN
            context = 'RDDTG1: Failed to find RATE column'
            GOTO 100
         ENDIF
         IF ( ftype .EQ. 1 ) THEN
            CALL ftgcve(unit, index, 1, 1, ndata, 1.e-32, phabuf,
     &                  qanyf, filestat)
         ELSEIF ( ftype .EQ. 2 ) THEN
            CALL ftgcve(unit, index, ispec, 1, ndata, 1.e-32, phabuf,
     &                  qanyf, filestat)
         ENDIF
         context = 'RDDTG1: Failed to read RATE data'
         IF (filestat .NE. 0) GOTO 100

      ENDIF

c  check for the Poisson error keyword

      CALL ftgkyl(unit, 'POISSERR', qpoiss, comment, filestat)
      IF (filestat .NE. 0) THEN
         CALL xwrite('POISSERR keyword not found, assuming FALSE', 
     &               10)
         qpoiss = .FALSE.
         filestat = 0
      ENDIF

c  if the Poisson error is set but the data was RATE then write a warning
c  message

      IF ( qpoiss .AND. itype .EQ. 1 ) THEN
         CALL xwrite('*Error*: POISSERR=T used with RATE column',10) 
         CALL xwrite('       Data as RATE requires errors in STAT_ERR',
     &               10)
      ENDIF

c  if the Poisson error keyword is set then errors are square roots
c  otherwise load them from the stat_err column

      IF (qpoiss) THEN

         DO ibin = 1, ndata
            phabuf(ibin,2) = SQRT(MAX(phabuf(ibin,1),0.))
         ENDDO

         CALL ftgcno(unit, .FALSE., 'STAT_ERR', index, filestat)
         IF ( filestat .EQ. 0 .AND. index .NE. 0 ) THEN
            WRITE(wrtstr,'(a,a,a)') 'Warning : ', 
     &        filenm(:lenact(filenm)), 
     &        ' has both POISSERR=T and STAT_ERR'
            CALL xwrite(wrtstr, 10)
            CALL xwrite('Using Poisson errors', 10)
         ENDIF
         filestat = 0

      ELSE

         CALL ftgcno(unit, .FALSE., 'STAT_ERR', index, filestat)

         IF (filestat .NE. 0) THEN

            filestat = 0
            CALL ftgkye(unit, 'STAT_ERR', phabuf(1,2), 
     &                  comment, filestat)
            context = 'RDDTG1: Failed to find STAT_ERR'
            IF ( filestat .NE. 0 ) GOTO 100
            DO ibin = 2, ndata
               phabuf(ibin,2) = phabuf(1,2)
            ENDDO

         ELSE

            IF ( ftype .EQ. 1 ) THEN
               CALL ftgcve(unit, index, 1, 1, ndata, 1.e-32,
     &                     phabuf(1,2), qanyf, filestat)
            ELSEIF ( ftype .EQ. 2 ) THEN
               CALL ftgcve(unit, index, ispec, 1, ndata, 1.e-32, 
     &                     phabuf(1,2), qanyf, filestat)
            ENDIF
            context = 'RDDTG1: Failed to read STAT_ERR data'
            IF (filestat .NE. 0) GOTO 100

            CALL ftgkye(unit, 'STAT_ERR', junk, comment, filestat)
            IF ( filestat .EQ. 0 ) THEN
               WRITE(wrtstr, '(a,i4,a)') 'Warning : file ', ifile,
     &              ' has both STAT_ERR keyword and column'
               CALL xwrite(wrtstr, 10)
               CALL xwrite('Using data from STAT_ERR column', 10)
            ENDIF
            filestat = 0

         ENDIF

      ENDIF

c  check for systematic errors

      qsysdt = .FALSE.

      CALL ftgcno(unit, .FALSE., 'SYS_ERR', index, filestat)

      IF ( filestat .NE. 0 ) THEN

         filestat = 0
         CALL ftgkye(unit, 'SYS_ERR', syserr, comment, filestat)
         IF (filestat .EQ. 0 .AND. syserr .NE. 0. ) THEN
            DO ibin = 1, ndata
               phabuf(ibin,3) = syserr
            ENDDO
            qsysdt = .TRUE.
         ENDIF
         filestat = 0

      ELSEIF ( filestat .EQ. 0 .AND. index .GT. 0 ) THEN

         IF ( ftype .EQ. 1 ) THEN
            CALL ftgcve(unit, index, 1, 1, ndata, 1.e-32, 
     &                  phabuf(1,3), qanyf, filestat)
         ELSEIF ( ftype .EQ. 2 ) THEN
            CALL ftgcve(unit, index, ispec, 1, ndata, 1.e-32, 
     &                  phabuf(1,3), qanyf, filestat)
         ENDIF
         context = 'RDDTG1: Failed to read SYS_ERR data'
         IF ( filestat .NE. 0 ) GOTO 100
         qsysdt = .TRUE.

         CALL ftgkye(unit, 'SYS_ERR', junk, comment, filestat)
         IF ( filestat .EQ. 0 ) THEN
            WRITE(wrtstr, '(a,i4,a)') 'Warning : file ', ifile,
     &           ' has both SYS_ERR keyword and column'
            CALL xwrite(wrtstr, 10)
            CALL xwrite('Using data from SYS_ERR column', 10)
         ENDIF
         filestat = 0

      ENDIF

c  read the BACKSCAL value(s). Note that the type II case is more
c  complicated now because we need to have the option of a vector
c  BACKSCAL on each row.

      CALL ftgcno(unit, .FALSE., 'BACKSCAL', index, filestat)

      IF ( filestat .NE. 0 ) THEN

         filestat = 0
         CALL ftgkye(unit, 'BACKSCAL', backscal, comment, filestat)
         IF (filestat .EQ. 0 .AND. backscal .NE. 0. ) THEN
            DO ibin = 1, ndata
               phabuf(ibin,4) = backscal
            ENDDO
         ELSE
            WRITE(wrtstr, '(a,i4,a)') 'Warning : file ', ifile,
     &           ' has no BACKSCAL keyword or column'
            CALL xwrite(wrtstr, 10)
            CALL xwrite('Assuming BACKSCAL = 1', 10)
            DO ibin = 1, ndata
               phabuf(ibin,4) = 1.0
            ENDDO
         ENDIF
         filestat = 0

      ELSEIF ( filestat .EQ. 0 .AND. index .GT. 0 ) THEN

         IF ( ftype .EQ. 1 ) THEN
            CALL ftgcve(unit, index, 1, 1, ndata, 1.e-32, 
     &                  phabuf(1,4), qanyf, filestat)
         ELSEIF ( ftype .EQ. 2 ) THEN
            CALL ftgdes(unit, index, ispec, nsize, ijunk, filestat)
            IF ( filestat .NE. 0 ) THEN
               filestat = 0
               CALL ftgtdm(unit, index, 1, ijunk, nsize, filestat)
            ENDIF
            IF ( nsize .GT. 1 ) THEN
               CALL ftgcve(unit, index, ispec, 1, ndata, 1.e-32, 
     &                     phabuf(1,4), qanyf, filestat)
            ELSE
               CALL ftgcve(unit, index, ispec, 1, 1, 1.e-32, 
     &                     backscal, qanyf, filestat)
               DO ibin = 1, ndata
                  phabuf(ibin,4) = backscal
               ENDDO
            ENDIF
         ENDIF
         context = 'RDDTG1: Failed to read BACKSCAL data'
         IF ( filestat .NE. 0 ) GOTO 100

         CALL ftgkye(unit, 'BACKSCAL', junk, comment, filestat)
         IF ( filestat .EQ. 0 ) THEN
            WRITE(wrtstr, '(a,i4,a)') 'Warning : file ', ifile,
     &           ' has both BACKSCAL keyword and column'
            CALL xwrite(wrtstr, 10)
            CALL xwrite('Using data from BACKSCAL column', 10)
         ENDIF
         filestat = 0

      ENDIF

c  read the AREASCAL value(s). Note that the type II case is more
c  complicated now because we need to have the option of a vector
c  AREASCAL on each row.

      CALL ftgcno(unit, .FALSE., 'AREASCAL', index, filestat)

      IF ( filestat .NE. 0 ) THEN

         filestat = 0
         CALL ftgkye(unit, 'AREASCAL', areascal, comment, filestat)
         IF (filestat .EQ. 0 .AND. areascal .NE. 0. ) THEN
            DO ibin = 1, ndata
               phabuf(ibin,5) = areascal
            ENDDO
         ELSE
            WRITE(wrtstr, '(a,i4,a)') 'Warning : file ', ifile,
     &           ' has no AREASCAL keyword or column'
            CALL xwrite(wrtstr, 10)
            CALL xwrite('Assuming AREASCAL = 1', 10)
            DO ibin = 1, ndata
               phabuf(ibin,5) = 1.0
            ENDDO
         ENDIF
         filestat = 0

      ELSEIF ( filestat .EQ. 0 .AND. index .GT. 0 ) THEN

         IF ( ftype .EQ. 1 ) THEN
            CALL ftgcve(unit, index, 1, 1, ndata, 1.e-32, 
     &                  phabuf(1,5), qanyf, filestat)
         ELSEIF ( ftype .EQ. 2 ) THEN
            CALL ftgdes(unit, index, ispec, nsize, ijunk, filestat)
            IF ( filestat .NE. 0 ) THEN
               filestat = 0
               CALL ftgtdm(unit, index, 1, ijunk, nsize, filestat)
            ENDIF
            IF ( nsize .GT. 1 ) THEN
               CALL ftgcve(unit, index, ispec, 1, ndata, 1.e-32, 
     &                     phabuf(1,5), qanyf, filestat)
            ELSE
               CALL ftgcve(unit, index, ispec, 1, 1, 1.e-32, 
     &                     areascal, qanyf, filestat)
               DO ibin = 1, ndata
                  phabuf(ibin,5) = areascal
               ENDDO
            ENDIF
         ENDIF
         context = 'RDDTG1: Failed to read AREASCAL data'
         IF ( filestat .NE. 0 ) GOTO 100

         CALL ftgkye(unit, 'AREASCAL', junk, comment, filestat)
         IF ( filestat .EQ. 0 ) THEN
            WRITE(wrtstr, '(a,i4,a)') 'Warning : file ', ifile,
     &           ' has both AREASCAL keyword and column'
            CALL xwrite(wrtstr, 10)
            CALL xwrite('Using data from AREASCAL column', 10)
         ENDIF
         filestat = 0

      ENDIF

      RETURN

c Check for any FITSIO errors

  100 CONTINUE
      IF (filestat .NE. 0) THEN
         wrtstr = filenm(:lenact(filenm))//' : '//context
         CALL xwrite(wrtstr, 10)
         WRITE(wrtstr,'(''   FITSIO = '', i4)') filestat
         CALL xwrite(wrtstr, 10)
      ENDIF

      RETURN

      END
