      SUBROUTINE rdfit(iunit, string, iparse, dset, dset_new)
c     to read fit packages from history file
c     fwj haberl march 1988
      INCLUDE 'xplot.inc'
      INCLUDE '../inc/phasf.inc'
      INCLUDE '../inc/historysf.inc'
      INTEGER*4 iunit, dset, dset_new, dset_add, iset
      CHARACTER*(*) string
      LOGICAL*4 qtime, qdone, qskip, qpar
      INTEGER*4 binrng
      PARAMETER (binrng=10)
      REAL*4 sum(binrng), sumer(binrng), dist, rdist, edelta, edeltasq
      REAL*8 sttime, endtime, sttime_s, endtime_s, t0
      INTEGER*4 key, year, day, hour, minute
      REAL*4 seconds
      INTEGER*4 timearray(5)
      INTEGER*4 fl_pkgx, fl_pkgy, nofp, eqw_pkgx, eqw_pkgy,
     &          ibin(2, binrng)
      INTEGER*4 ierrsf, maxsize, lenbuf, index, nsub, infoar(4), nstrd
      INTEGER*4 lenact, len, iparse, ibeg, iend, iflag, idelim, nret
      CHARACTER*12 pkgtyp
      CHARACTER*8 mname(40), pname(100)
      INTEGER*4 modind(100)
      INTEGER*4 ncom
      PARAMETER (ncom=11)
      CHARACTER*10 subcom(ncom)
      INTEGER*4 isub, is, ix, ipar(2), nbin
      INTEGER*4 iparx, ipary(binrng), icomx, icomy, irng, ip
      INTEGER*4 n_data, n_fit, n_error, n_flux
      INTEGER*4 n_eqw, n_spec, n_add
      PARAMETER (maxsize=max(len_his_data_files,len_his_fit,
     &           len_his_error,len_his_equiv_width))
      STRUCTURE /all_structure/
         UNION
            MAP
               BYTE buffer(maxsize)
            ENDMAP
            MAP
               RECORD /his_data_files/hisd
            ENDMAP
            MAP
               RECORD /his_fit/hisf
            ENDMAP
            MAP
               RECORD /his_error/hise
            ENDMAP
            MAP
               RECORD /his_flux/flu
            ENDMAP
            MAP
               RECORD /his_equiv_width/ewid
            ENDMAP
            MAP
               RECORD /his_spectrum/spec
            ENDMAP
         ENDUNION
      ENDSTRUCTURE
      RECORD /his_data_files/datfil
      RECORD /his_fit_sub1/hisf1
      RECORD /his_fit_sub2/hisf2
      RECORD /his_model/hmod
      RECORD /his_param_def/hpar
      RECORD /his_spectrum_sub/spsub
      RECORD /all_structure/xxx
      DATA subcom/'jdate', 'day_year', 'flux', 'chi_squa', 'eqwidth',
     &     'phacts', 'modcts', 'luminosity', 'phflux', 'seconds',
     &     'phase'/
      DATA iparx/1/, ipary/binrng*1/, icomx/1/, icomy/-1/
      DATA fl_pkgx/1/, fl_pkgy/1/, nofp/1/, ibin/binrng*1, binrng*1/
      DATA eqw_pkgx/1/, eqw_pkgy/1/
 
      CALL xgtarg(string, iparse, ibeg, iend, qskip, iflag, idelim)
      IF (iflag.LT.0) THEN
         GOTO 999
      ENDIF
      IF (iend.GE.ibeg) THEN
         CALL xmatch(string(ibeg:iend), subcom, ncom, icomx)
      ENDIF
      dset_add = 1
      iset = dset + dset_new + dset_add
      IF (icomx.EQ.1) THEN
c jdate
         xlab(iset) = 'Julian Date - 2440000'
      ELSEIF (icomx.EQ.2) THEN
c day_year
         CALL xgtint(string, iparse, 1, 'year', 1, year, nret, iflag,
     &               idelim)
         WRITE (xlab(iset), '(a,i4)') 'Day ', year
      ELSEIF (icomx.EQ.10) THEN
c seconds
         CALL xgtint(string, iparse, 1, 'year', 1, year, nret, iflag,
     &               idelim)
         CALL xgtint(string, iparse, 1, 'day', 1, day, nret, iflag,
     &               idelim)
         CALL xgtint(string, iparse, 1, 'hour', 1, hour, nret, iflag,
     &               idelim)
         CALL xgtint(string, iparse, 1, 'minute', 1, minute, nret,
     &               iflag, idelim)
         CALL xgtrl4(string, iparse, 1, 'seconds', 1, seconds, nret,
     &               iflag, idelim)
         timearray(1) = year
         timearray(2) = day
         timearray(3) = hour
         timearray(4) = minute
         timearray(5) = seconds
         CALL timak(timearray, key)
         t0 = dble(key) + dble(seconds) - dble(timearray(5))
         IF (dset_new.EQ.0) THEN
            WRITE (xlab(iset), '(a,i4,a,i3,a,i2,a,i2,a,f6.3)')
     &              'Seconds after ', year, ' ', day, ' ', hour, ':',
     &             minute, ':', seconds
         ELSE
            WRITE (xlab(iset), '(a)') 'Seconds'
         ENDIF
      ELSEIF (icomx.EQ.11) THEN
c phase from selectors package in PHA files
         xlab(iset) = 'Phase'
      ELSEIF (icomx.EQ.3) THEN
c flux
         CALL xgtint(string, iparse, 1, 'flux package no', 1, fl_pkgx,
     &               nret, iflag, idelim)
         xlab(iset) = 'Flux (ergs sec'//char(92)//'u-1'//char(92)//
     &                'd cm'//char(92)//'u-2'//char(92)//'d)'
      ELSEIF (icomx.EQ.4) THEN
c chi_squa
         CALL xgtint(string, iparse, 1, 'no free param', 1, nofp, nret,
     &               iflag, idelim)
         xlab(iset) = char(92)//'gx'//char(92)//'u2'//char(92)//
     &                'd'//char(92)//'dr'//char(92)//'u'
      ELSEIF (icomx.EQ.5) THEN
c equivalent width
         CALL xgtint(string, iparse, 1, 'equiv. width package no', 1,
     &               eqw_pkgx, nret, iflag, idelim)
         WRITE (xlab(iset), '(a)') 'Equivalent width (keV)'
      ELSEIF (icomx.EQ.6) THEN
c PHA counts
         CALL xgtint(string, iparse, 2*binrng, '<11 start, end bins', 1,
     &               ibin, nret, iflag, idelim)
         dset_add = nret/2
      ELSEIF (icomx.EQ.7) THEN
c model counts
         CALL xgtint(string, iparse, 2*binrng, '<11 start, end bins', 1,
     &               ibin, nret, iflag, idelim)
         dset_add = nret/2
      ELSEIF (icomx.EQ.8) THEN
c luminosity
         CALL xgtint(string, iparse, 1, 'flux package no', 1, fl_pkgx,
     &               nret, iflag, idelim)
         CALL xgtrl4(string, iparse, 1, 'distance in kpc', 1, dist,
     &               nret, iflag, idelim)
         rdist = dist*3.0856E4
         rdist = rdist*rdist*4.*3.1415
         xlab(iset) = 'Luminosity (10'//char(92)//'u34'//char(92)//
     &                'd ergs sec'//char(92)//'u-1'//char(92)//'d)'
      ELSEIF (icomx.EQ.9) THEN
c phflux
         CALL xgtint(string, iparse, 1, 'flux package no', 1, fl_pkgx,
     &               nret, iflag, idelim)
         xlab(iset) = 'Flux (photons sec'//char(92)//'u-1'//char(92)//
     &                'd cm'//char(92)//'u-2'//char(92)//'d)'
      ELSE
c parameter no
         IF (iend.GE.ibeg) THEN
            iparse = ibeg - 1
            CALL xgtint(string, iparse, 1, 'fit par', 1, iparx, nret,
     &                  iflag, idelim)
            if(iflag.eq.0)then
               icomx = -1
            ELSE
               icomx = 0
            ENDIF
         ENDIF
      ENDIF
      CALL xgtarg(string, iparse, ibeg, iend, qskip, iflag, idelim)
      IF (iflag.LT.0) THEN
         GOTO 999
      ENDIF
      IF (iend.GE.ibeg) CALL xmatch(string(ibeg:iend), subcom, ncom,
     &                              icomy)
      IF (icomy.EQ.1) THEN
c jdate
         ylab(iset) = 'Julian Date - 2440000'
      ELSEIF (icomy.EQ.2) THEN
c day_year
         CALL xgtint(string, iparse, 1, 'year', 1, year, nret, iflag,
     &               idelim)
         WRITE (ylab(iset), '(a,i4)') 'Day ', year
      ELSEIF (icomy.EQ.10) THEN
c seconds
         CALL xgtint(string, iparse, 1, 'year', 1, year, nret, iflag,
     &               idelim)
         CALL xgtint(string, iparse, 1, 'day', 1, day, nret, iflag,
     &               idelim)
         CALL xgtint(string, iparse, 1, 'hour', 1, hour, nret, iflag,
     &               idelim)
         CALL xgtint(string, iparse, 1, 'minute', 1, minute, nret,
     &               iflag, idelim)
         CALL xgtrl4(string, iparse, 1, 'seconds', 1, seconds, nret,
     &               iflag, idelim)
         timearray(1) = year
         timearray(2) = day
         timearray(3) = hour
         timearray(4) = minute
         timearray(5) = seconds
         CALL timak(timearray, key)
         t0 = key + seconds - timearray(5)
         IF (dset_new.EQ.0) THEN
            WRITE (ylab(iset), '(a,i4,a,i3,a,i2,a,i2,a,f6.3)')
     &              'Seconds after ', year, ' ', day, ' ', hour, ':',
     &             minute, ':', seconds
         ELSE
            WRITE (ylab(iset), '(a)') 'Seconds'
         ENDIF
      ELSEIF (icomy.EQ.11) THEN
c phase from selectors package in PHA files
         ylab(iset) = 'Phase'
      ELSEIF (icomy.EQ.3) THEN
c flux
         CALL xgtint(string, iparse, 1, 'flux package no', 1, fl_pkgy,
     &               nret, iflag, idelim)
         ylab(iset) = 'Flux (ergs sec'//char(92)//'u-1'//char(92)//
     &                'd cm'//char(92)//'u-2'//char(92)//'d)'
      ELSEIF (icomy.EQ.4) THEN
c chi_squa
         CALL xgtint(string, iparse, 1, 'no free param', 1, nofp, nret,
     &               iflag, idelim)
         ylab(iset) = char(92)//'gx'//char(92)//'u2'//char(92)//'d'//
     &                char(92)//'dr'//char(92)//'u'
      ELSEIF (icomy.EQ.5) THEN
c equivalent width
         CALL xgtint(string, iparse, 1, 'equiv. width package no', 1,
     &               eqw_pkgy, nret, iflag, idelim)
         WRITE (ylab(iset), '(a)') 'Equivalent width (keV)'
      ELSEIF (icomy.EQ.6) THEN
c PHA counts
         CALL xgtint(string, iparse, 2*binrng, '<11 start, end bin', 1,
     &               ibin, nret, iflag, idelim)
         dset_add = nret/2
      ELSEIF (icomy.EQ.7) THEN
c model counts
         CALL xgtint(string, iparse, 2*binrng, '<11 start, end bin', 1,
     &               ibin, nret, iflag, idelim)
         dset_add = nret/2
      ELSEIF (icomy.EQ.8) THEN
c luminosity
         CALL xgtint(string, iparse, 1, 'flux package no', 1, fl_pkgy,
     &               nret, iflag, idelim)
         CALL xgtrl4(string, iparse, 1, 'distance in kpc', 1, dist,
     &               nret, iflag, idelim)
         rdist = dist*3.0856E4
         rdist = rdist*rdist*4.*3.1415
         ylab(iset) = 'Luminosity (10'//char(92)//'u34'//char(92)//
     &                'd ergs sec'//char(92)//'u-1'//char(92)//'d)'
      ELSEIF (icomy.EQ.9) THEN
c phflux
         CALL xgtint(string, iparse, 1, 'flux package no', 1, fl_pkgy,
     &               nret, iflag, idelim)
         ylab(iset) = 'Flux (photons sec'//char(92)//'u-1'//char(92)//
     &                'd cm'//char(92)//'u-2'//char(92)//'d)'
      ELSE
c parameter no
         IF (iend.GE.ibeg) THEN
            iparse = ibeg - 1
            CALL xgtint(string, iparse, binrng, '<11 fit parameters', 
     &                  1, ipary, nret, iflag, idelim)
            if(iflag.gt.0)then
               dset_add = min(nret,binrng)
               icomy = -dset_add
            ELSE
               icomy = 0
            ENDIF
         ENDIF
      ENDIF
      qdone = .FALSE.
      n_data = nbine(dset+dset_new)
      n_fit = -1
      n_error = -1
      n_flux = -1
      n_spec = -1
      n_eqw = -1
      is = dset + dset_new + 1
      qpar = .FALSE.
      DO WHILE (.NOT.qdone)
         ierrsf = 0
         lenbuf = maxsize
         pkgtyp = ' '
         index = 0
         CALL nxpksf(iunit, pkgtyp, index, nsub, infoar, xxx.buffer,
     &               lenbuf, .FALSE., ierrsf)
         qdone = ierrsf.NE.0
         IF (qdone) THEN
            GOTO 900
         ENDIF
         IF (pkgtyp.EQ.'data files') THEN
            WRITE (*, '(4a)') indent, 'package ',
     &                        pkgtyp(:lenact(pkgtyp)), ' found'
            n_add = 1
            IF (n_fit.EQ.0) THEN
c no fit package
               IF (icomx.LT.0 .OR. icomy.LT.0) n_add = 0
            ENDIF
            IF (n_flux.EQ.0) THEN
c no flux package
               IF (icomx.EQ.3 .OR. icomy.EQ.3) n_add = 0
               IF (icomx.EQ.8 .OR. icomy.EQ.8) n_add = 0
               IF (icomx.EQ.9 .OR. icomy.EQ.9) n_add = 0
            ENDIF
            IF (n_eqw.EQ.0) THEN
c no equivalent width package
               IF (icomx.EQ.5 .OR. icomy.EQ.5) n_add = 0
            ENDIF
            IF (n_spec.EQ.0) THEN
c no spectrum package
               IF (icomx.EQ.6 .OR. icomy.EQ.6 .OR. icomx.EQ.7 .OR.
     &             icomy.EQ.7) n_add = 0
            ENDIF
            n_data = n_data + n_add
            n_fit = 0
            n_error = 0
            n_flux = 0
            n_eqw = 0
            n_spec = 0
            nstrd = 1
            DO isub = 1, nsub
               CALL rstrsf(iunit, datfil.charbuf, 1, nstrd, ierrsf)
               WRITE (*, '(2a)') indent,
     &                           datfil.fname(:lenact(datfil.fname))
               IF (icomx.EQ.1 .OR. icomx.EQ.2 .OR. icomx.EQ.10 .OR.
     &             icomy.EQ.1 .OR. icomy.EQ.2 .OR. icomy.EQ.10 .OR.
     &             icomx.EQ.11 .OR. icomy.EQ.11) THEN
                  IF (icomx.EQ.11 .OR. icomy.EQ.11) THEN
                     CALL phaphas(datfil.fname, sttime, endtime)
                  ELSEIF (icomx.EQ.10 .OR. icomy.EQ.10) THEN
                     CALL phatim(datfil.fname, sttime, endtime, .FALSE.)
                  ELSE
                     CALL phatim(datfil.fname, sttime, endtime, .TRUE.)
                  ENDIF
                  IF (isub.GT.1) THEN
                     IF (endtime.GT.endtime_s) THEN
                        WRITE (*, '(3a)') indent,
     &                'WARNING: data files contain different end times.'
     &                , ' Latest time will be used.'
                        endtime_s = endtime
                     ENDIF
                     IF (sttime.LT.sttime_s) THEN
                        WRITE (*, '(3a)') indent,
     &              'WARNING: data files contain different start times.'
     &              , ' Latest time will be used.'
                     ENDIF
                  ELSE
                     sttime_s = sttime
                     endtime_s = endtime
                  ENDIF
               ENDIF
            ENDDO
c Time from PHA files as jdate (1)
c                     or day number of given year (2)
c                     or seconds after reference time (10)
c Or phase from selctors package in PHA files (11)
            IF (icomx.EQ.11) THEN
               datax(n_data) = (sttime_s+endtime_s)*0.5D00
               errlx(n_data) = sttime_s
               errux(n_data) = endtime_s
            ELSEIF (icomy.EQ.11) THEN
               datay(n_data) = (sttime_s+endtime_s)*0.5D00
               errly(n_data) = sttime_s
               erruy(n_data) = endtime_s
            ELSEIF (icomx.EQ.10) THEN
               datax(n_data) = (sttime_s+endtime_s)*0.5D00 - t0
               errlx(n_data) = sttime_s - t0
               errux(n_data) = endtime_s - t0
            ELSEIF (icomy.EQ.10) THEN
               datay(n_data) = (sttime_s+endtime_s)*0.5D00 - t0
               errly(n_data) = sttime_s - t0
               erruy(n_data) = endtime_s - t0
            ELSEIF (icomx.EQ.1 .OR. icomx.EQ.2) THEN
               CALL timscl(icomx, t0, year)
               datax(n_data) = (sttime_s+endtime_s)*0.5D00 - t0
               errlx(n_data) = sttime_s - t0
               errux(n_data) = endtime_s - t0
            ELSEIF (icomy.EQ.1 .OR. icomy.EQ.2) THEN
               CALL timscl(icomy, t0, year)
               datay(n_data) = (sttime_s+endtime_s)*0.5D00 - t0
               errly(n_data) = sttime_s - t0
               erruy(n_data) = endtime_s - t0
            ENDIF
         ELSEIF (pkgtyp.EQ.'fit') THEN
            WRITE (*, '(4a)') indent, 'package ',
     &                        pkgtyp(:lenact(pkgtyp)), ' found'
            IF (n_fit.GE.1) WRITE (*, '(3a)') indent,
     &                     'more than one fit package for current data '
     &                     , 'file(s), last one will be used'
            IF (xxx.hisf.error.NE.0) THEN
               WRITE (*, '(2a,i3,a)') indent, 'Error flag ',
     &                                xxx.hisf.error,
     &                                ' at conclusion of fit'
            ELSE
               n_fit = n_fit + 1
               len = 4*xxx.hisf.nfit
               CALL rsubsf(iunit, hisf1, len, ierrsf)
               CALL rsubsf(iunit, hisf2, len, ierrsf)
               IF (icomx.LT.0) THEN
c parameter
                  datax(n_data) = hisf1.value(iparx)
                  IF (hisf2.esterr(iparx).LT.0.) hisf2.esterr(iparx)
     &                 = 0.
                  errlx(n_data) = datax(n_data) - hisf2.esterr(iparx)
                  IF (errlx(n_data).LT.0.) THEN
                     IF (qpar) THEN
                        IF (pname(iparx).EQ.'nH/10^22')THEN
                           errlx(n_data) = 0.
                        ENDIF
                     ENDIF
                  ENDIF
                  errux(n_data) = datax(n_data) + hisf2.esterr(iparx)
               ELSEIF (icomx.EQ.4) THEN
c chi_squa
                  datax(n_data) = xxx.hisf.chisqr/nofp
                  errlx(n_data) = datax(n_data)
                  errux(n_data) = datax(n_data)
               ENDIF
               IF (icomy.LT.0) THEN
c parameter
                  datay(n_data) = hisf1.value(ipary(1))
                  IF (hisf2.esterr(ipary(1)).LT.0.) 
     &                 hisf2.esterr(ipary(1)) = 0.
                  errly(n_data) = datay(n_data) - hisf2.esterr(ipary(1))
                  IF (errly(n_data).LT.0.) THEN
                     IF (qpar) THEN
                        IF (pname(ipary(1)).EQ.'nH/10^22')THEN
                           errly(n_data) = 0.
                        ENDIF
                     ENDIF
                  ENDIF
                  erruy(n_data) = datay(n_data) + hisf2.esterr(ipary(1))
c                 store other fit parameters in buf:
                  ip = n_data - nbine(dset+dset_new)
                  DO irng = 2, dset_add
                     is = (irng-1)*3 - 2
                     buf(ip, is) = hisf1.value(ipary(irng))
                     IF (hisf2.esterr(ipary(irng)).LT.0.) 
     &                    hisf2.esterr(ipary(irng)) = 0.
                     buf(ip, is+1) = buf(ip, is) - 
     &                               hisf2.esterr(ipary(irng))
                     IF (buf(ip, is+1).LT.0.) THEN
                        IF (qpar) THEN
                           IF (pname(ipary(irng)).EQ.'nH/10^22')THEN
                              buf(ip, is+1) = 0.
                           ENDIF
                        ENDIF
                     ENDIF
                     buf(ip, is+2) = buf(ip, is) + 
     &                               hisf2.esterr(ipary(irng))
                  ENDDO
               ELSEIF (icomy.EQ.4) THEN
c chi_squa
                  datay(n_data) = xxx.hisf.chisqr/nofp
                  errly(n_data) = datay(n_data)
                  erruy(n_data) = datay(n_data)
               ENDIF
            ENDIF
         ELSEIF (pkgtyp.EQ.'flux') THEN
            n_flux = n_flux + 1
            IF (icomx.EQ.3) THEN
c flux
               IF (n_flux.EQ.fl_pkgx) THEN
                  WRITE (*, '(a,i2,3a)') indent, n_flux, 'th package ',
     &                                   pkgtyp(:lenact(pkgtyp)),
     &                                   ' found'
                  datax(n_data) = xxx.flu.eflux
                  errlx(n_data) = xxx.flu.eflux*.99
                  errux(n_data) = xxx.flu.eflux*1.01
               ENDIF
            ELSEIF (icomx.EQ.9) THEN
c phflux
               IF (n_flux.EQ.fl_pkgx) THEN
                  WRITE (*, '(a,i2,3a)') indent, n_flux, 'th package ',
     &                                   pkgtyp(:lenact(pkgtyp)),
     &                                   ' found'
                  datax(n_data) = xxx.flu.pflux
                  errlx(n_data) = xxx.flu.pflux*.99
                  errux(n_data) = xxx.flu.pflux*1.01
               ENDIF
            ELSEIF (icomx.EQ.8) THEN
c luminosity
               IF (n_flux.EQ.fl_pkgx) THEN
                  WRITE (*, '(a,i2,3a)') indent, n_flux, 'th package ',
     &                                   pkgtyp(:lenact(pkgtyp)),
     &                                   ' found'
                  datax(n_data) = xxx.flu.eflux*rdist
                  errlx(n_data) = datax(n_data)*.99
                  errux(n_data) = datax(n_data)*1.01
               ENDIF
            ENDIF
            IF (icomy.EQ.3) THEN
c flux
               IF (n_flux.EQ.fl_pkgy) THEN
                  WRITE (*, '(a,i2,3a)') indent, n_flux, 'th package ',
     &                                   pkgtyp(:lenact(pkgtyp)),
     &                                   ' found'
                  datay(n_data) = xxx.flu.eflux
                  errly(n_data) = xxx.flu.eflux*.99
                  erruy(n_data) = xxx.flu.eflux*1.01
               ENDIF
            ELSEIF (icomy.EQ.9) THEN
c phflux
               IF (n_flux.EQ.fl_pkgy) THEN
                  WRITE (*, '(a,i2,3a)') indent, n_flux, 'th package ',
     &                                   pkgtyp(:lenact(pkgtyp)),
     &                                   ' found'
                  datay(n_data) = xxx.flu.pflux
                  errly(n_data) = xxx.flu.pflux*.99
                  erruy(n_data) = xxx.flu.pflux*1.01
               ENDIF
            ELSEIF (icomy.EQ.8) THEN
c luminosity
               IF (n_flux.EQ.fl_pkgy) THEN
                  WRITE (*, '(a,i2,3a)') indent, n_flux, 'th package ',
     &                                   pkgtyp(:lenact(pkgtyp)),
     &                                   ' found'
                  datay(n_data) = xxx.flu.eflux*rdist
                  errly(n_data) = datay(n_data)*.99
                  erruy(n_data) = datay(n_data)*1.01
               ENDIF
            ENDIF
         ELSEIF (pkgtyp.EQ.'equiv width') THEN
            n_eqw = n_eqw + 1
            IF (icomx.EQ.5) THEN
c equivalent width
               IF (n_eqw.EQ.eqw_pkgx) THEN
                  WRITE (*, '(a,i2,3a)') indent, n_eqw, 'th package ',
     &                                   pkgtyp(:lenact(pkgtyp)),
     &                                   ' found'
                  datax(n_data) = xxx.ewid.eqp
                  errlx(n_data) = xxx.ewid.eqp*.99
                  errux(n_data) = xxx.ewid.eqp*1.01
               ENDIF
            ENDIF
            IF (icomy.EQ.5) THEN
c equivalent width
               IF (n_eqw.EQ.eqw_pkgy) THEN
                  WRITE (*, '(a,i2,3a)') indent, n_eqw, 'th package ',
     &                                   pkgtyp(:lenact(pkgtyp)),
     &                                   ' found'
                  datay(n_data) = xxx.ewid.eqp
                  errly(n_data) = xxx.ewid.eqp*.99
                  erruy(n_data) = xxx.ewid.eqp*1.01
               ENDIF
            ENDIF
         ELSEIF (pkgtyp.EQ.'error') THEN
            WRITE (*, '(4a)') indent, 'package ',
     &                        pkgtyp(:lenact(pkgtyp)), ' found'
            IF (n_fit.EQ.0) THEN
               WRITE (*, '(2a)') indent,
     &                 'No useful fit package before this error package'
            ELSEIF (xxx.hise.errflag.NE.0) THEN
               WRITE (*, '(2a,i3,a)') indent, 'Error flag ',
     &                                xxx.hise.errflag, ' indicated'
            ELSE
               IF (xxx.hise.fitind.EQ.iparx) THEN
                  errlx(n_data) = xxx.hise.low
                  errux(n_data) = xxx.hise.high
               ELSEIF (xxx.hise.fitind.EQ.ipary(1)) THEN
                  errly(n_data) = xxx.hise.low
                  erruy(n_data) = xxx.hise.high
               ELSEIF (icomy.lt.-1) THEN
                  ip = n_data - nbine(dset+dset_new)
                  DO irng = 2, dset_add
                     IF (xxx.hise.fitind.EQ.ipary(irng)) THEN
                        is = (irng-1)*3 - 2
                        buf(ip, is+1) = xxx.hise.low
                        buf(ip, is+2) = xxx.hise.high
                     ENDIF
                  ENDDO
               ENDIF
            ENDIF
         ELSEIF (pkgtyp.EQ.'model') THEN
            WRITE (*, '(4a)') indent, 'package ',
     &                        pkgtyp(:lenact(pkgtyp)), ' found'
            len = len_his_model
            DO isub = 1, nsub
               CALL rsubsf(iunit, hmod, len, ierrsf)
               mname(isub) = hmod.mname
            ENDDO
         ELSEIF (pkgtyp.EQ.'param def') THEN
            WRITE (*, '(4a)') indent, 'package ',
     &                        pkgtyp(:lenact(pkgtyp)), ' found'
            len = len_his_param_def
            DO isub = 1, nsub
               CALL rsubsf(iunit, hpar, len, ierrsf)
               pname(isub) = hpar.pname
               modind(isub) = hpar.modind
            ENDDO
            qpar = .TRUE.
         ELSEIF (pkgtyp.EQ.'spectrum') THEN
c pha counts
            IF (icomx.EQ.6 .OR. icomy.EQ.6) THEN
               WRITE (*, '(4a)') indent, 'package ',
     &                           pkgtyp(:lenact(pkgtyp)), ' found'
               n_spec = n_spec + 1
               IF (n_spec.GT.1) WRITE (*, '(3a)') indent,
     &                'more than one spectrum package for current data '
     &                , 'file(s), last one will be used'
               len = len_his_spectrum_sub
               DO irng = 1, dset_add
                  sum(irng) = 0.
                  sumer(irng) = 0.
               ENDDO
               DO isub = 1, xxx.spec.nbin
                  CALL rsubsf(iunit, spsub, len, ierrsf)
                  DO irng = 1, dset_add
                     IF (isub.GE.ibin(1,irng) .AND. isub.LE.ibin(2,irng)
     &                   ) THEN
                        edelta = spsub.ehi - spsub.elow
                        edeltasq = edelta*edelta
                        sum(irng) = sum(irng) + spsub.phaobs*edelta
                        sumer(irng) = sumer(irng)
     &                           + spsub.phaerr*spsub.phaerr*edeltasq
                     ENDIF
                  ENDDO
               ENDDO
               DO irng = 1, dset_add
                  sumer(irng) = sqrt(sumer(irng))
               ENDDO
               IF (icomx.EQ.6) THEN
c pha counts
                  iset = dset + dset_new
                  IF (xxx.spec.sunits.EQ.'Counts/cm**2 sec keV' .or.
     &                xxx.spec.sunits.EQ.
     &                 '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')THEN
                     DO irng = 1, dset_add
                        iset = iset+1
                        xlab(iset) = 'Counts cm'//char(92)//'u-2'//
     &                               char(92)//'d sec'//char(92)//
     &                               'u-1'//char(92)//'d'
                     ENDDO
                  ELSE                  
                     DO irng = 1, dset_add
                        iset = iset+1
                        xlab(iset) = ' '
                     ENDDO
                  ENDIF
                  irng = 1
                  datax(n_data) = sum(irng)
                  errlx(n_data) = sum(irng) - sumer(irng)
                  errux(n_data) = sum(irng) + sumer(irng)
c                 store other channel ranges in buf:
                  ip = n_data - nbine(dset+dset_new)
                  DO irng = 2, dset_add
                     is = (irng-1)*3 - 2
                     buf(ip, is) = sum(irng)
                     buf(ip, is+1) = sum(irng) - sumer(irng)
                     buf(ip, is+2) = sum(irng) + sumer(irng)
                  ENDDO
               ENDIF
               IF (icomy.EQ.6) THEN
c pha counts
                  iset = dset + dset_new
                  IF (xxx.spec.sunits.EQ.'Counts/cm**2 sec keV' .or.
     &                xxx.spec.sunits.EQ.
     &                 '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')THEN
                     DO irng = 1, dset_add
                        iset = iset+1
                        ylab(iset) = 'Counts cm'//char(92)//'u-2'//
     &                               char(92)//'d sec'//char(92)//
     &                               'u-1'//char(92)//'d'
                     ENDDO
                  ELSE                  
                     DO irng = 1, dset_add
                        iset = iset+1
                        ylab(iset) = ' '
                     ENDDO
                  ENDIF
                  irng = 1
                  datay(n_data) = sum(irng)
                  errly(n_data) = sum(irng) - sumer(irng)
                  erruy(n_data) = sum(irng) + sumer(irng)
c                       store other channel ranges in buf:
                  ip = n_data - nbine(dset+dset_new)
                  DO irng = 2, dset_add
                     is = (irng-1)*3 - 2
                     buf(ip, is) = sum(irng)
                     buf(ip, is+1) = sum(irng) - sumer(irng)
                     buf(ip, is+2) = sum(irng) + sumer(irng)
                  ENDDO
               ENDIF
            ELSEIF (icomx.EQ.7 .OR. icomy.EQ.7) THEN
c model counts
               WRITE (*, '(4a)') indent, 'package ',
     &                           pkgtyp(:lenact(pkgtyp)), ' found'
               n_spec = n_spec + 1
               IF (n_spec.GT.1) WRITE (*, '(3a)') indent,
     &                'more than one spectrum package for current data '
     &                , 'file(s), last one will be used'
               len = len_his_spectrum_sub
               DO irng = 1, dset_add
                  sum(irng) = 0.
               ENDDO
               DO isub = 1, xxx.spec.nbin
                  CALL rsubsf(iunit, spsub, len, ierrsf)
                  DO irng = 1, dset_add
                     IF (isub.GE.ibin(1,irng) .AND. isub.LE.ibin(2,irng)
     &                   ) THEN
                        edelta = spsub.ehi - spsub.elow
                        sum(irng) = sum(irng) + spsub.phamod*edelta
                     ENDIF
                  ENDDO
               ENDDO
               IF (icomx.EQ.7) THEN
c model counts
                  iset = dset + dset_new
                  IF (xxx.spec.sunits.EQ.'Counts/cm**2 sec keV' .or.
     &                xxx.spec.sunits.EQ.
     &                 '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')THEN
                     DO irng = 1, dset_add
                        iset = iset+1
                        xlab(iset) = 'Counts cm'//char(92)//'u-2'//
     &                               char(92)//'d sec'//char(92)//
     &                               'u-1'//char(92)//'d'
                     ENDDO
                  ELSE                  
                     DO irng = 1, dset_add
                        iset = iset+1
                        xlab(iset) = ' '
                     ENDDO
                  ENDIF
                  irng = 1
                  datax(n_data) = sum(irng)
                  errlx(n_data) = sum(irng)
                  errux(n_data) = sum(irng)
c                 store other channel ranges in buf:
                  ip = n_data - nbine(dset+dset_new)
                  DO irng = 2, dset_add
                     is = (irng-1)*3 - 2
                     buf(ip, is) = sum(irng)
                     buf(ip, is+1) = sum(irng)
                     buf(ip, is+2) = sum(irng)
                  ENDDO
               ENDIF
               IF (icomy.EQ.7) THEN
c model counts
                  iset = dset + dset_new
                  IF (xxx.spec.sunits.EQ.'Counts/cm**2 sec keV' .or.
     &                xxx.spec.sunits.EQ.
     &                 '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')THEN
                     DO irng = 1, dset_add
                        iset = iset+1
                        ylab(iset) = 'Counts cm'//char(92)//'u-2'//
     &                               char(92)//'d sec'//char(92)//
     &                               'u-1'//char(92)//'d'
                     ENDDO
                  ELSE                  
                     DO irng = 1, dset_add
                        iset = iset+1
                        ylab(iset) = ' '
                     ENDDO
                  ENDIF
                  irng = 1
                  datay(n_data) = sum(irng)
                  errly(n_data) = sum(irng)
                  erruy(n_data) = sum(irng)
c                 store other channel ranges in buf:
                  ip = n_data - nbine(dset+dset_new)
                  DO irng = 2, dset_add
                     is = (irng-1)*3 - 2
                     buf(ip, is) = sum(irng)
                     buf(ip, is+1) = sum(irng)
                     buf(ip, is+2) = sum(irng)
                  ENDDO
               ENDIF
            ENDIF
         ENDIF
      ENDDO
 900  CONTINUE
      nbin = n_data - nbine(dset+dset_new)
      IF (nbin.LT.1) THEN
         GOTO 999
      ENDIF
      DO irng = 1, dset_add
         is = dset + dset_new + irng
         nbins(is) = nbine(is-1) + 1
         nbine(is) = nbine(is-1) + nbin
         WRITE (*, '(2a,i3,a,i4,a)') indent, 'Data set ', is, ';', nbin,
     &                               ' data points'
      ENDDO
      IF (icomx.EQ.6 .OR. icomx.EQ.7) THEN
         n_data = nbine(dset+dset_new+1)
         DO irng = 2, dset_add
            is = (irng-1)*3 - 2
            DO ip = 1, nbin
               n_data = n_data + 1
               datax(n_data) = buf(ip, is)
               errlx(n_data) = buf(ip, is+1)
               errux(n_data) = buf(ip, is+2)
               datay(n_data) = datay(n_data-nbin)
               errly(n_data) = errly(n_data-nbin)
               erruy(n_data) = erruy(n_data-nbin)
            ENDDO
         ENDDO
      ENDIF
      IF (icomy.EQ.6 .OR. icomy.EQ.7 .OR. icomy.LT.0) THEN
         n_data = nbine(dset+dset_new+1)
         DO irng = 2, dset_add
            is = (irng-1)*3 - 2
            DO ip = 1, nbin
               n_data = n_data + 1
               datay(n_data) = buf(ip, is)
               errly(n_data) = buf(ip, is+1)
               erruy(n_data) = buf(ip, is+2)
               datax(n_data) = datax(n_data-nbin)
               errlx(n_data) = errlx(n_data-nbin)
               errux(n_data) = errux(n_data-nbin)
            ENDDO
         ENDDO
      ENDIF
      IF (qpar) THEN
         IF (icomx.LT.0) THEN
            iset = dset + dset_new
            DO irng = 1, dset_add
               iset = iset+1
               xlab(iset) = mname(modind(iparx))
     &                      //' '//pname(iparx)
            ENDDO
            iset = dset + dset_new + 1
            DO irng = 2, dset_add
               iset = iset+1
               ylab(iset) = ylab(iset-1)
            ENDDO
         ENDIF
         IF (icomy.LT.0) THEN
            iset = dset + dset_new
            DO irng = 1, dset_add
               iset = iset+1
               ylab(iset) = mname(modind(ipary(irng)))
     &                      //' '//pname(ipary(irng))
            ENDDO
            iset = dset + dset_new + 1
            DO irng = 2, dset_add
               iset = iset+1
               xlab(iset) = xlab(iset-1)
            ENDDO
         ENDIF
      ELSE
         IF (icomx.LT.0) THEN
            iset = dset + dset_new
            DO irng = 1, dset_add
               iset = iset+1
               WRITE (xlab(iset), '(a,i2)') 'parameter ', iparx
            ENDDO
            iset = dset + dset_new + 1
            DO irng = 2, dset_add
               iset = iset+1
               ylab(iset) = ylab(iset-1)
            ENDDO
         ENDIF
         IF (icomy.LT.0) THEN
            iset = dset + dset_new
            DO irng = 1, dset_add
               iset = iset+1
               WRITE (ylab(iset), '(a,i2)') 'parameter ', ipary(irng)
            ENDDO
            iset = dset + dset_new
            DO irng = 2, dset_add + 1
               iset = iset+1
               xlab(iset) = xlab(iset-1)
            ENDDO
         ENDIF
      ENDIF
      dset_new = dset_new + dset_add
 999  CONTINUE
      RETURN
      END
