      SUBROUTINE addat(string, iparse, dset_new, dset)
c        adds two data sets
c        fwj haberl may 1988
      INCLUDE 'xplot.inc'
      INCLUDE '../inc/writebuf.inc'
      CHARACTER*(*) string
      INTEGER*4 iparse, dset, iset(2), dset_new, jset
      INTEGER*4 nret, iflag, idelim, i, i1, i2, mbin
      INTEGER*4 ncom
      PARAMETER (ncom=2)
      CHARACTER*10 subcom(ncom)
      INTEGER*4 ibeg, iend, icom
      REAL*4 e1, e2, const, cerrl, cerru
      LOGICAL*4 qskip
      DATA subcom/'constant', 'x'/
 
      CALL xgtint(string, iparse, 1, 'dataset 1', 1, iset(1), nret,
     &            iflag, idelim)
      IF (iflag.LT.0) THEN
         GOTO 900
      ENDIF
      jset = dset + dset_new
      IF (iset(1).LT.1 .OR. iset(1).GT.jset) THEN
         WRITE (*, '(a,i3)') ' Data set no ', iset(1), ' not available'
         GOTO 900
      ENDIF
      i1 = nbine(iset(1)) - nbins(iset(1)) + 1
      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), subcom, ncom, icom)
      ENDIF
      const = 0.
      cerrl = 0.
      cerru = 0.
      IF (icom.EQ.1) THEN
         CALL xgtrl4(string, iparse, 1, 'const', 1, const, nret, iflag,
     &               idelim)
         CALL xgtrl4(string, iparse, 1, 'error lo', 1, cerrl, nret,
     &               iflag, idelim)
         CALL xgtrl4(string, iparse, 1, 'error hi', 1, cerru, nret,
     &               iflag, idelim)
      ELSE
         IF (iend.GE.ibeg) THEN
            CALL xptarg(string, iparse, string(ibeg:iend), iflag, 
     &                  idelim)
            CALL xgtint(string, iparse, 1, 'dataset 2', 1, iset(2),
     &                  nret, iflag, idelim)
            IF (iset(2).LT.1 .OR. iset(2).GT.jset) THEN
               WRITE (*, '(a,i3)') ' Data set no ', iset(2),
     &                             ' not available'
               GOTO 900
            ELSE
               i2 = nbine(iset(2)) - nbins(iset(2)) + 1
               IF (i1.NE.i2) THEN
                  WRITE (*, '(a)')
     &                 ' The 2 data sets have different numbers of bins'
                  GOTO 900
               ENDIF
            ENDIF
         ENDIF
      ENDIF
      dset_new = dset_new + 1
      jset = dset + dset_new
      xlab(jset) = xlab(iset(1))
      nbins(jset) = nbine(jset-1) + 1
      nbine(jset) = nbine(jset-1) + i1
      IF (icom.EQ.1) THEN
         i1 = nbins(iset(1))
         DO i = nbins(jset), nbine(jset)
            datax(i) = datax(i1)
            errlx(i) = errlx(i1)
            errux(i) = errux(i1)
            datay(i) = datay(i1) + const
            e1 = datay(i1) - errly(i1)
            e2 = const - cerrl
            e1 = sqrt(e1*e1+e2*e2)
            errly(i) = datay(i) - e1
            e1 = erruy(i1) - datay(i1)
            e2 = cerru - const
            e1 = sqrt(e1*e1+e2*e2)
            erruy(i) = datay(i) + e1
            modely(i) = modely(i1) + const
            i1 = i1 + 1
         ENDDO
         WRITE (ylab(jset), '(a,i3,a,g15.4)') 
     &                  'data set', iset(1), ' + ', const
      ELSE
         i1 = nbins(iset(1))
         i2 = nbins(iset(2))
         DO i = nbins(jset), nbine(jset)
            IF (datax(i1).NE.datax(i2)) THEN
               WRITE (*, '(a)') ' Warning: data has different x-values'
            ENDIF
            datax(i) = datax(i1)
            errlx(i) = errlx(i1)
            errux(i) = errux(i1)
            datay(i) = datay(i1) + datay(i2)
            e1 = datay(i1) - errly(i1)
            e2 = datay(i2) - errly(i2)
            e1 = sqrt(e1*e1+e2*e2)
            errly(i) = datay(i) - e1
            e1 = erruy(i1) - datay(i1)
            e2 = erruy(i2) - datay(i2)
            e1 = sqrt(e1*e1+e2*e2)
            erruy(i) = datay(i) + e1
            modely(i) = modely(i1) + modely(i2)
            i1 = i1 + 1
            i2 = i2 + 1
         ENDDO
         WRITE (ylab(jset), '(2(a,i3))') 'data set', iset(1),
     &                                 '+ data set', iset(2)
      ENDIF
      mbin = nbine(jset) - nbins(jset) + 1
      WRITE (wrtstr, '(2a,i3,a,i4,a)') indent, 'Data set ', jset,
     &                                    ';', mbin, ' bins'
      CALL xwrite(wrtstr, 10)
 900  CONTINUE
      RETURN
 
      ENTRY muldat(string, iparse, dset_new, dset)
c        multiplies two data sets
c        fwj haberl   6-APR-1990 11:18:46
 
      CALL xgtint(string, iparse, 1, 'dataset 1', 1, iset(1), nret,
     &            iflag, idelim)
      IF (iflag.LT.0) THEN
         GOTO 900
      ENDIF
      jset = dset + dset_new
      IF (iset(1).LT.1 .OR. iset(1).GT.jset) THEN
         WRITE (*, '(a,i3)') ' Data set no ', iset(1), ' not available'
         GOTO 900
      ENDIF
      i1 = nbine(iset(1)) - nbins(iset(1)) + 1
      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), subcom, ncom, icom)
      ENDIF
      const = 1.
      cerrl = 1.
      cerru = 1.
      IF (icom.EQ.1) THEN
c        multiplication with given value (including errors)
         CALL xgtrl4(string, iparse, 1, 'const', 1, const, nret, iflag,
     &               idelim)
         CALL xgtrl4(string, iparse, 1, 'error lo', 1, cerrl, nret,
     &               iflag, idelim)
         CALL xgtrl4(string, iparse, 1, 'error hi', 1, cerru, nret,
     &               iflag, idelim)
      ELSEIF (icom.EQ.2) THEN
c        multiplication with x values (no errors on x)
      ELSE
         IF (iend.GE.ibeg) THEN
            CALL xptarg(string, iparse, string(ibeg:iend), iflag,
     &                  idelim)
            CALL xgtint(string, iparse, 1, 'dataset 2', 1, iset(2),
     &                  nret, iflag, idelim)
            IF (iset(2).LT.1 .OR. iset(2).GT.jset) THEN
               WRITE (*, '(a,i3)') ' Data set no ', iset(2),
     &                             ' not available'
               GOTO 900
            ELSE
               i2 = nbine(iset(2)) - nbins(iset(2)) + 1
               IF (i1.NE.i2) THEN
                  WRITE (*, '(a)')
     &                 ' The 2 data sets have different numbers of bins'
                  GOTO 900
               ENDIF
            ENDIF
         ENDIF
      ENDIF
      dset_new = dset_new + 1
      jset = dset + dset_new
      xlab(jset) = xlab(iset(1))
      nbins(jset) = nbine(jset-1) + 1
      nbine(jset) = nbine(jset-1) + i1
      IF (icom.EQ.1) THEN
         i1 = nbins(iset(1))
         DO i = nbins(jset), nbine(jset)
            datax(i) = datax(i1)
            errlx(i) = errlx(i1)
            errux(i) = errux(i1)
            datay(i) = datay(i1)*const
            e1 = (datay(i1)-errly(i1))*const
            e2 = (const-cerrl)*datay(i1)
            e1 = sqrt(e1*e1+e2*e2)
            errly(i) = datay(i)*(1.-e1)
            e1 = (erruy(i1)-datay(i1))*const
            e2 = (cerru-const)*datay(i1)
            e1 = sqrt(e1*e1+e2*e2)
            erruy(i) = datay(i)*(1.+e1)
            modely(i) = modely(i1)*const
            i1 = i1 + 1
         ENDDO
         WRITE (ylab(jset), '(a,i3,a,g15.4)') 
     &         'data set', iset(1), ' * ', const
      ELSEIF (icom.EQ.2) THEN
         i1 = nbins(iset(1))
         DO i = nbins(jset), nbine(jset)
            datax(i) = datax(i1)
            errlx(i) = errlx(i1)
            errux(i) = errux(i1)
            datay(i) = datay(i1)*datax(i1)
            errly(i) = errly(i1)*datax(i1)
            erruy(i) = erruy(i1)*datax(i1)
            modely(i) = modely(i1)*datax(i1)
            i1 = i1 + 1
         ENDDO
         WRITE (ylab(jset), '(a,i3,a)') 
     &         'data set', iset(1), 'x*y'
      ELSE
         i1 = nbins(iset(1))
         i2 = nbins(iset(2))
         DO i = nbins(jset), nbine(jset)
            IF (datax(i1).NE.datax(i2)) THEN
               WRITE (*, '(a)') ' Warning: data has different x-values'
            ENDIF
            datax(i) = datax(i1)
            errlx(i) = errlx(i1)
            errux(i) = errux(i1)
            datay(i) = datay(i1)*datay(i2)
            e1 = (datay(i1)-errly(i1))*datay(i2)
            e2 = (datay(i2)-errly(i2))*datay(i1)
            e1 = sqrt(e1*e1+e2*e2)
            errly(i) = datay(i)*(1.-e1)
            e1 = (erruy(i1)-datay(i1))*datay(i2)
            e2 = (erruy(i2)-datay(i2))*datay(i1)
            e1 = sqrt(e1*e1+e2*e2)
            erruy(i) = datay(i)*(1.+e1)
            modely(i) = modely(i1)*modely(i2)
            i1 = i1 + 1
            i2 = i2 + 1
         ENDDO
         WRITE (ylab(jset), '(2(a,i3))') 'data set', iset(1),
     &                                 '+ data set', iset(2)
      ENDIF
      mbin = nbine(jset) - nbins(jset) + 1
      WRITE (wrtstr, '(2a,i3,a,i4,a)') indent, 'Data set ', jset,
     &                                    ';', mbin, ' bins'
      CALL xwrite(wrtstr, 10)
 901  CONTINUE
      RETURN
 
      ENTRY divdat(string, iparse, dset_new, dset)
c        divides two data sets
c        fwj haberl march 1988
 
c      real*4 loglim
c      data loglim/-1.e30/
 
      CALL xgtint(string, iparse, 2, 'dataset 1', 1, iset, nret, iflag,
     &            idelim)
      IF (iflag.LT.0) THEN
         GOTO 902
      ENDIF
      jset = dset + dset_new
      IF (iset(1).LT.1 .OR. iset(1).GT.jset) THEN
         WRITE (*, '(a,i3,a)') ' Data set no ', iset(1),
     &                         ' not available'
         GOTO 902
      ELSEIF (iset(2).LT.1 .OR. iset(2).GT.jset) THEN
         WRITE (*, '(a,i3,a)') ' Data set no ', iset(2),
     &                         ' not available'
         GOTO 902
      ELSE
         i1 = nbine(iset(1)) - nbins(iset(1)) + 1
         i2 = nbine(iset(2)) - nbins(iset(2)) + 1
         IF (i1.NE.i2) THEN
            WRITE (*, '(a)')
     &                 ' The 2 data sets have different numbers of bins'
            GOTO 902
         ENDIF
         dset_new = dset_new + 1
         jset = dset + dset_new
         xlab(jset) = xlab(iset(1))
         nbins(jset) = nbine(jset-1) + 1
         nbine(jset) = nbine(jset-1) + i1
         i1 = nbins(iset(1))
         i2 = nbins(iset(2))
         DO i = nbins(jset), nbine(jset)
            IF (datax(i1).NE.datax(i2)) THEN
               WRITE (*, '(a)')
     &                     ' Warning: data sets have different x values'
            ENDIF
            datax(i) = datax(i1)
            errlx(i) = errlx(i1)
            errux(i) = errux(i1)
            IF (datay(i2).NE.0. .AND. datay(i1).NE.0.) THEN
               datay(i) = datay(i1)/datay(i2)
               e1 = (datay(i1)-errly(i1))/datay(i1)
               e2 = (datay(i2)-errly(i2))/datay(i2)
               e1 = sqrt(e1*e1+e2*e2)
               errly(i) = datay(i)*(1.-e1)
               e1 = (erruy(i1)-datay(i1))/datay(i1)
               e2 = (erruy(i2)-datay(i2))/datay(i2)
               e1 = sqrt(e1*e1+e2*e2)
               erruy(i) = datay(i)*(1.+e1)
            ELSE
               datay(i) = 0.
               errly(i) = 0.
               erruy(i) = 0.
            ENDIF
            IF (abs(modely(i2)).GT.1.e-10) THEN
               modely(i) = modely(i1)/modely(i2)
            ELSE
               modely(i) = 0.
            ENDIF
c           this is for calculating optical depth
c           if(errly(i).gt.0.)then
c              errly(i)=alog(errly(i))
c           else
c              errly(i)=loglim
c           endif
c           if(erruy(i).gt.0.)then
c              erruy(i)=alog(erruy(i))
c           else
c              erruy(i)=loglim
c           endif
c           if(datay(i).gt.0.)then
c              datay(i)=alog(datay(i))
c           else
c              datay(i)=loglim
c           endif
            i1 = i1 + 1
            i2 = i2 + 1
         ENDDO
         WRITE (ylab(jset), '(2(a,i3))') 'data set', iset(1),
     &                                 '/ data set', iset(2)
c        ylab(jset)='additional optical depth'
      ENDIF
      mbin = nbine(jset) - nbins(jset) + 1
      WRITE (wrtstr, '(2a,i3,a,i4,a)') indent, 'Data set ', jset,
     &                                    ';', mbin, ' bins'
      CALL xwrite(wrtstr, 10)
 902  CONTINUE
      RETURN
      END
