      SUBROUTINE contpl(string, iparse, chi, stepsubx, stepsuby, dset,
     &                  dset_new, qini)
 
c	fwj haberl   15-MAR-1989 14:38:22
c  Plotting of chi-square contours using XSPEC history package 'steppar'
 
      INCLUDE 'xplot.inc'
      INCLUDE '../inc/historysf.inc'
      INTEGER*4 conmax
      PARAMETER (conmax=10)
      REAL*4 chi, chimin(conmax)
      REAL*4 xsq(256, 256)
      REAL*4 trmat(6), clev(20, conmax), baselev(20)
      INTEGER*4 nlev(conmax), dset, dset_new, iset, lenact
      INTEGER*4 nret, iflag, idelim, iparse, nx, ny, i, j, ij, ix, iy
      INTEGER*4 ilev, divxy(2)
      CHARACTER*(*) string
      LOGICAL*4 qini
      CHARACTER*64 lab
      REAL*4 xmin(conmax), xmax(conmax), ymin(conmax), ymax(conmax)
      REAL*4 anxmin, anxmax, anymin, anymax, x, a, b, c, plim(2, 2)
      CHARACTER*12 box_x(xplotmax, yplotmax), box_y(xplotmax, yplotmax)
      RECORD /his_steppar_sub1/stepsubx, stepsuby
      DATA trmat/0., 1., 0., 0., 0., 1./
      DATA nlev/conmax*4/
      DATA baselev/1., 2.71, 4.61, 9.21, 16*1.e4/
      DATA divxy/1, 1/
c
c  STATEMENT FUNCTIONS.
c
      REAL winlim
c	qint(x, a, b)= .8*(x-a)/(b-a) + .1
      winlim(x, a, b, c) = 1. + (c-1.)*(x-a)/(b-a)
 
      IF (dset.EQ.0 .AND. dset_new.EQ.0) THEN
         WRITE (*, '(a)') ' No data to plot read in yet'
         GOTO 1000
      ENDIF
c plot device:
      IF (.NOT.qini) THEN
         WRITE (*, *) ' Plot device not defined yet'
         RETURN
      ENDIF
      CALL pgbegin(0, device, 1, 1)
      IF (dset_new.GT.0) THEN
         iset = dset + dset_new
         IF (iset.GT.conmax) THEN
            WRITE (*, '(a,i2,a)') ' Cannot handle more than ', conmax,
     &                            ' contour plots'
            GOTO 1000
         ENDIF
         chimin(iset) = chi
         CALL xgtrl4(string, iparse, 1, 'chi square', 1, chimin(iset),
     &               nret, iflag, idelim)
         IF (iflag.LT.0) THEN
            GOTO 1000
         ENDIF
         CALL xgtint(string, iparse, 1, 'no of contours', 1, nlev(iset),
     &               nret, iflag, idelim)
         IF (iflag.LT.0) THEN
            GOTO 1000
         ENDIF
         CALL xgtrl4(string, iparse, nlev(iset), 'levels', 1, baselev,
     &               nret, iflag, idelim)
         IF (iflag.LT.0) THEN
            GOTO 1000
         ENDIF
         CALL xgtint(string, iparse, 2, 'plot position', 1, divxy, nret,
     &               iflag, idelim)
         IF (iflag.LT.0) THEN
            GOTO 1000
         ENDIF
         ix = divxy(1)
         iy = divxy(2)
         pos(1, iset) = ix
         pos(2, iset) = iy
         DO ilev = 1, nlev(iset)
            clev(ilev, iset) = chimin(iset) + baselev(ilev)
         ENDDO
         IF (stepsubx.log) qxlog(ix) = .TRUE.
         IF (stepsuby.log) qylog(iy) = .TRUE.
         xmin(iset) = stepsubx.lowpar
         xmax(iset) = stepsubx.hipar
         ymin(iset) = stepsuby.lowpar
         ymax(iset) = stepsuby.hipar
         IF (.NOT.qslimx(ix)) THEN
            IF (limitx(1,ix).EQ.limitx(2,ix)) THEN
               limitx(1, ix) = xmin(iset)
               limitx(2, ix) = xmax(iset)
            ELSE
               IF (xmin(iset).LT.limitx(1,ix)) THEN
                  limitx(1, ix) = xmin(iset)
               ENDIF
               IF (xmax(iset).GT.limitx(2,ix)) THEN
                  limitx(2, ix) = xmax(iset)
               ENDIF
            ENDIF
         ENDIF
         IF (.NOT.qslimx(ix)) THEN
            IF (limity(1,iy).EQ.limity(2,iy)) THEN
               limity(1, iy) = ymin(iset)
               limity(2, iy) = ymax(iset)
            ELSE
               IF (ymin(iset).LT.limity(1,iy)) THEN
                  limity(1, iy) = ymin(iset)
               ENDIF
               IF (ymax(iset).GT.limity(2,iy)) THEN
                  limity(2, iy) = ymax(iset)
               ENDIF
            ENDIF
         ENDIF
         WRITE (lab, '(a,i2)') ' Fit parameter ', stepsubx.fitind
         xlab(ix) = lab
         WRITE (lab, '(a,i2)') ' Fit parameter ', stepsuby.fitind
         ylab(iy) = lab
         dset = iset
         dset_new = 0
      ENDIF
      DO iset = 1, dset
         ix = pos(1, iset)
         iy = pos(2, iset)
         nx = nbinz(iset, 1)
         ny = nbinz(iset, 2)
         ij = 0
         DO i = 1, iset - 1
            ij = ij + nbinz(i, 1)*nbinz(i, 2)
         ENDDO
         DO j = 1, ny
            DO i = 1, nx
               ij = ij + 1
               xsq(i, j) = datax(ij)
            ENDDO
         ENDDO
         IF (qxlog(ix)) THEN
            plim(1, 1) = alog10(limitx(1,ix))
            plim(2, 1) = alog10(limitx(2,ix))
            DO i = 1, divy
               IF (i.EQ.1) THEN
                  box_x(ix, i) = 'BCNSTL'
               ELSE
                  box_x(ix, i) = 'BCSTL'
               ENDIF
            ENDDO
            anxmin = winlim(plim(1,1), alog10(xmin(iset)),
     &               alog10(xmax(iset)), float(nx))
            anxmax = winlim(plim(2,1), alog10(xmin(iset)),
     &               alog10(xmax(iset)), float(nx))
         ELSE
            plim(1, 1) = limitx(1, ix)
            plim(2, 1) = limitx(2, ix)
            DO i = 1, divy
               IF (i.EQ.1) THEN
                  box_x(ix, i) = 'BCNST'
               ELSE
                  box_x(ix, i) = 'BCST'
               ENDIF
            ENDDO
            anxmin = winlim(plim(1,1), xmin(iset), xmax(iset), float(nx)
     &               )
            anxmax = winlim(plim(2,1), xmin(iset), xmax(iset), float(nx)
     &               )
         ENDIF
         IF (qylog(iy)) THEN
            plim(1, 2) = alog10(limity(1,iy))
            plim(2, 2) = alog10(limity(2,iy))
            DO i = 1, divx
               IF (i.EQ.1) THEN
                  box_y(i, iy) = 'BCNSTL'
               ELSE
                  box_y(i, iy) = 'BCSTL'
               ENDIF
            ENDDO
            anymin = winlim(plim(1,2), alog10(ymin(iset)),
     &               alog10(ymax(iset)), float(ny))
            anymax = winlim(plim(2,2), alog10(ymin(iset)),
     &               alog10(ymax(iset)), float(ny))
         ELSE
            plim(1, 2) = limity(1, iy)
            plim(2, 2) = limity(2, iy)
            DO i = 1, divx
               IF (i.EQ.1) THEN
                  box_y(i, iy) = 'BCNST'
               ELSE
                  box_y(i, iy) = 'BCST'
               ENDIF
            ENDDO
            anymin = winlim(plim(1,2), ymin(iset), ymax(iset), float(ny)
     &               )
            anymax = winlim(plim(2,2), ymin(iset), ymax(iset), float(ny)
     &               )
         ENDIF
         CALL pgvport(framex1(ix), framex2(ix), framey1(iy), framey2(iy)
     &                )
         CALL pgwindow(plim(1,1), plim(2,1), plim(1,2), plim(2,2))
         CALL pgsci(colind(4))
         CALL pgbox(box_x(ix,iy), 0.0, 0, box_y(ix,iy), 0.0, 0)
         IF (qabc) CALL pgmtext('T', -1.5, 0.95, 1.0, abc(ix,iy))
         DO i = 1,10
            CALL pgtext (vpos(i,1), vpos(i,2), vlab(i))
         ENDDO
         IF (ix.EQ.1) THEN
            lab = ylab(iy)
            IF (lenact(lab).GT.0) CALL pgmtext('L', 3., 0.5, 0.5, lab(:
     &                                         lenact(lab)))
         ENDIF
         IF (iy.EQ.divy) THEN
            lab = phead(ix)
            IF (lenact(lab).GT.0) CALL pgmtext('T', 1.5, 0.5, 0.5, lab(:
     &                                         lenact(lab)))
            lab = tllab(ix)
            IF (lenact(lab).GT.0) CALL pgmtext('T', 1.5, 0.0, 0.0, lab(:
     &                                         lenact(lab)))
            lab = trlab(ix)
            IF (lenact(lab).GT.0) CALL pgmtext('T', 1.5, 1.0, 1.0, lab(:
     &                                         lenact(lab)))
         ENDIF
         IF (iy.EQ.1) THEN
            lab = xlab(ix)
            IF (lenact(lab).GT.0) CALL pgmtext('B', 3., 0.5, 0.5, lab(:
     &                                         lenact(lab)))
         ENDIF
         CALL pgwindow(anxmin, anxmax, anymin, anymax)
         CALL pgsci(colind(2))
         CALL pgcont(xsq, 256, 256, 1, nx, 1, ny, clev(1,iset),
     &               nlev(iset), trmat)
      ENDDO
      CALL pgend
 1000 CONTINUE
      RETURN
      END
