      SUBROUTINE rdspec(string, iparse, iunit, dset_new, dset, mlabel)
c        reads spectra from a SF history file created by XSPEC
c        fwj haberl march 1988
      INCLUDE '../inc/historysf.inc'
      INCLUDE 'xplot.inc'
      INCLUDE '../inc/writebuf.inc'
      INTEGER*4 iparse, iunit, dset_new, newspec, dset
      CHARACTER*(*) string, mlabel
      INTEGER*4 infoar(4), nsub, ierrsf, index, lenact, nstrd, maxsize,
     &          nfree
      INTEGER*4 ibeg, iend, iflag, idelim, icom, len, ij, ix, iy, kbin,
     &          mbin, nset, iset
      REAL*4 xsavl, xsavh, chi, chir
      CHARACTER*12 pkgtyp
      CHARACTER*80 modstrg, bestmod
      LOGICAL*4 qdone, qfit, qskip
      PARAMETER (maxsize=max(len_his_data_files,len_his_fit,
     &           len_his_model,len_his_error,len_his_equiv_width))
      STRUCTURE /all_structure/
         UNION
            MAP
               BYTE buffer(maxsize)
            ENDMAP
            MAP
               RECORD /his_spectrum/spec
            ENDMAP
            MAP
               RECORD /his_fit/hfit
            ENDMAP
            MAP
               RECORD /his_channels/hchan
            ENDMAP
            MAP
               RECORD /his_param_range/parran
            ENDMAP
         ENDUNION
      ENDSTRUCTURE
      RECORD /his_spectrum_sub/spsub
      RECORD /his_model/hmod
      RECORD /all_structure/xxx
c entry rdmod
      RECORD /his_photon/pmod
      RECORD /his_photon_sub/pmodsub

      CALL xgtarg(string, iparse, ibeg, iend, qskip, iflag, idelim)
      IF (iflag.LT.0) THEN
         GOTO 900
      ENDIF
      IF (iend.GE.ibeg) THEN
         CALL xmatch(string(ibeg:iend), 'bestfit', 1, icom)
      ELSE
         icom = 0
      ENDIF
      IF (icom.EQ.1) THEN
         chi = 1.e30
         qfit = .FALSE.
         DO WHILE (.NOT.qfit)
            qdone = .FALSE.
            DO WHILE (.NOT.qdone)
               index = 0
               len = maxsize
               pkgtyp = ' '
               CALL nxpksf(iunit, pkgtyp, index, nsub, infoar,
     &                     xxx.buffer, len, .FALSE., ierrsf)
               IF (ierrsf.EQ.7 .OR. pkgtyp.EQ.'data files') THEN
                  qdone = .TRUE.
                  qfit = .TRUE.
               ENDIF
               IF (pkgtyp.EQ.'fit') THEN
                  qdone = .TRUE.
               ELSEIF (pkgtyp.EQ.'channels') THEN
                  len = len_his_channels
                  kbin = 0
                  DO WHILE (len.NE.0)
                     CALL rsubsf(iunit, xxx, len, ierrsf)
                     IF (len.NE.0) kbin = kbin + xxx.hchan.nset
                  ENDDO
               ELSEIF (pkgtyp.EQ.'model') THEN
                  len = len_his_model
                  modstrg = ' '
                  DO ij = 1, nsub
                     CALL rsubsf(iunit, hmod, len, ierrsf)
                     modstrg = modstrg(1:lenact(modstrg))
     &                         //' '//hmod.mname
                  ENDDO
                  pkgtyp = 'param def'
                  len = len_his_param_def
                  CALL nxpksf(iunit, pkgtyp, index, nsub, infoar, xxx,
     &                        len, .FALSE., ierrsf)
                  IF (ierrsf.NE.0) THEN
                     WRITE (*, '(a)') ' ERROR: No parameter definition '
     &                                //'package after model package'
                     GOTO 900
                  ENDIF
                  pkgtyp = 'param range'
                  CALL nxpksf(iunit, pkgtyp, index, nsub, infoar, xxx,
     &                        len, .FALSE., ierrsf)
                  IF (ierrsf.NE.0) THEN
                     WRITE (*, '(a)') ' ERROR: No parameter range '//
     &                                'package after model package'
                     GOTO 900
                  ENDIF
                  nfree = nsub
                  len = len_his_param_range
                  DO ij = 1, nsub
                     CALL rsubsf(iunit, xxx, len, ierrsf)
                     IF (xxx.parran.delta.LE.0.) THEN
                        nfree = nfree - 1
                     ENDIF
                  ENDDO
               ENDIF
            ENDDO
            IF (pkgtyp.EQ.'fit') THEN
               chir = xxx.hfit.chisqr/(kbin-nfree)
               WRITE (wrtstr, '(a)') indent//modstrg
               CALL xwrite(wrtstr, 10)
               WRITE (wrtstr, '(a,g15.5,a,g15.5,a)')
     &                 indent//'Chi square (reduced) of fit ',
     &                xxx.hfit.chisqr, ' (', chir, ' )'
               CALL xwrite(wrtstr, 10)
               IF (chir.LT.chi) THEN
c                 better fit found
                  bestmod = modstrg
                  mlabel = bestmod
                  chi = chir
                  index = 0
                  len = len_his_spectrum
                  pkgtyp = 'spectrum'
                  CALL nxpksf(iunit, pkgtyp, index, nsub, infoar, xxx,
     &                        len, .FALSE., ierrsf)
                  qdone = ierrsf.NE.0
                  IF (.NOT.qdone) THEN
                     len = len_his_spectrum_sub
                     xsavl = -2.
                     xsavh = -1.
                     newspec = 0
                     iset = dset + dset_new + newspec + 1
                     mbin = nbine(iset-1)
                     nbins(iset) = mbin + 1
                     DO ij = 1, xxx.spec.nbin
                        CALL rsubsf(iunit, spsub, len, ierrsf)
                        IF (spsub.elow.gt.1.0001*xsavh .OR.
     &                      spsub.elow.lt.0.9999*xsavh) THEN
                           IF (newspec.GT.0) THEN
                              nbine(iset) = mbin
                              nbins(iset+1) = mbin + 1
                              xlab(iset) = 'Channel Energy ('
     &                          //xxx.spec.eunits
     &                          (:lenact(xxx.spec.eunits))//')'
                              IF (xxx.spec.sunits.EQ.
     &                         'Counts/cm**2 sec keV') THEN
                                ylab(iset) = 'Counts cm'//char(92)//
     &                           'u-2'//char(92)//'d sec'//char(92)//
     &                           'u-1'//char(92)//'d keV'//char(92)//
     &                           'u-1'//char(92)//'d'
                              ELSE
                                ylab(iset) = xxx.spec.sunits
                              ENDIF
                           ENDIF
                           newspec = newspec + 1
                        ENDIF
                        mbin = mbin + 1
                        iset = dset + dset_new + newspec
                        datax(mbin) = 0.5*(spsub.elow+spsub.ehi)
                        errlx(mbin) = spsub.elow
                        errux(mbin) = spsub.ehi
                        datay(mbin) = spsub.phaobs
                        errly(mbin) = spsub.phaobs - spsub.phaerr
                        erruy(mbin) = spsub.phaobs + spsub.phaerr
                        modely(mbin) = spsub.phamod
                        xsavh = spsub.ehi
                        xsavl = spsub.elow
                     ENDDO
                     nbine(iset) = mbin
                     xlab(iset) = 'Channel Energy ('//xxx.spec.eunits
     &                          (:lenact(xxx.spec.eunits))//')'
                     IF (xxx.spec.sunits.EQ.'Counts/cm**2 sec keV') THEN
                        ylab(iset) = 'Counts cm'//char(92)//'u-2'//
     &                    char(92)//'d sec'//char(92)//'u-1'//
     &                    char(92)//'d keV'//char(92)//'u-1'//
     &                    char(92)//'d'
                     ELSE
                      ylab(iset) = xxx.spec.sunits
                     ENDIF
                  ENDIF
               ENDIF
            ENDIF
         ENDDO
         dset_new = dset_new + newspec
         WRITE (wrtstr, '(a,g15.5)') indent//
     &                               'Best fit with reduced chi square '
     &                               , chi
         CALL xwrite(wrtstr, 5)
         IF (lenact(bestmod).GT.0) THEN
            WRITE (wrtstr, '(a)') indent//'for model '//bestmod
         ELSE
            WRITE (wrtstr, '(a)') indent//
     &                          'no info on model given in history file'
         ENDIF
         CALL xwrite(wrtstr, 5)
      ELSE
         qdone = .FALSE.
         DO WHILE (.NOT.qdone)
            pkgtyp = ' '
            index = 0
            len = maxsize
            CALL nxpksf(iunit, pkgtyp, index, nsub, infoar, xxx, len,
     &                  .FALSE., ierrsf)
            qdone = ierrsf.EQ.7
            newspec = 0
            IF (.NOT.qdone) THEN
               IF (pkgtyp.EQ.'model') THEN
                  len = len_his_model
                  modstrg = ' '
                  DO ij = 1, nsub
                     CALL rsubsf(iunit, hmod, len, ierrsf)
                     modstrg = modstrg(1:lenact(modstrg))
     &                         //' '//hmod.mname
                  ENDDO
               ELSEIF (pkgtyp.EQ.'spectrum') THEN
                  len = len_his_spectrum_sub
                  xsavl = -2.
                  xsavh = -1.
                  iset = dset + dset_new + 1
                  mbin = nbine(iset-1)
                  nbins(iset) = mbin + 1
                  DO ij = 1, xxx.spec.nbin
                     CALL rsubsf(iunit, spsub, len, ierrsf)
                     IF (spsub.elow.gt.1.0001*xsavh .OR.
     &                   spsub.elow.lt.0.9999*xsavh) THEN
                        IF (newspec.GT.0) THEN
                           nbine(iset) = mbin
                           nbins(iset+1) = mbin + 1
                        ENDIF
                        dset_new = dset_new + 1
                        newspec = newspec + 1
                     ENDIF
                     mbin = mbin + 1
                     iset = dset + dset_new
                     datax(mbin) = 0.5*(spsub.elow+spsub.ehi)
                     errlx(mbin) = spsub.elow
                     errux(mbin) = spsub.ehi
                     datay(mbin) = spsub.phaobs
                     errly(mbin) = spsub.phaobs - spsub.phaerr
                     erruy(mbin) = spsub.phaobs + spsub.phaerr
                     modely(mbin) = spsub.phamod
                     xsavh = spsub.ehi
                     xsavl = spsub.elow
                  ENDDO
                  nbine(iset) = mbin
                  qdone = .TRUE.
               ENDIF
            ENDIF
         ENDDO
         IF (newspec.EQ.0) THEN
            WRITE (wrtstr, '(a)') ' No more spectrum in history file'
            CALL xwrite(wrtstr, 15)
            mlabel = ' '
            GOTO 900
         ELSE
            mlabel = modstrg
         ENDIF
      ENDIF
      DO ij = iset - newspec + 1, iset
         mbin = nbine(ij) - nbins(ij) + 1
         WRITE (wrtstr, '(2a,i3,a,i4,a)') indent, 'Data set ', ij, 
     &                                    ';', mbin, ' bins'
         CALL xwrite(wrtstr, 10)
         xlab(ij) = 'Channel Energy ('//xxx.spec.eunits
     &              (:lenact(xxx.spec.eunits))//')'
         IF (xxx.spec.sunits.EQ.'Counts/cm**2 sec keV') THEN
            ylab(ij) = 'Counts cm'//char(92)//'u-2'//char(92)//
     &      'd sec'//char(92)//'u-1'//char(92)//'d keV'//char(92)//
     &      'u-1'//char(92)//'d'
         ELSE
            ylab(ij) = xxx.spec.sunits
         ENDIF
      ENDDO
 900  CONTINUE
      RETURN
 
      ENTRY rdmod(iunit, dset_new, dset)
c		read photon model spectra from a SF history file created
c		by XSPEC dump command
c	fwj haberl   15-APR-1989
 
      pkgtyp = 'photon'
      index = 0
      len = len_his_photon
      CALL nxpksf(iunit, pkgtyp, index, nsub, infoar, pmod, len, .TRUE.,
     &            ierrsf)
      qdone = ierrsf.EQ.7
      IF (.NOT.qdone) THEN
         nset = dset + dset_new
         iset = nset + 1
         len = len_his_photon_sub
         nbins(iset) = nbine(nset) + 1
         nbine(iset) = nbine(nset) + pmod.nerg
         DO ij = nbins(iset), nbine(iset)
            CALL rsubsf(iunit, pmodsub, len, ierrsf)
            datax(ij) = 0.5*(pmodsub.elow+pmodsub.ehi)
            errlx(ij) = pmodsub.elow
            errux(ij) = pmodsub.ehi
            datay(ij) = pmodsub.photflux
            errly(ij) = pmodsub.photflux
            erruy(ij) = pmodsub.photflux
            modely(ij) = pmodsub.photflux
         ENDDO
      ELSE
         WRITE (*, '(a)')
     &                  ' No more photon model spectrum in history file'
         GOTO 910
      ENDIF
      WRITE (wrtstr, '(2a,i3,a,i4,a)') indent, 'Data set ', iset, ';',
     &                                 pmod.nerg, ' energy ranges'
      dset_new = dset_new + 1
      xlab(dset+dset_new) = 
     &        'Energy ('//pmod.eunits(:lenact(pmod.eunits))//')'
      IF (pmod.munits.EQ.'Photons/cm**2 s keV') THEN
         ylab(dset+dset_new) = 'Photons cm'//char(92)//'u-2'//
     &     char(92)//'d s'//char(92)//'u-1'//char(92)//'d keV'//
     &     char(92)//'u-1'//char(92)//'d'
      ELSE
         ylab(dset+dset_new) = pmod.munits
      ENDIF
 910  CONTINUE
      RETURN

      END
