C*****************************************************************************
C TASK NAME: econvpha
C
C FILE NAME: econvpha.f
C
C DESCRIPTION: FTOOL to convert EGRET .spec (or .specin) files into FITS
C     files in XSPEC format.  The .spec files are generated by the EGRET
C     data analysis codes SKYUTIL and LIKE.
C
C AUTHOR/DATE: Peter J.T. Leonard, HSTX, 01/97 and 03/97
C              Jeffrey M.S. Silvis, HSTX, 03/97
C
C NOTES:
C
C USAGE:
C     Host: econvpha
C     IRAF: econvpha
C
C ROUTINES IMPLEMENTED IN THIS FILE:
C     subroutine econva - top level subroutine, called by IRAF or host
C                         C wrapper in heconvpha.c
C     subroutine readat - reads in EGRET file from SKYUTIL and LIKE
C     subroutine wftspf - writes out FITS file in XSPEC format
C
C MODIFICATION HISTORY:
C     Ning Gan  1-July-1998: Changed the comments for DATEx keywords.
C*****************************************************************************

C*****************************************************************************
C SUBROUTINE: econva
C
C DESCRIPTION: Program to convert EGRET .spec (or .specin) files into FITS
C     files in XSPEC format.  The .spec files are generated by the EGRET data
C     analysis codes SKYUTIL and LIKE.
C
C AUTHOR/DATE: Keith A. Arnaud, 02/93
C              John R. Mattox, 07/94
C
C NOTES: Based on Keith A. Arnaud's subroutine P2FSUB, but uses it in reverse.
C     Consequently, phafil is poorly named; the input file is actually a .spec
C     file, not a .pha file.  It is the output file that is the .pha file.
C
C ARGUMENTS:
C     phafil - SKYUTIL and LIKE .spec filename
C     outfil - XSPEC .pha filename
C
C PRIMARY LOCAL VARIABLES:
C
C CALLED ROUTINES:
C     subroutine readat - reads in EGRET file from SKYUTIL and LIKE 
C     subroutine wftspf - writes out FITS file in XSPEC format
C
C MODIFICATION HISTORY:
C     01/97, cosmetic changes, Peter J.T. Leonard, HSTX
C     03/97, viewing period information, PJTL and JMSS, HSTX
C*****************************************************************************

      subroutine econva ()

      implicit none
 
      INTEGER*4 ierr, status, MAXCHN, MAXHST
      INTEGER*4 dtype, nchan, i, nhist, lenact
      PARAMETER(MAXCHN=4096, MAXHST=50)
      INTEGER*2 group(MAXCHN), qualty(MAXCHN), chan(MAXCHN)
      INTEGER*2 ipha(MAXCHN)
      EXTERNAL lenact

      REAL*4 pha(MAXCHN), errors(MAXCHN), sysfrc(MAXCHN)
      REAL*4 equinx, ranom, decnom, roll, texpos
      REAL*4 ascale, bscale, cscale
 
      CHARACTER phafil*255, outfil*255
      CHARACTER*255 assfls(4)
      CHARACTER*80 sdsstr, edsstr, shsstr, ehsstr, objnam
      CHARACTER*80 tlscpe, instrm, filter
      CHARACTER*80 hist(MAXHST)
      CHARACTER*70 charlines(100)
 
      LOGICAL qerror, qsys, qqual, qgroup

C Setup common block for taskname--necessary for fcerr.
      CHARACTER*80 msg
      CHARACTER*40 taskname
      common /task/taskname
      taskname = 'econvpha'

      hist(1) = 'XSPEC PHA file produced from EGRET LIKElihood analysis'
      nhist = 1

C Note that phafil is poorly named.  The input file is actually a .spec
C file, not a .pha file.  It is the output file that is the .pha file.
 
      status = 0
      call uclgst ('phafil', phafil, status)
      IF (status.NE.0) THEN
         msg = 'Problem getting name of input file!'
         call fcerr (msg)
         msg = 'Aborting program.'
         call fcerr (msg)
         call exit (1)
      END IF
 
      status = 0
      call uclgst ('outfil', outfil, status)
      IF (status.NE.0) THEN
         msg = 'Problem getting name of output file!'
         call fcerr (msg)
         msg = 'Aborting program.'
         call fcerr (msg)
         call exit (1)
      END IF

C Read the EGRET .spec file produced by SKYUTIL and LIKE.
      open (77, file=phafil, status='old')
         call readat (77, nchan, pha, errors, ranom, decnom, charlines)
      close (77)

      status = 0
      call uclgst ('objnam', objnam, status)
      IF (status.NE.0) THEN
         msg = 'Problem getting name of observed object!'
         call fcerr (msg)
         msg = 'Aborting program.'
         call fcerr (msg)
         call exit (1)
      END IF

C Read the grouping card and set the quality and grouping flags.
      DO i = 1, nchan
         qualty(i) = 0
      ENDDO
      qqual = .FALSE.
      qgroup = .FALSE.

C Set the array giving the channel numbers.  These are assumed to be
C monotonic starting at one.
      DO i = 1, nchan
         chan(i) = i
      ENDDO

      tlscpe = 'EGRET'
      instrm = 'EGRET'
      equinx = 2000.0
      qerror = .true.
C Systematic errors not included.
      qsys = .false.
      filter = 'none'
C Roll angle irrelevant.
      roll = 0.0

C The EGRET data for spectral analysis are likelihood counts estimates
C and uncertainties of these estimates.  However, XSPEC does not allow
C specification of a counts uncertainty, but it does allow for specification
C of a counts rate uncertainty.  Thus, we have to do some pretending...
 
C Say it is counts/sec.
      dtype = 2
 
C Also, lie about the exposure time to get the flux right in the end.
      texpos = 1.0
	
C EGRET does not have associated files (back, RMF, corr, ARF).
      DO i = 1, 4
         assfls(i) = 'none'
      ENDDO

C EGRET spectra are often obtained from multiple observations, while XSPEC
C accomodates only a single start and stop.
      sdsstr = ' '
      shsstr = ' '
      edsstr = ' '
      ehsstr = ' '

C Area, background, and correction scale factors.
      ascale = 1.0
      bscale = 1.0
      cscale = 0.0

C Write the FITS file.
      status = 0
      call wftspf (outfil, phafil, hist, nhist, assfls, equinx,
     &             sdsstr, shsstr, edsstr, ehsstr, objnam, tlscpe,
     &             instrm, filter, ranom, decnom, roll, texpos,
     &             ascale, bscale, cscale, qgroup, group, qqual,
     &             qualty, nchan, dtype, chan, ipha, pha, qerror,
     &             errors, qsys, sysfrc, charlines, ierr, status)
      IF (status.NE.0) THEN
         msg = 'Problem with writing output file!'
         call fcerr (msg)
         msg = 'Aborting program.'
         call fcerr (msg)
         call exit (1)
      END IF

      RETURN
 
      END
 
C*****************************************************************************
C SUBROUTINE: readat
C
C DESCRIPTION: Reads in EGRET .spec file produced by SKYUTIL and LIKE.
C
C AUTHOR/DATE: Patrick L. Nolan, 04/92
C              John R. Mattox, 07/94
C
C NOTES:
C
C ARGUMENTS:
C     lusel      -  fortran unit number
C     nregns     -  number of acceptance cones
C     nmeas1     -  counts in cones
C     sigma1     -  uncertainty in counts in cones
C     ra, dec    -  r.a. and dec. of pointing axis
C     charlines  -  viewing period information
C
C PRIMARY LOCAL VARIABLES:
C
C CALLED ROUTINES: None.
C
C MODIFICATION HISTORY:
C     01/97, cosmetic changes, Peter J.T. Leonard, HSTX
C     03/97, viewing period information, PJTL and JMSS, HSTX
C******************************************************************************
 
      subroutine readat (lusel, nregns, nmeas1, sigma1, ra, dec,
     &                   charlines)

      implicit none

      integer*4 first, last, spacing, jmax
C First bin in incident spectrum.
      parameter (first=-9)
C Last bin in incident spectrum.
      parameter (last=260)
C Incident bins between 20 MeV and 20 GeV.
      parameter (spacing=200)
C Maximum number of observed bins.
      parameter (jmax=15)
C Number of acceptance cones.
      integer*4 nregns
      integer*4 lusel, n, irec
 
      real*4 nmeas1(0:jmax), sigma1(0:jmax)
      real*4 ra, dec
      real*4 emin, emax, counts, err

      character*70 charlines(100)
      character*70 charline
      character*8 keywd
      equivalence (keywd, charline)

C Read header lines from file and compare with header stored in matrix file.
 
C Get a line from data file.
      read (lusel, '(a)') charline
      irec = 1
      charlines(irec) = charline
      do while (keywd.ne.'.......e'.and.keywd.ne.'End time')
         if (keywd.eq.'Source p') then
            read (charline(28:46), *) ra, dec
         end if
	 irec = irec + 1
         read (lusel, '(a)') charline
         charlines(irec) = charline
      end do
      charlines(irec-1) = 'the end'

C First value is at beginning of array for XSPEC.
      n = -1
C No class C.  May fix this some day.
      nmeas1(0) = 0.0
      sigma1(0) = 0.0
C Read data lines from file.
      do while (.true.)
         read (lusel, *, end=999) emin, emax, counts, err
         if (emin.le.0.) go to 999
         if (n.ge.jmax) go to 999
         n = n + 1
         sigma1(n) = err
         nmeas1(n) = counts
      end do
999   continue

      nregns = n + 1
 
      return
 
      end
 
C*****************************************************************************
C SUBROUTINE: wftspf
C
C DESCRIPTION: Writes out FITS file in XSPEC format.
C
C AUTHOR/DATE: Keith A. Arnaud, 199?
C
C NOTES:
C
C ARGUMENTS:
C     outfil     -  FITS filename
C     phafil     -  input PHA filename
C     sddstr     -  start day of observation
C     shhstr     -  start time of observation
C     eddstr     -  end day of observation
C     ehhstr     -  end time of observation
C     objnam     -  object observed
C     tlscpe     -  telescope/satellite name
C     instrm     -  instrument name
C     filter     -  instrument filter
C     hist(*)    -  history records
C     assfls(4)  -  associated files (back, RMF, corr, ARF)
C     nhist      -  number of history records
C     nchan      -  number of channels
C     dtype      -  date type  1=counts, 2=counts/sec
C     group(*)   -  channel group
C     qualty(*)  -  channel quality
C     chan(*)    -  channel numbers
C     ipha(*)    -  data in counts
C     pha(*)     -  data in counts/sec
C     errors(*)  -  errors in counts/sec
C     sysfrc(*)  -  fractional systematics
C     equinx     -  observation epoch
C     ranom      -  observation RA (decimal degrees)
C     decnom     -  observation DEC (decimal degrees)
C     roll       -  observation roll angle (decimal degrees)
C     texpos     -  exposure time
C     ascale     -  area scaling factor
C     bscale     -  background scaling factor
C     cscale     -  correction scaling factor
C     qgroup     -  data is grouped
C     qqual      -  data has quality flags
C     qerror     -  statistical errors included
C     charlines  -  viewing period information
C     qsys       -  systematic errors included
C     ierr       -  error
C                      1 = failed to open FITS file
C                      2 = failed to write primary header
C                      3 = failed to set up SPECTRUM extension
C                      4 = failed to write XSPEC keywords
C                      5 = failed to write extra extension keywords
C                      6 = failed to write SPECTRUM data
C                      7 = failed to close FITS file
C     status     -  FITSIO status
C
C PRIMARY LOCAL VARIABLES:
C
C CALLED ROUTINES:
C
C MODIFICATION HISTORY:
C     01/97, cosmetic changes, Peter J.T. Leonard, HSTX
C     03/97, viewing period information, PJTL and JMSS, HSTX
C*****************************************************************************
 
      SUBROUTINE wftspf (outfil, phafil, hist, nhist, assfls, equinx,
     &                   sdsstr, shsstr, edsstr, ehsstr, objnam, tlscpe,
     &                   instrm, filter, ranom, decnom, roll, texpos, 
     &                   ascale, bscale, cscale, qgroup, group, qqual, 
     &                   qualty, nchan, dtype, chan, ipha, pha, qerror,
     &                   errors, qsys, sysfrc, charlines, ierr, status)

      implicit none

      INTEGER*4 nhist, nchan, dtype, ierr, status
      INTEGER*4 ounit, i, tfields, colnum, lenact
      INTEGER*2 group(*), qualty(*), chan(*), ipha(*)
 
      REAL*4 pha(*), errors(*), sysfrc(*)
      REAL*4 equinx, ranom, decnom, roll, texpos
      REAL*4 ascale, bscale, cscale

      CHARACTER*(*) outfil, phafil
      CHARACTER*(*) sdsstr, shsstr, edsstr, ehsstr
      CHARACTER*(*) objnam, tlscpe, instrm, filter
      CHARACTER*(*) hist(*), assfls(4)
      CHARACTER*(*) charlines(100)
      CHARACTER*20 ttype(6), tunit(6)
      CHARACTER*5  tform(6)
      CHARACTER*7 keywd
 
      LOGICAL qgroup, qqual, qerror, qsys

C Open the output FITS file.
      CALL getlun (ounit)
      CALL ftinit (ounit, outfil, 1, status) 
      IF (status .NE. 0) THEN
         ierr = 1
         RETURN
      ENDIF

C Write the primary header.
      CALL ftpdef (ounit, 8, 0, 0, 0, 1, status)

C Write the basic primary array keywords.
      CALL ftphpr (ounit, .TRUE., 8, 0, 0, 0, 1, .TRUE., status)

C Write out the additional keywords about the creation of the FITS file.
      CALL ftpkys (ounit, 'CONTENT', 'SPECTRUM',
     &        'spectrum file contains time intervals and event', status)
 
      CALL ftpkys (ounit, 'FILENAME', phafil(:lenact(phafil)),
     &             'File that FITS was produced from', status)
 
      CALL ftpkys (ounit, 'ORIGIN', 'NASA/GSFC',
     &             'origin of FITS file', status)

C Write the TELESCOP and INSTRUME keywords.
      CALL ftpkys (ounit,'TELESCOP', tlscpe,
     &             'Telescope (mission) name', status)
 
      CALL ftpkys (ounit,'INSTRUME', instrm,
     &             'Instrument name', status)

C Write the object name.
      CALL ftpkys (ounit, 'OBJECT', objnam,
     &             'Name of observed object', status)

C Write the pointing data.
      CALL ftpkye (ounit,'RA-NOM', ranom, 4,
     &             'Right Ascension of target (deci. deg)', status)
 
      CALL ftpkye (ounit,'DEC-NOM', decnom, 4,
     &             'Declination of target (deci. deg)', status)
 
      CALL ftpkye (ounit,'DROLLANG', roll, 4,
     &             'Mean roll angle (deci. deg)', status)

C Write the start and end times.
      CALL ftpkys (ounit, 'DATE-OBS', sdsstr,
     &             'Date observations were made (yyyy-mm-dd)', status)
 
      CALL ftpkys (ounit, 'TIME-OBS', shsstr,
     &             'Time observations were made (hh:mm:ss)', status)
 
      CALL ftpkys (ounit, 'DATE-END', edsstr,
     &             'Date observations were made (yyyy-mm-dd)', status)
 
      CALL ftpkys (ounit, 'TIME-END', ehsstr,
     &             'Time observations were made (hh:mm:ss)', status)
 
C Write out viewing period information.
      i = 5
      keywd = 'not end'
      DO WHILE (keywd .ne. 'the end')
         CALL ftphis (ounit, charlines(i), status)
         write (keywd,100) charlines(i+1)
100      format (a7)
         i = i + 1
      END DO

C Write any history records to the FITS file as history.
      DO i = 1, nhist
         CALL ftphis (ounit, hist(i), status)
      ENDDO

      IF (status .NE. 0) THEN
         ierr = 2
         RETURN
      ENDIF

C Write the data extension with the SPECTRUM data.

C Create the data extension.
      CALL ftcrhd (ounit, status)

C Set up the header keywords for a binary extension.
      tfields = 2
      ttype(1) = 'CHANNEL'
      tform(1) = 'I'
      tunit(1) = ' '
 
      IF (dtype .EQ. 1) THEN
         ttype(2) = 'COUNTS'
         tform(2) = 'I'
         tunit(2) = 'counts'
      ELSEIF (dtype .EQ. 2) THEN
         ttype(2) = 'RATE'
         tform(2) = 'E'
         tunit(2) = 'counts/sec'
      ENDIF
 
      IF (qerror) THEN
         tfields = tfields + 1
         ttype(tfields) = 'STAT_ERR'
         tform(tfields) = 'E'
         tunit(tfields) = 'counts/sec'
      ENDIF
 
      IF (qsys) THEN
         tfields = tfields + 1
         ttype(tfields) = 'SYS_ERR'
         tform(tfields) = 'E'
         tunit(tfields) = ' '
      ENDIF
 
      IF (qqual) THEN
         tfields = tfields + 1
         ttype(tfields) = 'QUALITY'
         tform(tfields) = 'I'
         tunit(tfields) = ' '
      ENDIF
 
      IF (qgroup) THEN
         tfields = tfields + 1
         ttype(tfields) = 'GROUPING'
         tform(tfields) = 'I'
         tunit(tfields) = ' '
      ENDIF

C Write the main header keywords.
      CALL ftphbn (ounit, nchan, tfields, ttype, tform, tunit,
     &             'SPECTRUM', 0, status)

C Write the additional keywords describing the data stored.
      IF (dtype .EQ. 1) THEN
         CALL ftpkyl (ounit, 'POISSERR', .TRUE.,
     &                'Poissonian errors to be assumed', status)
      ENDIF
 
      IF (.NOT. qerror) THEN
         CALL ftpkyj (ounit, 'STAT_ERR', 0,
     &                'no statistical error specified', status)
      ENDIF
 
      IF (.NOT. qsys) THEN
         CALL ftpkyj (ounit, 'SYS_ERR', 0,
     &                'no systematic error specified', status)
      ENDIF
 
      IF (.NOT. qgroup) THEN
         CALL ftpkyj (ounit, 'GROUPING', 0,
     &                'no grouping of the data has been defined', 
     &                status)
      ENDIF
 
      IF (.NOT. qqual) THEN
         CALL ftpkyj (ounit, 'QUALITY', 0,
     &                'no data quality information specified', status)
      ENDIF

      IF (status .NE. 0) THEN
         ierr = 3
         RETURN
      ENDIF

C Write the XSPEC mandatory keywords.

C Write the TELESCOP and INSTRUME keywords.
      CALL ftpkys (ounit, 'TELESCOP', tlscpe,
     &             'Telescope (mission) name', status)
 
      CALL ftpkys (ounit,'INSTRUME', instrm,
     &             'Instrument name', status)

C Write the filter keyword.
      CALL ftpkys (ounit, 'FILTER', filter,
     &             'Instrument filter in use', status)

C Write the exposure keyword.
      CALL ftpkye (ounit, 'EXPOSURE', texpos, 4,
     &             'Exposure time', status)

C Write the area scaling factor.
      CALL ftpkye (ounit, 'AREASCAL', ascale, 4,
     &             'nominal effective area', status)

C Write the background scaling factor.
      CALL ftpkye (ounit, 'BACKSCAL', bscale, 4,
     &             'background scale factor', status)

C Write the correlation scaling factor.
      CALL ftpkye (ounit, 'CORRSCAL', cscale, 4,
     &             'Correction scale factor', status)

C Write the background, correction, RMF and ARF filenames.
      CALL ftpkys (ounit, 'BACKFILE', assfls(1),
     &             'background FITS file for this object', status)
 
      CALL ftpkys (ounit, 'CORRFILE', assfls(3),
     &             'correction FITS file for this object', status)
 
      CALL ftpkys (ounit, 'RESPFILE', assfls(2),
     &             'redistribution matrix (RMF)', status)
 
      CALL ftpkys (ounit, 'ANCRFILE', assfls(4),
     &             'ancillary response matrix (ARF)', status)

C Write any XSPEC filter keywords.

C Write the channel type.
      CALL ftpkys (ounit, 'CHANTYPE', 'PHA',
     &             'Channels assigned by detector electronics', status)

C Write the number of detector channels.
      CALL ftpkyj (ounit, 'DETCHANS', nchan,
     &  	   'Total no. detector channels available', status)

C Write the FITS format used.
      CALL ftpkys (ounit, 'PHAVERSN', '1992a',
     &             'OGIP classification of FITS format style', status)

      IF (status .NE. 0) THEN
         ierr = 4
         RETURN
      ENDIF

C Write out optional keywords - these are mostly repeats of those 
C in the primary header.
      CALL ftpkys (ounit, 'OBJECT', objnam,
     &             'Name of observed object', status)

      CALL ftpkys (ounit, 'FILENAME', phafil(:lenact(phafil)),
     &             'File that FITS was produced from', status)

      CALL ftpkys (ounit, 'ORIGIN', 'NASA/GSFC',
     &             'origin of FITS file', status)

      CALL ftpdat (ounit, status)

      CALL ftpkys (ounit,'TELESCOP', tlscpe,
     &             'Telescope (mission) name', status)

      CALL ftpkys (ounit,'INSTRUME', instrm,
     &             'Instrument name', status)

      CALL ftpkye (ounit,'RA-NOM', ranom, 4,
     &             'Right Ascension of target (deci. deg)', status)

      CALL ftpkye (ounit,'DEC-NOM', decnom, 4,
     &             'Declination of target (deci. deg)', status)

      CALL ftpkye (ounit,'DROLLANG', roll, 4,
     &             'Mean roll angle (deci. deg)', status)

      CALL ftpkyf (ounit, 'EQUINOX', equinx, 1,
     &             'Equinox of celestial coord system', status)

      CALL ftpkys (ounit, 'DATE-OBS', sdsstr,
     &             'Date observations were made (yyyy-mm-dd)', status)
 
      CALL ftpkys (ounit, 'TIME-OBS', shsstr,
     &             'Time observations were made (hh:mm:ss)', status)

      CALL ftpkys (ounit, 'DATE-END', edsstr,
     &             'Date observations were made (yyyy-mm-dd)', status)
 
      CALL ftpkys (ounit, 'TIME-END', ehsstr,
     &             'Time observations were made (hh:mm:ss)', status)

      IF (status .NE. 0) THEN
         ierr = 5
         RETURN
      ENDIF

C Write out the SPECTRUM data.

C Define the binary table.
      CALL ftbdef (ounit, tfields, tform, 0, nchan, status)

C Write the elements into the table.
      CALL ftpcli (ounit, 1, 1, 1, nchan, chan, status)

      IF (dtype .EQ. 1) THEN
         CALL ftpcli (ounit, 2, 1, 1, nchan, ipha, status)
      ELSEIF (dtype .EQ. 2) THEN
         CALL ftpcle (ounit, 2, 1, 1, nchan, pha, status)
      ENDIF

      colnum = 2
      IF (qerror) THEN
         colnum = colnum + 1
         CALL ftpcle (ounit, colnum, 1, 1, nchan, errors, status)
      ENDIF
      IF (qsys) THEN
         colnum = colnum + 1
         CALL ftpcle (ounit, colnum, 1, 1, nchan, sysfrc, status)
      ENDIF
      IF (qqual) THEN
         colnum = colnum + 1
         CALL ftpcli (ounit, colnum, 1, 1, nchan, qualty, status)
      ENDIF
      IF (qgroup) THEN
         colnum = colnum + 1
         CALL ftpcli (ounit, colnum, 1, 1, nchan, group, status)
      ENDIF

      IF (status .NE. 0) THEN
         ierr = 6
         RETURN
      ENDIF
      
C Close the FITS file.
      CALL ftclos (ounit, status)

      IF (status .NE. 0) THEN
         ierr = 7
         RETURN
      ENDIF

      RETURN
 
      END
