

C **********************************************************************
C SUBROUTINE:
C      readlc
C
C DESCRIPTION:
C      Read a  Standard OGIP Timing FITS files.
C
C AUTHOR:
C      James Lochner 4/29/94
C
C MODIFICATION HISTORY:
C      10/3/94 - allow user to specify names of time, rate and error
C                 columns.  Defaults remain TIME, RATE, ERROR
C      10/7/94 - implemented abort option when RATE (and ERROR) columns
C                 not found
C      12/19/1997 Jeff Guerber - Check for nrows > maxsize.
C
C NOTES:
C
C USAGE:
C      subroutine readlc(iunit, infile, inopen, block, nrows,
C                  mintime, maxtime, timename, ratename, errname,
C                  timecol, time, rate, error, timedel, runit,
C                  timeunit, tzero, maxsize)
C
C ARGUMENTS:
C      iunit    - logical unit for input file
C      infile   - name of input file
C      inopen   - logical indicating open status of input file
C      block    -
C      nrows    - number of rows in input file (= number of data points)
C      mintime  - absolute minimum time in the data (in same units as TIMEZERO)
C      maxtime  - absolute maximum time in the data (in same units as TIMEZERO)
C      timename - alternate name for time column in infile
C      ratename - alternate name for intensity column in infile
C      errname  - alternate name for error column in infile
C      timecol  - flag indicating presence of a time column in input file
C      time     - array of input time values (in units of timeunit)
C      rate     - array of input rate (intensity) values
C      error    - array of input error values
C      timedel  - bin size of the data (on output, in units of timeunit)
C      runit    - units associated with the RATE column
C      timeunit - units of the time array
C      tzero    - MJDREF + TIMEZERO value
C      maxsize  - maximum size for output arrays
C      abort    - status flag for non-fitsio errors
C
C PRIMARY LOCAL VARIABLES:
C      ttype    - array of column names
C      tunit    - array of column units
C      spanunit - value of the TIMEUNIT keyword - input units for timedel
C                  and timezero
C      mjdref   - mjdref value from input extenson
C      timecol  - logical indicating presence of time column
C      ratecol  - logical indicating presence of rate column
C      errcol   - logical indicating presence of error column
C      timepos  - column number of the time column
C      ratepos  - column number of the rate column
C      errpos   - column number of the error column
C      ftstatus - fitsio status flag
C      context  - error messages generated by this subroutine
C
C CALLED ROUTINES:
C      subroutine fcpars - parse the input file name and extension number
C      function   fcstln - position of last non-space in a string
C      subroutine ftopen - opens the input FITS file
C      subroutine ftmrhd - relative move to FITS extension
C      subroutine ftgkns -
C      subroutine ftgkyd - get real*8 valued keyword
C      subroutine ftgkyj - getinteger valued keyword
C      subroutine ftgcvd - read real*8 values from table column
C      subroutine ftgcve - read real*4 values from table column
C
C *********************************************************************************
	subroutine readlc(iunit, infile, inopen, block, nrows,
     &     mintime, maxtime, timename, ratename, errname, timecol,
     &     time, rate, error, timedel, runit, timeunit, tzero, maxsize,
     &     abort)


C data declarations
	integer maxsize
	integer iunit, block, nrows

	character *(*) infile, runit, timeunit
        character *(*) timename, ratename, errname
 	double precision time(maxsize), timedel, mintime, maxtime
	real rate(maxsize), error(maxsize)
	logical abort, inopen, timecol

	logical anynull, ratecol, errcol
	character*8 keyroot, ttype(10), tunit(10), spanunit
	character*160 filename
	character*80 errstr
 	character*70 context, comment
	double precision nlvald, timezf, tzero, mjdref, mjdreff
        real nlvale
        integer timepos, ratepos, errpos, mjdrefi
	integer startno, maxkeys, ncol, mcol, timezi
	integer j, ftstatus, frow, felem, extnum, htype
	integer fcstln

        startno = 1
        maxkeys = 10
        anynull = .false.
        nlvald = -100.
        nlvale = -100.
        frow = 1
        felem = 1
        timecol = .false.
        ratecol = .false.
	ftstatus=0
	errcol = .false.

c get the filename and extension number
c  ? What does this do that the par file didn't ?
	call fcpars(infile,filename,extnum,ftstatus)

c open the input FITS file
	call ftopen(iunit,filename,0,block,ftstatus)
	if (ftstatus .ne. 0) then
	    context = 'unable to open infile'
	    call fcerr(context)
	    goto 999
	endif
 	inopen = .true.

c move to the extension
	call ftmrhd(iunit,1,htype,ftstatus)
	if (ftstatus .ne. 0) then
	  errstr = 'error moving to extension number '
	  call fcerr(context)
	  goto 999
	endif

c      Assumed a BINNED light curve rather than an EVENT list

c      Get info on number of columns and their content
        keyroot = 'TTYPE'
        call ftgkns(iunit,keyroot,startno,maxkeys,ttype,ncol,ftstatus)
        if (ftstatus .ne. 0) then
           context = 'error reading TTYPE keywords'
           call fcerr(context)
           goto 999
        endif

        keyroot = 'TUNIT'
        call ftgkns(iunit,keyroot,startno,maxkeys,tunit,mcol,ftstatus)
        if (ftstatus .ne. 0) then
           context = 'error reading TUNIT keywords'
           call fcerr(context)
           goto 999
        endif

c      get the number of rows from keyword
        call ftgkyj(iunit,'naxis2',nrows,comment,ftstatus)
        if (ftstatus .ne. 0) then
           errstr = 'error reading nrows keyword from FITS file'
           call fcerr(errstr)
           goto 999
        endif
	if (nrows .gt. maxsize) then
	    write (errstr, '(3a,i6)') 'Warning: ',
     &          infile(1:fcstln(infile)),
     &          ' has too many rows, truncating at ', maxsize
	    call fcerr(errstr)
	    nrows = maxsize
	endif

C      get TIMEUNIT value (units for TIMEZER and TIMEDEL)
       call ftgkys(iunit,'TIMEUNIT',spanunit,comment,ftstatus)
       if (ftstatus .ne. 0) then
          errstr = 'error getting TIMEUNIT value'
          call fcerr(errstr)
          goto 999
       endif

c      get the bin size of the data (in TIMEUNIT units)
       call ftgkyd(iunit,'TIMEDEL',timedel,comment,ftstatus)
       if (ftstatus .ne. 0) then
          errstr = 'error getting TIMEDEL value'
          call fcerr(errstr)
          goto 999
       endif

C      get MJDREF value (if it doesn't exist simply keep it at 0.0)
c      Read the integer and fractional pieces of MJDREF.
c      If they do not exist (ftstatus .eq. 202), then look for the
c      MJDREF keyword.
       mjdref = 0.0
       call ftgkyj(iunit,'MJDREFI',mjdrefi,comment,ftstatus)
       if (ftstatus .eq. 202) then
          ftstatus = 0
          call ftgkyd(iunit,'MJDREF',mjdref,comment,ftstatus)
          if (ftstatus .ne. 0 .and. ftstatus .ne. 202) then
             errstr = 'error getting MJDREF value'
             call fcerr(errstr)
             goto 999
          endif
       else
          if (ftstatus .ne. 0) then
             errstr = 'error getting MJDREFI value'
             call fcerr(errstr)
             goto 999
          endif
          call ftgkyd(iunit,'MJDREFF',mjdreff,comment,ftstatus)
          if (ftstatus .ne. 0) then
             errstr = 'error getting MJDREFF value'
             call fcerr(errstr)
             goto 999
          endif
          mjdref = mjdrefi + mjdreff
       endif

c      Read the integer and fractional pieces of TIMEZER.
c      If they do not exist (ftstatus .eq. 202), then look for the
c      TIMEZERO keyword.
         call ftgkyj(iunit,'TIMEZERI',timezi,comment,ftstatus)
         if (ftstatus .eq. 202) then
            ftstatus = 0
            call ftgkyd(iunit,'TIMEZERO',tzero,comment,ftstatus)
            if (ftstatus .ne. 0) then
               errstr = 'error getting TIMEZERO value'
               call fcerr(errstr)
               goto 999
            endif
         else
            if (ftstatus .ne. 0) then
               errstr = 'error getting TIMEZERI value'
               call fcerr(errstr)
               goto 999
            endif
            call ftgkyd(iunit,'TIMEZERF',timezf,comment,ftstatus)
            if (ftstatus .ne. 0) then
               errstr = 'error getting TIMEZERF value'
               call fcerr(errstr)
               goto 999
            endif
            tzero = timezi + timezf
         endif

C      Reckon TZERO value to include MJDREF, converting TZERO to days
c      if necessary
         if (spanunit .eq. 's') then
            tzero = mjdref + tzero / 86400.
         else
            tzero = mjdref + tzero
         endif

c      Determine column numbers of TIME, RATE and ERROR columns, if present
c        For RATE, get the units on the RATE column as well
         do j = 1,ncol
            if (TTYPE(j) .eq. timename) then
               timecol = .true.
               timepos = j
               timeunit = TUNIT(j)
            end if
            if (TTYPE(j) .eq. ratename) then
               ratecol = .true.
               ratepos = j
               runit = TUNIT(j)
            endif
            if (TTYPE(j) .eq. errname) then
               errcol = .true.
               errpos = j
            endif
         end do
c MJT 15July96 (g77/linux) change to .eqv./.neqv. from .eq./.ne.
c     16Aug  -- treating logicals in a more sensible way
         if (.not. ratecol) then
            abort = .true.
            errstr = 'Requested Intensity column not found'
            call fcerr(errstr)
            go to 999
         endif
c MJT 15July96 (g77/linux) change to .eqv./.neqv. from .eq./.ne.
c     16Aug  -- treating logicals in a more sensible way
         if ((errname .ne. 'ERROR') .and. (.not. errcol)) then
            abort = .true.
            errstr = 'Requested Error column not found'
            call fcerr(errstr)
            go to 999
         endif

c      Read file with TIME column present
         if (timecol) then

c      Reckon bin size into units of the time array
           if (spanunit .ne. timeunit) then
              if (timeunit .eq. 's') timedel = 86400. * timedel
              if (timeunit .eq. 'd') timedel = timedel / 86400.
           endif

c      read the time, counts and error
c      (assume TIME in column 1, RATE in column 2, and ERROR (if present)
c           in column 3)
           call ftgcvd(iunit,timepos,frow,felem,nrows,nlvald,time,
     +              anynull,ftstatus)
           call ftgcve(iunit,ratepos,frow,felem,nrows,nlvale,rate,
     +              anynull,ftstatus)
           if (errcol) call ftgcve(iunit,errpos,frow,felem,nrows,
     +        nlvale,error,anynull,ftstatus)


c      Read FITS file having implicit TIME column
        else

c      Construct the Time array
         do j = 1, nrows
            time(j) = (j-1) * timedel
         end do
         timeunit = spanunit

c      Read the RATE and ERROR columns
c      (Assume RATE is in column 1, ERROR (if present) in col 2)
         call ftgcve(iunit,ratepos,frow,felem,nrows,nlvale,rate,anynull,
     &        ftstatus)
         if (errcol)
     &        call ftgcve(iunit,errpos,frow,felem,nrows,nlvale,error,
     &        anynull,ftstatus)
      endif


c      Assign absolute min & max times
      if (timeunit .eq. 's') then
         mintime = tzero + time(felem) / 86400.
         maxtime = tzero + time(nrows) / 86400.
      else
         mintime = tzero + time(felem)
         maxtime = tzero + time(nrows)
      endif

c  Exit subroutine
999	continue
	if (ftstatus .ne. 0) call fcerrm(ftstatus)

1000	format(A34,I3)
	return
	end
