SUBROUTINE wrtfak(phaobs, varobs, region, areasc, qbckfil)

      IMPLICIT NONE
 
      INCLUDE '../../inc/xspec.inc'
 
      REAL phaobs(*), varobs(*), region(*), areasc(*)
      LOGICAL qbckfil

C
 
c      wrtfak	burton 15 April 94
c 		XSPEC subroutine to write out the accumulated fake data
c		routines, using FITS format
c      3.0 - rewritten to write out data into FITS files
 
c  	phaobs   r    i: The count rate to be written out if counting statistics was set. 
c                        else an integer number of actual counts is written out.
c       varobs   r    i: The errors on the count rate if Poisson statistics were
c                        not assumed for the background.
c       region   r    i: The BACKSCAL values
c       areasc   r    i: The AREASCAL values
c       qbckfil  l    i: If true then write out the background file
c
c   7/21/97 - correctly write the keywords in the case of POISSERR = .FALSE.  
c   4/30/99 - change handling for pre-existing files: now dealt with before this function
c               is called. B. Dorman   
c   4/29/02 - added option to write out background files
c

      INTEGER icounts, ierror, ichannl, iqualty, igroup

      REAL grasp, backscal, areascal
 
      INTEGER istat, maxchn, iunit, ifile, i, bitpix, ichan, ibin
      INTEGER status, nrows, nfield, pcount, block, chan0, icol
      
      
      integer ridum
      character*20 rhdum
      character*3 rrchtyp
      real rrdum
 
      CHARACTER*11 cdate
      INTEGER yr,mo,day
      CHARACTER filenm*255, name*255
      CHARACTER contxt*72, extens*10
      CHARACTER cgroup*1
      CHARACTER*20 tform(6), ttype(6), tunit(6), extnam
      LOGICAL simple, extend, qqual, qgroup, qpoiss, qback, qarea
      LOGICAL qexist
 
      REAL dgtime, dgbtme, dgcrnm
      INTEGER dgndst, lenact, dgungc, dgnbnb, dgnbne, dgchn0, ftgsdt
      CHARACTER fgflnm*255, dgtlsc*20
      CHARACTER fgauxf*255, dginst*20, dggcrd*1
      LOGICAL dgqnot, dgqfst, dgdtps
      EXTERNAL lenact, fgflnm, fgauxf, dgtlsc, dgdtps
      EXTERNAL dgtime, dgbtme, dgcrnm, dgndst, dgungc, dginst
      EXTERNAL dgnbnb, dgnbne, dgchn0, dggcrd, dgqnot, ftgsdt
 
      status = 0
 
c Find largest number of channels in a data set.

      maxchn = 0
      DO ifile = 1, dgndst()
         maxchn = MAX(maxchn, dgungc(ifile))
      ENDDO

c Allocate memory for workspace arrays.

      CALL udmget( maxchn, 6, icounts, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     &     ' WARNING: Memory allocation failed in WRTFAK.', 2)      

      CALL udmget( maxchn, 6, ierror, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     &     ' WARNING: Memory allocation failed in WRTFAK.', 2)      

      CALL udmget( maxchn, 4, ichannl, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     &     ' WARNING: Memory allocation failed in WRTFAK.', 2)      

      CALL udmget( maxchn, 4, iqualty, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     &     ' WARNING: Memory allocation failed in WRTFAK.', 2)      

      CALL udmget( maxchn, 4, igroup, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     &     ' WARNING: Memory allocation failed in WRTFAK.', 2)      

c
c Define values that are constant
c
      bitpix = 8
      block = 2880
      status =  ftgsdt(day,mo,yr,status)
      write(cdate,'(i4,"-",i2.2,"-",i2.2)') yr,mo,day
      status = 0
c
c Loop around files to be written
c
      CALL getlun(iunit)
      DO ifile = 1, dgndst()

         IF ( qbckfil .AND. dgbtme(ifile) .EQ. 0. ) GOTO 300

         qpoiss = .TRUE.
         IF ( dgbtme(ifile) .GT. 0. .AND. .NOT.dgdtps(ifile,'b') ) 
     &      qpoiss = .FALSE.

c Create new file. fakeit has already checked whether the user
c wishes to overwrite this file, so if it does already exist, delete
C it.

         filenm = fgflnm(ifile)
         IF ( qbckfil ) THEN
            extens = filenm(index(filenm,'.')+1:)
            filenm = filenm(:index(filenm,'.')-1)//'_bkg.'//extens         
         ENDIF

         INQUIRE(file=filenm, exist=qexist)
         IF ( qexist ) CALL xdelfil(filenm, status)

         CALL ftinit(iunit, filenm, block, status)
         contxt = 'WRTFAK : Failed to create '//filenm(:lenact(filenm))
         IF( status.NE.0 )GOTO 100
c
c
         simple = .TRUE.
         extend = .TRUE.
 
c
c Define primary header
 
         CALL ftpdef(iunit, bitpix, 0, 0, 0, 1, status)
         contxt = 'WRTFAK : Failed to define primary header'
         IF( status.NE.0 )GOTO 100
 
c
c Write primary array keywords
c
         CALL ftphpr(iunit, simple, bitpix, 0, 0, 0, 1, extend, status)
         contxt = 'WRTFAK : Failed to write primary array keywords'
         IF( status.NE.0 )GOTO 100
c
c Create extension
c
         CALL ftcrhd(iunit, status)
         contxt = 'WRTFAK : Failed to create extension'
         IF( status.NE.0 )GOTO 100
 
c Find out whether the quality and grouping columns need to be
c written and if so set the arrays
 
         qqual = .FALSE.
         qgroup = .FALSE.
         DO 150 i = 1, dgungc(ifile)
            memi(iqualty+i-1) = 0
            memi(igroup+i-1) = 1
            cgroup = dggcrd(ifile, i)
            IF( cgroup.EQ.'-' )THEN
               memi(igroup+i-1) = -1
               qgroup = .TRUE.
            ELSEIF( cgroup.EQ.'*' )THEN
               memi(iqualty+i-1) = 2
               qqual = .TRUE.
            ELSEIF( cgroup.EQ.' ' )THEN
               memi(iqualty+i-1) = 1
               qqual = .TRUE.
            ENDIF
 150     CONTINUE

c Find out whether the BACKSCAL array needs to be written or whether we can
c get away with a scalar keyword

         qback = .FALSE.
         backscal = region(dgnbnb(ifile))
         DO i = dgnbnb(ifile)+1, dgnbne(ifile)
            IF ( ABS(backscal-region(i))/backscal .GT. 1.e-6 ) THEN
               qback = .TRUE.
            ENDIF
         ENDDO

c Do the same for the AREASCAL

         qarea = .FALSE.
         areascal = areasc(dgnbnb(ifile))
         DO i = dgnbnb(ifile)+1, dgnbne(ifile)
            IF ( ABS(areascal-areasc(i))/areascal .GT. 1.e-6 ) THEN
               qarea = .TRUE.
            ENDIF
         ENDDO

c Write required header keywords for extension

         nrows = dgungc(ifile)
         ttype(1) = 'CHANNEL'
         tform(1) = 'I'
         tunit(1) = ' '
         IF(qpoiss)THEN  
            ttype(2) = 'COUNTS'
            tunit(2) = 'COUNTS'
            tform(2) = 'E'
            nfield = 2
         ELSE
            ttype(2) = 'RATE'
            tunit(2) = 'cts/s'
            tform(2) = 'E'
            ttype(3) = 'STAT_ERR'
            tunit(3) = 'cts/s'
            tform(3) = 'E'
            nfield = 3
         ENDIF 
         IF( qqual )THEN
            nfield = nfield + 1
            ttype(nfield) = 'QUALITY'
            tform(nfield) = 'I'
            tunit(nfield) = ' '
         ENDIF
         IF( qgroup )THEN
            nfield = nfield + 1
            ttype(nfield) = 'GROUPING'
            tform(nfield) = 'I'
            tunit(nfield) = ' '
         ENDIF
         IF( qback )THEN
            nfield = nfield + 1
            ttype(nfield) = 'BACKSCAL'
            tform(nfield) = 'E'
            tunit(nfield) = ' '
         ENDIF
         IF( qarea )THEN
            nfield = nfield + 1
            ttype(nfield) = 'AREASCAL'
            tform(nfield) = 'E'
            tunit(nfield) = ' '
         ENDIF
         extnam = 'SPECTRUM'
         pcount = 0

         CALL ftphbn(iunit, nrows, nfield, ttype, tform, tunit, extnam, 
     &               pcount, status)
         contxt = 'WRTFAK : Failed to write extension header keywords'
         IF( status.NE.0 )GOTO 100

         CALL rdrinf(fgauxf(ifile,'r'),ridum,rhdum,ridum,rhdum,rhdum,
     &          rhdum,rrchtyp,ridum,ridum,ridum,rrdum,
     &          ridum,rrdum,ridum,ridum,istat)

c Write additional header keywords for extension
 
         CALL ftpkys(iunit, 'HDUCLASS', 'OGIP', 
     &               'format conforms to OGIP/GSFC conventions', status)
         contxt = 'WRTFAK : Failed to write HDUCLASS keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkys(iunit, 'HDUCLAS1', 'SPECTRUM', 
     &               'extension contains a spectrum', status)
         contxt = 'WRTFAK : Failed to write HDUCLAS1 keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkys(iunit, 'HDUCLAS2', 'TOTAL', 
     &               'spectrum is gross PHA spectrum', status)
         contxt = 'WRTFAK : Failed to write HDUCLAS2 keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkys(iunit, 'HDUCLAS3', 'COUNT', 'spectrum is counts', 
     &               status)
         contxt = 'WRTFAK : Failed to write HDUCLAS3 keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkys(iunit, 'HDUVERS1', '1.1.0', 
     &               'version number of the format', status)
         contxt = 'WRTFAK : Failed to write HDUCLAS3 keyword'
         IF( status.NE.0 )GOTO 100

         name = 'XSPEC '//VERSION//' fakeit' 
         CALL ftpkys(iunit, 'CREATOR', name, 
     &               'program that created this file', status)
         contxt = 'WRTFAK : Failed to write CREATOR keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkys(iunit, 'FILENAME', filenm, 'source file for fit', 
     &               status)
         contxt = 'WRTFAK : Failed to write FILENAME keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkys(iunit, 'DATE    ', cdate, 
     &               'FITS creation date (dd/mm/yyyy)', status)
         contxt = 'WRTFAK : Failed to write DATE keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkys(iunit, 'TELESCOP', dgtlsc(ifile), 
     &               'telescope (mission) name', status)
         contxt = 'WRTFAK : Failed to write TELESCOP keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkys(iunit, 'INSTRUME', dginst(ifile), 
     &               'instrument name', status)
         contxt = 'WRTFAK : Failed to write INSTRUME keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkys(iunit, 'FILTER  ', 'none   ', 
     &               'instrument filter used', status)
         contxt = 'WRTFAK : Failed to write FILTER keyword'
         IF( status.NE.0 )GOTO 100
 
         IF ( qbckfil ) THEN
            CALL ftpkye(iunit, 'EXPOSURE', dgbtme(ifile), 4, 
     &                  'Exposure time', status)
         ELSE
            CALL ftpkye(iunit, 'EXPOSURE', dgtime(ifile), 4, 
     &                  'Exposure time', status)
         ENDIF
         contxt = 'WRTFAK : Failed to write EXPOSURE keyword'
         IF( status.NE.0 )GOTO 100
 
c If we aren't writing a backscal column then put it in the keyword

         IF ( .NOT.qback ) THEN
            CALL ftpkye(iunit, 'BACKSCAL', backscal, 4, 
     &                  'Background scale factor', status)
            contxt = 'WRTFAK : Failed to write BACKSCAL keyword'
            IF( status.NE.0 )GOTO 100
         ENDIF

c If we aren't writing a areascal column then put it in the keyword

         IF ( .NOT.qarea ) THEN
            CALL ftpkye(iunit, 'AREASCAL', areascal, 4, 
     &                  'Area scale factor', status)
            contxt = 'WRTFAK : Failed to write AREASCAL keyword'
            IF( status.NE.0 )GOTO 100
         ENDIF

         CALL ftpkye(iunit, 'CORRSCAL', dgcrnm(ifile), 4, 
     &               'Correction scale factor', status)
         contxt = 'WRTFAK : Failed to write CORRSCAL keyword'
         IF( status.NE.0 )GOTO 100
 
         name = ' '
         CALL ftpkys(iunit, 'OBJECT  ', name, 'object observed', status)
         contxt = 'WRTFAK : Failed to write OBJECT keyword'
         IF( status.NE.0 )GOTO 100
 
         IF ( qbckfil ) THEN
            CALL ftpkys(iunit, 'BACKFILE', 'NONE', 
     &                  'bgrd FITS file', status)
         ELSE
            filenm = fgflnm(ifile)
            extens = filenm(index(filenm,'.')+1:)
            IF ( index(filenm,'.') .GT. 0 .AND. 
     &           dgbtme(ifile) .NE. 0. ) THEN         
               filenm = filenm(:index(filenm,'.')-1)//'_bkg.'//extens
               CALL ftpkys(iunit, 'BACKFILE', filenm, 
     &                     'bgrd FITS file', status)
            ELSE
               CALL ftpkys(iunit, 'BACKFILE', 'NONE', 
     &                     'bgrd FITS file', status)
            ENDIF
         ENDIF
         contxt = 'WRTFAK : Failed to write BACKFILE keyword'
         IF( status.NE.0 )GOTO 100
 
         IF ( qbckfil ) THEN
            CALL ftpkys(iunit, 'CORRFILE', 'NONE', 
     &                  'Correction FITS file', status)
         ELSE
            CALL ftpkys(iunit, 'CORRFILE', fgauxf(ifile,'c'), 
     &                  'Correction FITS file', status)
         ENDIF 
         contxt = 'WRTFAK : Failed to write CORFFILE keyword'
         IF( status.NE.0 )GOTO 100

         CALL ftpkys(iunit, 'RESPFILE', fgauxf(ifile,'r'), 
     &               'Redistribution matrix file (RMF)', status)
         contxt = 'WRTFAK : Failed to write RESPFILE keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkys(iunit, 'ANCRFILE', fgauxf(ifile,'a'), 
     &               'Anc. Resp. file (ARF)', status)
         contxt = 'WRTFAK : Failed to write ANCRFILE keyword'
         IF( status.NE.0 )GOTO 100
         
           
         CALL ftpkys(iunit, 'CHANTYPE', rrchtyp, 
     &               'Chan. assigned by det. electronics', status)
         contxt = 'WRTFAK : Failed to write CHANTYPE keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkyj(iunit, 'DETCHANS', dgungc(ifile), 
     &               'Total no. detector channels available', status)
         contxt = 'WRTFAK : Failed to write DETCHANS keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkyj(iunit, 'TLMIN1', dgchn0(ifile), 
     &               'Minimum value legally allowed in column 1', 
     &               status)
         contxt = 'WRTFAK : Failed to write TLMIN1 keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkyj(iunit, 'TLMAX1', dgungc(ifile)-(1-dgchn0(ifile)), 
     &               'Maximum value legally allowed in column 1', 
     &               status)
         contxt = 'WRTFAK : Failed to write TLMIN1 keyword'
         IF( status.NE.0 )GOTO 100

         IF ( qpoiss ) THEN
            CALL ftpkyj(iunit, 'STAT_ERR', 0, 'No stat error spcf.', 
     &                  status)
            contxt = 'WRTFAK : Failed to write STAT_ERR keyword'
            IF( status.NE.0 )GOTO 100
         ENDIF
 
         CALL ftpkyl(iunit, 'POISSERR', qpoiss, 'Pois. err. assumed ?', 
     &               status)

         contxt = 'WRTFAK : Failed to write POISSERR keyword'
         IF( status.NE.0 )GOTO 100
 
         CALL ftpkyj(iunit, 'SYS_ERR', 0, 'No sys. err spcf.', status)
         contxt = 'WRTFAK : Failed to write SYS_ERR keyword'
         IF( status.NE.0 )GOTO 100
 
         IF( .NOT.qgroup )THEN
            CALL ftpkyj(iunit, 'GROUPING', 0, 'No grpng of data spcf.', 
     &                  status)
            contxt = 'WRTFAK : Failed to write GROUPING keyword'
            IF( status.NE.0 )GOTO 100
         ENDIF
 
         IF( .NOT.qqual )THEN
            CALL ftpkyj(iunit, 'QUALITY', 0, 'No data qual info spcf.', 
     &                  status)
            contxt = 'WRTFAK : Failed to write QUALITY keyword'
            IF( status.NE.0 )GOTO 100
         ENDIF
 
c Define extension data structure
 
         CALL ftbdef(iunit, nfield, tform, pcount, nrows, status)
         contxt = 'WRTFAK : Failed to define extension data structure'
         IF( status.NE.0 )GOTO 100
c
c Write data to extension
c
         chan0 = dgchn0(ifile)
         DO 200 i = 1, nrows
            memi(ichannl+i-1) = (i - 1) + chan0
 200     CONTINUE
         CALL ftpclj(iunit, 1, 1, 1, nrows, memi(ichannl), status)
         contxt = 'WRTFAK : Failed to write CHANNEL data to extension'
         IF( status.NE.0 )GOTO 100
 
c Set up the COUNTS/RATE&ERROR arrays with zero values included where a channel
c was ignored or grouped up.

         IF ( qbckfil ) THEN
            grasp = dgbtme(ifile)
         ELSE
            grasp = dgtime(ifile)
         ENDIF

         ibin = dgnbnb(ifile) - 1
         ichan = 0
         DO 250 i = 1, nrows
            memr(icounts+i-1) = 0.
            memr(ierror+i-1) = 0.
            IF( (dggcrd(ifile,i).EQ.'+' .OR. dggcrd(ifile,i).EQ.'*') )
     &          THEN
               ichan = ichan + 1
               IF( dgqnot(ifile,ichan) )THEN
                  ibin = ibin + 1
                  IF ( qpoiss ) THEN
                     memr(icounts+i-1) = NINT(grasp*areasc(ibin)
     &                                             *phaobs(ibin))
                  ELSE
                     memr(icounts+i-1) = areasc(ibin)*phaobs(ibin)
                     memr(ierror+i-1) = areasc(ibin)*sqrt(varobs(ibin))
                  ENDIF
               ENDIF
            ENDIF
 250     CONTINUE
 
         CALL ftpcle(iunit, 2, 1, 1, nrows, memr(icounts), status)

         IF(qpoiss)THEN
            contxt='WRTFAK : Failed to write COUNTS data to extension'
         ELSE
            contxt='WRTFAK : Failed to write RATE data to extension'
         ENDIF  
         IF( status.NE.0 )GOTO 100
         icol = 2

         IF ( .NOT. qpoiss ) THEN
            icol = icol + 1
            CALL ftpcle(iunit, icol, 1, 1, nrows, memr(ierror), status)
            contxt='WRTFAK : Failed to write STAT_ERR data to extension'
            IF( status.NE.0 )GOTO 100
         ENDIF
 
         IF( qqual )THEN
            icol = icol + 1
            CALL ftpclj(iunit, icol, 1, 1, nrows, memi(iqualty), status)
            contxt = 
     &              'WRTFAK : Failed to write QUALITY data to extension'
            IF( status.NE.0 )GOTO 100
         ENDIF
 
         IF ( qgroup )THEN
            icol = icol + 1
            CALL ftpclj(iunit, icol, 1, 1, nrows, memi(igroup), status)
            contxt = 
     &             'WRTFAK : Failed to write GROUPING data to extension'
            IF( status.NE.0 )GOTO 100
         ENDIF

c If we need to write BACKSCAL as a vector then do so. As with the actual
c data we need to put values in for channels that have been ignored or
c grouped up. For these channels we set the BACKSCAL to that for the last
c valid bin. This is not perfect but will have to do.

         IF ( qback ) THEN
            icol = icol + 1
            ibin = dgnbnb(ifile) - 1
            ichan = 0
            DO i = 1, nrows
               IF ( ibin .LT. dgnbnb(ifile) ) THEN
                  memr(icounts+i-1) = region(dgnbnb(ifile))
               ELSE
                  memr(icounts+i-1) = region(ibin)
               ENDIF
               IF( (dggcrd(ifile,i).EQ.'+' .OR. 
     &              dggcrd(ifile,i).EQ.'*') ) THEN
                  ichan = ichan + 1
                  IF( dgqnot(ifile,ichan) )THEN
                     ibin = ibin + 1
                     memr(icounts+i-1) = region(ibin)
                  ENDIF
               ENDIF
            ENDDO
 
            CALL ftpcle(iunit, icol, 1, 1, nrows, memr(icounts), 
     &                  status)
            contxt = 
     &             'WRTFAK : Failed to write BACKSCAL data to extension'
            IF( status.NE.0 )GOTO 100
         ENDIF

c If we need to write AREASCAL as a vector then do so. As with the actual
c data we need to put values in for channels that have been ignored or
c grouped up. For these channels we set the AREASCAL to that for the last
c valid bin. This is not perfect but will have to do.

         IF ( qarea ) THEN
            icol = icol + 1
            ibin = dgnbnb(ifile) - 1
            ichan = 0
            DO i = 1, nrows
               IF ( ibin .LT. dgnbnb(ifile) ) THEN
                  memr(icounts+i-1) = areasc(dgnbnb(ifile))
               ELSE
                  memr(icounts+i-1) = areasc(ibin)
               ENDIF
               IF( (dggcrd(ifile,i).EQ.'+' .OR. 
     &              dggcrd(ifile,i).EQ.'*') ) THEN
                  ichan = ichan + 1
                  IF( dgqnot(ifile,ichan) )THEN
                     ibin = ibin + 1
                     memr(icounts+i-1) = areasc(ibin)
                  ENDIF
               ENDIF
            ENDDO
 
            CALL ftpcle(iunit, icol, 1, 1, nrows, memr(icounts), 
     &                  status)
            contxt = 
     &             'WRTFAK : Failed to write AREASCAL data to extension'
            IF( status.NE.0 )GOTO 100
         ENDIF

c Write out error message if necessary
 
 100     CONTINUE
         IF( status.NE.0 )THEN
            CALL xwrite(contxt, 10)
            WRITE(contxt, '(a, i4, a, i4)') 'FITSIO error = ', status,
     &           'for file number ', ifile
            CALL xwrite(contxt, 10)
         ENDIF
c
c Close file
c
         CALL ftclos(iunit, status)
 
 300     CONTINUE

      ENDDO

      CALL frelun(iunit)
 
      CALL udmfre(icounts, 6, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     &     ' WARNING: Memory allocation failed in WRTFAK.', 2)      

      CALL udmfre(ierror, 6, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     &     ' WARNING: Memory allocation failed in WRTFAK.', 2)      

      CALL udmfre(ichannl, 4, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     &     ' WARNING: Memory allocation failed in WRTFAK.', 2)      

      CALL udmfre(iqualty, 4, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     &     ' WARNING: Memory allocation failed in WRTFAK.', 2)      

      CALL udmfre(igroup, 4, istat )
      IF( istat .NE. 0 ) CALL xwrite(
     &     ' WARNING: Memory allocation failed in WRTFAK.', 2)      

      RETURN
      END