SUBROUTINE pltfld(pchsen, pchsn2, ixopt, fluxar, ear, nadd, 
     &                  ntchan, iphacmp, y, iery, npts, nvec, cmd, ncmd)

      REAL y(*)
      REAL fluxar(*), ear(0:*)

      INTEGER iery(*), npts, nvec, ncmd, nadd, ntchan
      INTEGER ixopt, iphacmp

      CHARACTER cmd(*)*(*), pchsen*(*), pchsn2*(*)


c---
c XSPEC plotting subroutine to plot the folded model and/or observed data
c---
c arguments
c pchsen    c	 type of plot chosen. One of data, 
c                 counts, log data, residuals, 
c                 contributions to chi-squared,
c                 delta chi, data/model, and integrated counts.
c pchsn2    c    type of 2nd plot chosen. One of
c                 none, residuals, contributions to 
c                 chi-squared, delta chi, data/model
c ixopt     i    x-axis option : 1 - channels
c                                2 - energy
c                                3 - wavelength
c fluxar    r    model fluxes in energy space
c nadd      i    number of additive model components - if there is only one
c                then nadd was set to zero
c ntchan    i    number of bins allocated for each component
c phacmp    r    folded additive model components
c---

      REAL NO
      PARAMETER (NO=-1.2e-34)

      REAL conv, msyst, obs, err, modl, totcou
      REAL xmin, xmax, ymin, ymax
      INTEGER i, j, ioff, ncou, k, status, labno, ntest
      LOGICAL qcount
      CHARACTER xlabel*72, ylabel*72, tlabel*72

      INTEGER dgnplg, cgncmp, dgdpnt
      REAL xgmsys, gammln, xpglkt, xpglem, xpglz
      CHARACTER*11 xgstat
      LOGICAL xpgqid
      EXTERNAL dgnplg, cgncmp, xgstat, xgmsys, gammln, dgdpnt
      EXTERNAL xpgqid, xpglkt, xpglem, xpglz

      IF ( pchsen .NE. 'data' .AND. pchsen .NE. 'counts' .AND.
     &     pchsen .NE. 'ldata' .AND. pchsen .NE. 'residuals' .AND.
     &     pchsen .NE. 'chisq' .AND. pchsen .NE. 'delchi' .AND.
     &     pchsen .NE. 'ratio' .AND. pchsen .EQ. 'foldmodel' .AND.
     &     pchsen .NE. 'icounts' .AND.
     &     pchsn2 .NE. 'none' .AND. pchsn2 .NE. 'residuals' .AND. 
     &     pchsn2 .NE. 'chisq' .AND. pchsn2 .NE. 'delchi' .AND. 
     &     pchsn2 .NE. 'ratio' ) RETURN

c  Set qcount if counts are to be plotted

      qcount = .FALSE.
      IF ( pchsen .EQ. 'counts' .OR. pchsen .EQ. 'icounts' ) 
     &   qcount = .TRUE.

c  If there are folded additive model components to be processed then
c  do so. Load results into y(i+7*npts and higher).

      IF ( ( pchsen .EQ. 'data' .OR. pchsen .EQ. 'counts' .OR.
     &       pchsen .EQ. 'ldata' .OR. pchsen .EQ. 'icounts' ) 
     &      .AND. nadd .GT. 1 ) THEN
         DO j = 1, nadd
            ioff = iphacmp + (j-1)*ntchan
            CALL ldparr(ioff, fluxar, ear, ixopt, y, npts, cmd, ncmd, 
     &                  qcount)
            DO i = 1, npts
               y(i+(j+6)*npts) = y(i+4*npts)
            ENDDO
         ENDDO
      ENDIF

c  Load the plot data into the y array

      CALL ldparr(dgdpnt('mod'), fluxar, ear, ixopt, y, npts, cmd, 
     &            ncmd, qcount)


c  Load the residuals and errors into columns 6 and 7.

      DO i = 1, npts
         IF (y(i) .NE. NO) THEN
            y(i+5*npts) = y(i+2*npts) - y(i+4*npts)
            y(i+6*npts) = y(i+3*npts)
         ELSE
            y(i+5*npts) = NO
            y(i+6*npts) = NO
         ENDIF
      ENDDO

c  The contents of the Y array are now 
c       Y(i)         -  channels/energy
c       Y(i+npts)    -  delta on channels/energy
c       Y(i+2*npts)  -  data
c       Y(i+3*npts)  -  errors on data
c       Y(i+4*npts)  -  folded model
c       Y(i+5*npts)  -  residuals
c       Y(i+6*npts)  -  errors on residuals
c       Y(i+7*npts)  -  folded model component
c    ...Y(i+(nadd+6))

c  Change the y array for those plots which are not simple data,
c  error and folded model.

c  **Case of integral counts

      IF ( pchsen .EQ. 'icounts' ) THEN


         y(1+2*npts) = y(1+2*npts)*2*y(1+npts)
         y(1+4*npts) = y(1+4*npts)*2*y(1+npts)
         ncou = 1
         DO j = 1, nadd
            y(1+(j+6)*npts) = y(1+(j+6)*npts)*2*y(1+npts)
         ENDDO
         DO i = 2, npts
            IF (y(i-1) .EQ. NO) THEN
               y(i+2*npts) = y(i+2*npts)*2*y(i+npts)
               y(i+4*npts) = y(i+4*npts)*2*y(i+npts)
               DO j = 1, nadd
                  y(i+(j+6)*npts) = y(i+(j+6)*npts)*2*y(i+npts)
               ENDDO
               ncou = 1
            ELSEIF (y(i) .NE. NO) THEN
               ncou = ncou + 1
               y(i+2*npts) = y(i+2*npts)*2*y(i+npts) + y(i-1+2*npts)
               y(i+4*npts) = y(i+4*npts)*2*y(i+npts) + y(i-1+4*npts)
               DO j = 1, nadd
                  y(i+(j+6)*npts) = y(i+(j+6)*npts)*2*y(i+npts)
     &                              + y(i-1+(j+6)*npts)
               ENDDO
            ELSEIF (y(i) .EQ. NO) THEN
               totcou = y(i-1+2*npts)
               DO k = i-ncou, i-1
                  y(k+2*npts) = y(k+2*npts)/totcou
                  y(k+4*npts) = y(k+4*npts)/totcou
                  DO j = 1, nadd
                     y(k+(j+6)*npts) = y(k+(j+6)*npts)/totcou
                  ENDDO
               ENDDO
            ENDIF
         ENDDO
         totcou = y(3*npts)
         IF ( totcou .NE. 0. ) THEN
            DO k = npts-ncou+1, npts
               y(k+2*npts) = y(k+2*npts)/totcou
               y(k+4*npts) = y(k+4*npts)/totcou
               DO j = 1, nadd
                  y(k+(j+6)*npts) = y(k+(j+6)*npts)/totcou
               ENDDO
            ENDDO
         ENDIF
         DO i = 1, npts
            IF ( y(i) .NE. NO ) THEN
               y(i+5*npts) = y(i+2*npts) - y(i+4*npts)
            ENDIF
         ENDDO

c  **Case of just want the folded model

      ELSEIF ( pchsen .EQ. 'foldmodel' ) THEN

         DO i = 1, npts
            IF (y(i) .NE. NO) THEN
               y(i+2*npts) = y(i+4*npts)
            ENDIF
         ENDDO
         DO j = 1, nadd
            DO i = 1, npts
               y(i+(j+2)*npts) = y(i+(j+6)*npts)
            ENDDO
         ENDDO

c  **Case of residuals

      ELSEIF ( pchsen .EQ. 'residuals' ) THEN

         DO i = 1, npts
            IF (y(i) .NE. NO) THEN
               y(i+2*npts) = y(i+5*npts)
               y(i+3*npts) = y(i+6*npts)
            ENDIF
         ENDDO

c  **Case of contributions to statistic

      ELSEIF ( pchsen .EQ. 'chisq' ) THEN

         msyst = xgmsys()**2

         IF ( xgstat() .EQ. 'Chi-Squared') THEN
            DO i = 1, npts
               IF (y(i) .NE. NO) THEN
                  obs = y(i+2*npts)
                  err = y(i+3*npts)
                  modl = y(i+4*npts)
                  y(i+2*npts) = SIGN( (obs-modl)**2, obs-modl)
                  IF ( msyst .GT. 0. ) THEN
                     y(i+2*npts) = y(i+2*npts)/
     &                 (err**2+msyst*modl**2)
                  ELSE
                     y(i+2*npts) = y(i+2*npts)/err**2
                  ENDIF
               ENDIF
            ENDDO

c  for the max-likelihood case what we actually plot is the delta likelihood
c  between the current observation and an observation that matched the model
c  precisely.

         ELSEIF ( xgstat() .EQ. 'C-statistic') THEN
            DO i = 1, npts
               IF (y(i) .NE. NO) THEN
                  obs = y(i+2*npts)
                  err = y(i+3*npts)
                  modl = y(i+4*npts)
                  IF ( obs .GT. 0) THEN
                     conv = obs/err**2
                     obs = obs*conv
                     modl = modl*conv
                     IF (modl .GT. 0) THEN
                        y(i+2*npts) = SIGN(
     &                     2*(modl-obs*(1-(log(obs)-log(modl)))), 
     &                        obs-modl)
                     ELSE
                        y(i+2*npts) = NO
                     ENDIF
                  ELSE
                     y(i+2*npts) = 2*modl/err
                  ENDIF
               ENDIF
            ENDDO
         ENDIF

c  **Case of delta chi

      ELSEIF ( pchsen .EQ. 'delchi' ) THEN

         IF ( xgstat() .EQ. 'C-statistic') RETURN

         msyst = xgmsys()**2

         DO i = 1, npts
            IF (y(i) .NE. NO) THEN
               obs = y(i+2*npts)
               err = y(i+3*npts)
               modl = y(i+4*npts)
               IF ( msyst .GT. 0. ) THEN
                  y(i+2*npts) = (obs-modl)/SQRT(err**2 + msyst*modl**2)
               ELSE
                  y(i+2*npts) = (obs-modl)/err
               ENDIF
               y(i+3*npts) = 1.
            ENDIF
         ENDDO

c  **Case of ratio

      ELSEIF ( pchsen .EQ. 'ratio' ) THEN

         DO i = 1, npts
            IF (y(i) .NE. NO .AND. y(i+4*npts) .GT. 1.0E-34) THEN
               IF ( y(i+2*npts) .NE. 0. ) 
     &                 y(i+2*npts) = y(i+2*npts)/y(i+4*npts)
               IF ( y(i+3*npts) .NE. 0. ) 
     &                 y(i+3*npts) = y(i+3*npts)/y(i+4*npts)
            ELSE
               y(i+2*npts) = NO
               y(i+3*npts) = NO
            ENDIF
         ENDDO

      ENDIF

c  If no lower pane is required then shift down all the folded model components

      IF ( pchsn2 .EQ. 'none' .AND. nadd .GT. 1 ) THEN
         DO j = 1, nadd
            DO i = 1, npts
               y(i+(j+4)*npts) = y(i+(j+6)*npts)
            ENDDO
         ENDDO
      ENDIF

c  Now reset y(i+5*npts) and y(i+6*npts) for cases when the second panel is
c  not simple residuals

c  **Case of contributions to chi^2 as a lower panel

      IF ( pchsn2 .EQ. 'chisq' ) THEN

         IF ( xgstat() .EQ. 'Chi-Squared') THEN
            msyst = xgmsys()**2
            DO i = 1, npts
               IF (y(i) .NE. NO) THEN
                  obs = y(i+2*npts)
                  err = y(i+3*npts)
                  modl = y(i+4*npts)
                  y(i+5*npts) = SIGN( (obs-modl)**2, obs-modl)
                  IF ( msyst .LE. 0 ) THEN
                     y(i+5*npts) = y(i+5*npts)/err**2
                  ELSE
                     y(i+5*npts) = y(i+5*npts)/(err**2+msyst*modl**2)
                  ENDIF
                  y(i+6*npts) = 0.
               ENDIF
            ENDDO

c  for the max-likelihood case what we actually plot is the delta likelihood
c  between the current observation and an observation that matched the model
c  precisely.

         ELSEIF ( xgstat() .EQ. 'C-statistic') THEN
            DO i = 1, npts
               IF (y(i) .NE. NO) THEN
                  obs = y(i+2*npts)
                  err = y(i+3*npts)
                  modl = y(i+4*npts)
                  IF ( obs .GT. 0) THEN
                     conv = obs/err**2
                     obs = obs*conv
                     modl = modl*conv
                     IF (modl .GT. 0) THEN
                        y(i+5*npts) = SIGN(
     &                     2*(modl-obs*(1-(log(obs)-log(modl)))), 
     &                        obs-modl)
                     ELSE
                        y(i+5*npts) = NO
                     ENDIF
                  ELSE
                     y(i+5*npts) = 2*modl/err
                  ENDIF
               ENDIF
            ENDDO
         ENDIF

c  **Case of delta chi in the lower pane

      ELSEIF ( pchsn2 .EQ. 'delchi' ) THEN

         IF ( xgstat() .EQ. 'C-statistic') RETURN

         msyst = xgmsys()**2
         DO i = 1, npts
            IF (y(i) .NE. NO) THEN
               obs = y(i+2*npts)
               err = y(i+3*npts)
               modl = y(i+4*npts)
               IF ( msyst .GT. 0. ) THEN
                  y(i+5*npts) = (obs-modl)/SQRT(err**2 + msyst*modl**2)
               ELSE
                  y(i+5*npts) = (obs-modl)/err
               ENDIF
               y(i+6*npts) = 1.
            ENDIF
         ENDDO

c  **Case of ratio in the lower pane

      ELSEIF ( pchsn2 .EQ. 'ratio' ) THEN

         DO i = 1, npts
            IF (y(i) .NE. NO .AND. y(i+4*npts) .GT. 1.0E-34) THEN
               IF ( y(i+2*npts) .NE. 0. ) THEN
                  y(i+5*npts) = y(i+2*npts)/y(i+4*npts)
               ELSE
                  y(i+5*npts) = 0.
               ENDIF
               IF ( y(i+3*npts) .NE. 0. ) THEN
                  y(i+6*npts) = y(i+3*npts)/y(i+4*npts)
               ELSE
                  y(i+6*npts) = 0.
               ENDIF
            ELSE
               y(i+5*npts) = NO
               y(i+6*npts) = NO
            ENDIF
         ENDDO

      ENDIF

c  Keep track of the numbered labels

      labno = 0

c  For data or counts plots make sure that the y-axis starts at zero

      IF ( pchsen .EQ. 'data' .OR. pchsen .EQ. 'counts' .OR.
     &     pchsen .EQ. 'icounts' ) THEN
         cmd(ncmd+1) = 'r y 0'
         ncmd = ncmd + 1
      ENDIF

c  Draw the horizontal line through zero or one for the different
c  types of residual plots

      IF ( pchsen .EQ. 'ratio' ) THEN
         WRITE (cmd(ncmd+1), 510) labno + 1
         ncmd = ncmd + 1
         labno = labno + 1
 510     FORMAT ('la', i3, ' pos .001,1 li 0,100. co 3 " "')
      ELSEIF ( pchsen .EQ. 'residuals' .OR. pchsen .EQ. 'chisq' .OR. 
     &         pchsen .EQ. 'delchi' ) THEN
         WRITE (cmd(ncmd+1), 511) labno + 1
         ncmd = ncmd + 1
         labno = labno + 1
 511     FORMAT ('la', i3, ' pos .001,0 li 0,100. co 3 " "')
      ENDIF

c  Set up the top label

      IF ( pchsen .EQ. 'delchi' ) THEN
         tlabel = 'label t delta '//char(92)//'gx'
      ELSEIF ( pchsen .EQ. 'ratio' ) THEN
         tlabel = 'label t data/model'
      ELSEIF ( pchsen .EQ. 'chisq' ) then
         IF ( xgstat() .EQ. 'Chi-Squared') THEN
            tlabel = 'label t contributions to '//char(92)
     &                    //'gx'//char(92)//'u2'
         ELSE
            tlabel = 'label t contributions to C-statistic'
         ENDIF
      ELSEIF ( pchsen .EQ. 'residuals' ) then
         tlabel = 'label t residuals'
      ELSEIF ( pchsen .EQ. 'data' .OR. pchsen .EQ. 'counts' .OR.
     &         pchsen .EQ. 'ldata' ) THEN
         tlabel = 'label t data and folded model'
      ELSEIF ( pchsen .EQ. 'icounts' ) THEN
         tlabel = 'label t integral counts and folded model'
      ELSEIF ( pchsen .EQ. 'foldmodel' ) THEN
         tlabel = 'label t folded model'
      ENDIF

c Set up the y axis label

      IF ( pchsen .EQ. 'delchi' ) THEN
         ylabel = 'la y '//char(92)//'gx'
      ELSEIF ( pchsen .EQ. 'ratio' ) THEN
         ylabel = 'la y ratio'
      ELSEIF ( pchsen .EQ. 'chisq' ) then
         IF ( xgstat() .EQ. 'Chi-Squared') THEN
            ylabel = 'la y sign(data-model) * '//char(92)//'gx'
     &                    //char(92)//'u2'
         ELSE
            ylabel = 'la y sign(data-model) * '//char(92)//'gD'
     &                    //'C-statistic'
         ENDIF
      ELSEIF ( pchsen .EQ. 'counts' ) THEN
         IF ( ixopt .EQ. 1 ) THEN
            ylabel = 'la y counts/chan'
         ELSEIF ( ixopt .EQ. 2 ) THEN
            ylabel = 'la y counts/keV'
         ELSEIF ( ixopt .EQ. 3 ) THEN
            ylabel = 'la y counts/Angstrom'
         ENDIF
      ELSEIF ( pchsen .EQ. 'icounts' ) THEN
         ylabel = 'la y integral counts'
      ELSE
         IF ( ixopt .EQ. 1 ) THEN
            ylabel = 'la y normalized counts/sec/chan'
         ELSEIF ( ixopt .EQ. 2 ) THEN
            ylabel = 'la y normalized counts/sec/keV'
         ELSEIF ( ixopt .EQ. 3 ) THEN
            ylabel = 'la y normalized counts/sec/Angstrom'
         ENDIF
      ENDIF

c Set up the x axis label

      IF ( ixopt .EQ. 1 ) THEN
         xlabel = 'la x channels'
      ELSEIF ( ixopt .EQ. 2 ) THEN
         xlabel = 'la x channel energy (keV)'
      ELSEIF ( ixopt .EQ. 3 ) THEN
         xlabel = 'la x channel wavelength (Angstrom)'
      ENDIF

c If x-axis is energy or wavelength then plot as log

      IF ( ixopt .EQ. 2 .OR. ixopt .EQ. 3 ) THEN
         cmd(ncmd+1) = 'log x'
         ncmd = ncmd + 1
      ENDIF

c If a second pane has been requested then set up the vertical plot.

      IF ( ( pchsen .EQ. 'data' .OR. pchsen .EQ. 'counts' .OR. 
     &       pchsen .EQ. 'ldata' .OR. pchsen .EQ. 'icounts' ) 
     &     .AND. pchsn2 .NE. 'none' ) THEN

c The top window so include the data and the folded model.

         cmd(ncmd+1) = 'window 1'
         cmd(ncmd+2) = 'location 0,0.302,1,0.9'
         cmd(ncmd+3) = 'lab nx off'
         IF ( pchsen .NE. 'icounts' ) cmd(ncmd+4) = 'gap error'
         ncmd = ncmd + 4
         DO i = 1, dgnplg()
            IF ( nadd .GT. 1 ) THEN
               WRITE(cmd(ncmd+1),
     &            '(''yplot on'',i4,'','',i4,'','',i4,''..'',i4)') 
     &               1+(i-1)*(nadd+3), 2+(i-1)*(nadd+3), 
     &               4+(i-1)*(nadd+3), nadd+3+(i-1)*(nadd+3)
            ELSE
               WRITE(cmd(ncmd+1),
     &            '(''yplot on'',i4,'','',i4)') 1+(i-1)*3, 2+(i-1)*3
            ENDIF
            ncmd = ncmd + 1
            WRITE(cmd(ncmd+1),'(''line step on '',i4)') 
     &               2+(i-1)*(nadd+3)
            WRITE(cmd(ncmd+2),'(''color '', i4, '' on '',i4)') 
     &                        MOD(i-1,15)+1, 1+(i-1)*(nadd+3)
            WRITE(cmd(ncmd+3),'(''color '', i4, '' on '',i4)') 
     &                        MOD(i-1,15)+1, 2+(i-1)*(nadd+3)
            ncmd = ncmd + 3
            DO j = 1, nadd
               WRITE(cmd(ncmd+1),'(''line step on '',i4)') 
     &               j+3+(i-1)*(nadd+3)
               WRITE(cmd(ncmd+2),'(''color '', i4, '' on '',i4)') 
     &               MOD(i-1,15)+1, j+3+(i-1)*(nadd+3)
               WRITE(cmd(ncmd+3),'(''lstyle '', i4, '' on '',i4)') 
     &               j+1, j+3+(i-1)*(nadd+3)
               ncmd = ncmd + 3
            ENDDO
         ENDDO
         cmd(ncmd+1) = ylabel
         cmd(ncmd+2) = tlabel
         ncmd = ncmd + 2

c If the ldata option was selected then set log on the y-axis

         IF ( pchsen .EQ. 'ldata' ) THEN
            cmd(ncmd+1) = 'log y on 1'
            ncmd = ncmd + 1
         ENDIF

c The bottom window for the various types of residuals

         cmd(ncmd+1) = 'window 2'
         cmd(ncmd+2) = 'location 0,0.1,1,0.39'
         IF ( pchsen .NE. 'icounts' ) cmd(ncmd+3) = 'gap error'
         cmd(ncmd+4) = 'mmaster 1 x 2'
         ncmd = ncmd + 4
         DO i = 1, dgnplg()
            WRITE(cmd(ncmd+1),'(''yplot on'',i4)') 3+(i-1)*(nadd+3)
            ncmd = ncmd + 1
            IF ( pchsn2 .EQ. 'chisq' ) THEN
               WRITE(cmd(ncmd+1),'(''line step on '',i4)') 
     &                      3+(i-1)*(nadd+3)
               ncmd = ncmd + 1
            ENDIF
            WRITE(cmd(ncmd+1),'(''color '', i4, '' on '',i4)') 
     &                        MOD(i-1,15)+1, 3+(i-1)*(nadd+3)
            ncmd = ncmd + 1
         ENDDO
         IF ( pchsn2 .EQ. 'residuals' ) THEN
            cmd(ncmd+1) = 'la y residuals'
         ELSEIF ( pchsn2 .EQ. 'chisq' ) THEN
            IF ( xgstat() .EQ. 'Chi-Squared') THEN
               cmd(ncmd+1) = 'la y sign(d-m)*'//char(92)//'gx'//
     &                       char(92)//'u2'
            ELSE
               cmd(ncmd+1) = 'la y sign(d-m)*C'
            ENDIF
         ELSEIF ( pchsn2 .EQ. 'delchi' ) THEN
            cmd(ncmd+1) = 'la y '//char(92)//'gx'
         ELSEIF ( pchsn2 .EQ. 'ratio' ) THEN
            cmd(ncmd+1) = 'la y ratio'
         ENDIF
         ncmd = ncmd + 1
         IF ( pchsn2 .EQ. 'ratio' ) THEN
            WRITE(cmd(ncmd+1),'(a,i3,a)') 'la ', labno+1,
     &        '  pos .001,1 li 0,100. co 1 lst 2 " "'
         ELSE
            WRITE(cmd(ncmd+1),'(a,i3,a)') 'la ', labno+1,
     &        '  pos .001,0 li 0,100. co 1 lst 2 " "'
         ENDIF
         cmd(ncmd+2) = xlabel
         cmd(ncmd+3) = 'window 1'
         cmd(ncmd+4) = 'lab nx on 2'
         ncmd = ncmd + 4
         labno = labno + 1

c The standard single pane plot so write the labels and set the log y axis
c if requested

      ELSE

         cmd(ncmd+1) = ylabel
         cmd(ncmd+2) = tlabel
         cmd(ncmd+3) = xlabel
         ncmd = ncmd + 3
         IF ( pchsen .EQ. 'ldata' ) THEN
            cmd(ncmd+1) = 'log y'
            ncmd = ncmd + 1
         ENDIF

c If folded model then turn on the stepped line and colors.

         IF ( ( pchsen .EQ. 'data' .OR. pchsen .EQ. 'counts' .OR.
     &          pchsen .EQ. 'ldata' .OR. pchsen .EQ. 'foldmodel' .OR.
     &          pchsen .EQ. 'icounts' ) 
     &          .AND. cgncmp() .GT. 0 ) THEN

            DO i = 1, dgnplg()
               WRITE(cmd(ncmd+1),'(''line step on '',i4)') 
     &                        2+(i-1)*(nadd+2)
               WRITE(cmd(ncmd+2),'(''color '', i4, '' on '',i4)')
     &                        MOD(i-1,15)+1, 1+(i-1)*(nadd+2)
               WRITE(cmd(ncmd+3),'(''color '', i4, '' on '',i4)')
     &                        MOD(i-1,15)+1, 2+(i-1)*(nadd+2)
               ncmd = ncmd + 3
               DO j = 1, nadd
                  WRITE(cmd(ncmd+1),'(''line step on '',i4)') 
     &                  j+2+(i-1)*(nadd+2)
                  WRITE(cmd(ncmd+2),'(''color '', i4, '' on '',i4)') 
     &                  MOD(i-1,15)+1, j+2+(i-1)*(nadd+2)
                  WRITE(cmd(ncmd+3),'(''lstyle '', i4, '' on '',i4)') 
     &                  j+1, j+2+(i-1)*(nadd+2)
                  ncmd = ncmd + 3
               ENDDO
            ENDDO

c If contributions to chi-squared then turn on the stepped line

         ELSEIF ( pchsen .EQ. 'chisq' ) THEN

            DO i = 1, dgnplg()
               WRITE(cmd(ncmd+1),'(''line step on '',i4)') i
               ncmd = ncmd + 1
            ENDDO

         ENDIF

      ENDIF

c If integral counts then turn the errors off and set the line step on

      IF ( pchsen .EQ. 'icounts' ) THEN
         cmd(ncmd+1) = 'error off'
         ncmd = ncmd + 1
         IF ( cgncmp() .GT. 0 ) THEN
            DO i = 1, dgnplg()
               WRITE(cmd(ncmd+1),'(''line step on '',i4)') 
     &                   1+(i-1)*(nadd+3)
               ncmd = ncmd + 1
            ENDDO
         ELSE
            DO i = 1, dgnplg()
               WRITE(cmd(ncmd+1),'(''line step on '',i4)') i
               ncmd = ncmd + 1
            ENDDO
         ENDIF
      ENDIF

c Temporary fix for PLT bug

      IF ( pchsen .NE. 'chisq' .AND. pchsen .NE. 'icounts' ) THEN
         cmd(ncmd+1) = 'line off 1'
         ncmd = ncmd + 1
      ENDIF

c Set the number of vectors and the array that specifies whether an error
c is given. Notice that each vector includes all the plot groups so the number
c of vectors are the number of different quantities plotted. If no model is
c defined (as specified by cgncmp()=0) then only data and errors are plotted.

      iery(1) = 1
      iery(2) = 1
      iery(3) = 0
      iery(4) = 1
      IF ( ( pchsen .EQ. 'data' .OR. pchsen .EQ. 'counts' .OR. 
     &       pchsen .EQ. 'ldata' .OR. pchsen .EQ. 'icounts' ) 
     &     .AND. nadd .GT. 1 ) THEN
         IF ( pchsn2 .EQ. 'none' ) THEN
            DO j = 1, nadd
               iery(j+3) = 0
            ENDDO
         ELSE
            DO j = 1, nadd
               iery(j+4) = 0
            ENDDO
         ENDIF
      ENDIF
      IF ( pchsen .EQ. 'foldmodel' ) THEN
         iery(2) = 0
         DO j = 1, nadd
            iery(j+2) = 0
         ENDDO
      ENDIF

      IF ( pchsen .EQ. 'chisq' ) iery(2) = 0

      IF ( ( pchsen .EQ. 'data' .OR. pchsen .EQ. 'counts' .OR. 
     &       pchsen .EQ. 'ldata' .OR. pchsen .EQ. 'icounts' ) 
     &      .AND. cgncmp() .GT. 0 ) THEN
         nvec = 3
         IF ( nadd .GT. 1 ) nvec = nvec + nadd
      ELSE
         nvec = 2
      ENDIF

      IF ( pchsen .EQ. 'foldmodel' .AND. nadd .GT. 1 ) THEN
         nvec = nvec + nadd
      ENDIF
      ntest = nvec

      IF ( ( pchsen .EQ. 'data' .OR. pchsen .EQ. 'counts' .OR. 
     &       pchsen .EQ. 'ldata' .OR. pchsen .EQ. 'icounts' ) 
     &     .AND. pchsn2 .NE. 'none' ) THEN
         nvec = nvec + 1
      ENDIF

c Find min and max for the X and Y axes

      xmin = 1.e32
      xmax = -1.e32
      DO i = 1, npts
         IF ( y(i) .NE. NO ) THEN
            xmin = MIN(xmin, y(i)-y(i+npts))
            xmax = MAX(xmax, y(i)+y(i+npts))
         ENDIF
      ENDDO

      ymin = 1.e32
      ymax = -1.e32
      ioff = 2*npts
      DO j = 2, ntest
         DO i = 1, npts
            IF ( y(i+ioff) .NE. NO ) THEN
               ymin = MIN(ymin, y(i+ioff))
               ymax = MAX(ymax, y(i+ioff))
            ENDIF
         ENDDO
         ioff = ioff + (iery(j)+1)*npts
      ENDDO

c  For data or counts plots the y-axis minimum is forced to be zero

      IF ( pchsen .EQ. 'data' .OR. pchsen .EQ. 'counts' .OR.
     &     pchsen .EQ. 'icounts' ) THEN
         ymin = 0.
      ENDIF

c Rescale the x-axis (this should override the PLT/PGPLOT default for log
c axes which can leave blank space at each end)

      WRITE(cmd(ncmd+1), *) 'R X ', xmin, xmax
      ncmd = ncmd + 1

c If required make the labels to plot line IDs. Note that we can only do
c this if we have energy or wavelength on the x-axis.

      IF ( xpgqid() .AND. ixopt .NE. 1 ) THEN
         CALL mklnid(xmin, xmax, ymin, ymax, xpglkt(), xpglem(), 
     &               xpglz(), 1000, cmd, ncmd, labno, status)
      ENDIF

      RETURN
      END