 
      SUBROUTINE pltmod(ear, fluxar, y, iery, npts, nvec, cmd, 
     &                  ncmd, esavar, ixopt, nkeV)
 
      INCLUDE '../../inc/xspec.inc'
 
      INTEGER iery(*), npts, nvec, ncmd, ixopt, nkeV
      REAL ear(0:*), fluxar(*)
      REAL esavar(*)
      REAL y(*)
      CHARACTER cmd(*)*(*)
 
C---
C XSPEC subroutine to plot the current model
C---
C EAR       I    Array of energy bin boundarys
C FLUXAR    I    Total number of counts in energy range
C Y         I/O  Array used by PLT.
C IERY      I/O  Used by PLT.
C CMD       I/O  Command array used by PLT.
C ESAVAR    I	 Array of calculated model components
C ixopt     I    If 1 or 2 then plot vs energy, if 3 vs wavelength
C nkeV      I    Number of keV or A to multiply by
C---
C 9 Jan 1983  - rashafer
C 1.0:	Plot only the total model (as contained in the fluxar array)
C 2.0:	Also plots additive components
 
C---
      REAL NO, KEVTOA
      PARAMETER(NO = -1.2E-34, KEVTOA = 12.398419)
 
      REAL param(MAXMPR)
      INTEGER ie, icomp, ikey, ncomp, jkey, idset, iy, impar
      INTEGER ielt, nmdtg, ncstrt, modbgn, ies, ieear
      INTEGER iye, istat, labno, ioff
      REAL xmin, xmax, ymin, ymax, cnorm
 
      REAL cgnorm, xpglkt, xpglem, xpglz
      DOUBLE PRECISION mgpval
      INTEGER cgmgrp, dgndst, rgnenr, rgdtse, rgcmpe
      INTEGER rgbear, ggbcmp, ggecmp, dgdtgr, cgbmpr, cgempr
      LOGICAL cgisad, cgismu, cgisco, xpgqid
      EXTERNAL cgnorm, cgisad, cgismu, cgmgrp, dgndst, rgnenr, cgisco
      EXTERNAL rgdtse, rgcmpe, rgbear, ggbcmp, ggecmp, dgdtgr
      EXTERNAL cgbmpr, cgempr, mgpval, xpgqid, xpglkt, xpglem, xpglz
 
c Turn skip on to separate the datagroups
 
      cmd(ncmd + 1) = 'SKIP ON'
      ncmd = ncmd + 1
 
      ies = 0
 
c  Calculate the total number of X-axis points in the plot
 
      npts = 0
      DO idset = 1, dgndst() - 1
         npts = npts + rgnenr(idset) + 1
      ENDDO
      npts = npts + rgnenr(dgndst())
 
C Loop over datasets

      xmin = 1.e30
      xmax = -1.e30
 
      DO idset = 1, dgndst()
 
         ncomp = 3
 
c Set the X-axis array
 
         IF( ixopt.EQ.1 .OR. ixopt.EQ.2 )THEN
            DO ie = 1, rgnenr(idset)
               iy = ie + ies
               ieear = ie + rgbear(idset)
               y(iy) = (ear(ieear) + ear(ieear-1))/2.
               y(iy + npts) = (ear(ieear) - ear(ieear-1))/2.
               xmin = MIN(xmin, y(iy)-y(iy+npts))
               xmax = MAX(xmax, y(iy)+y(iy+npts))
            ENDDO
         ELSEIF( ixopt.EQ.3 )THEN
            DO ie = 1, rgnenr(idset)
               iy = ie + ies
               ieear = ie + rgbear(idset)
               IF ( ear(ieear-1) .LE. 0. .OR. ear(ieear) .LE. 0. ) THEN
                  y(iy) = NO
                  y(iy + npts) = NO
               ELSE
                  y(iy) = ((KEVTOA/ear(ieear))+(KEVTOA/ear(ieear-1)))/2.
                  y(iy + npts) = ((KEVTOA/ear(ieear-1)) 
     &                              - (KEVTOA/ear(ieear)))/2.
               ENDIF
               xmin = MIN(xmin, y(iy)-y(iy+npts))
               xmax = MAX(xmax, y(iy)+y(iy+npts))
            ENDDO
         ENDIF
 
C ... the total flux.
 
         ielt = rgdtse(idset)
         DO ie = 1, rgnenr(idset)
            iy = ie + ies
            ieear = ie + rgbear(idset)
            IF( y(iy+npts).EQ.0. )THEN
               IF( iy.GT.ies + 1 )THEN
                  y(iy + (ncomp-1)*npts) = y(iy - 1 + (ncomp-1)*npts)
               ELSE
                  y(iy + (ncomp-1)*npts) = 0.
               ENDIF
            ELSE
               y(iy + (ncomp-1)*npts) = fluxar(ie + ielt)
            ENDIF
         ENDDO
 
         ncstrt = ncomp + 1
 
C Set up the components.
 
         nmdtg = ggecmp(dgdtgr(idset))
         modbgn = ggbcmp(dgdtgr(idset))
 
c Loop round the components for this dataset
 
         DO ikey = modbgn, nmdtg
 
c If it's additive then make another vector
 
            IF( cgisad(ikey) )THEN
 
               ncomp = ncomp + 1
               ielt = rgcmpe(idset, ikey)
               cnorm = cgnorm(ikey)
               DO ie = 1, rgnenr(idset)
                  iy = ie + ies
                  ieear = ie + rgbear(idset)
                  IF( y(iy+npts).EQ.0. )THEN
                     y(iy + (ncomp-1)*npts) = 0.
                  ELSE
                     y(iy + (ncomp-1)*npts) = esavar(ie + ielt)*cnorm
                  ENDIF
               ENDDO
 
c Now apply all multiplicative and convolution components associated
c with this additive component, either as components of this group
c or overall components.
 
               IF( cgmgrp(ikey).NE. - 1 )THEN
                  DO jkey = nmdtg, modbgn, -1
                     IF( cgmgrp(ikey).EQ.cgmgrp(jkey) )THEN
                        ielt = rgcmpe(idset, jkey)
                        IF( cgismu(jkey) )THEN
 
                           DO ie = 1, rgnenr(idset)
                              iy = ie + ies
                              y(iy + (ncomp-1)*npts)
     &                           = y(iy + (ncomp-1)*npts)
     &                           *esavar(ie + ielt)
                           ENDDO
 
                        ELSEIF( cgisco(jkey) )THEN
 
                           DO impar = cgbmpr(jkey), cgempr(jkey)
                              param(impar - cgbmpr(jkey) + 1)
     &                           = sngl(mgpval(impar,'v'))
                           ENDDO
                           CALL udmget( npts, 6, iye, istat )
			   IF( istat .NE. 0 ) CALL xwrite(
     &        ' WARNING: Memory allocation failed in PLTMOD.', 2)

                           CALL conmod(jkey, ear(rgbear(idset)), 
     &                             rgnenr(idset), param, idset, 
     &                             y(ies+1+(ncomp-1)*npts), memr(iye))
 
                           CALL udmfre( iye, 6, istat )     ! error not used
                        ENDIF
                     ENDIF
                  ENDDO
 
                  DO jkey = nmdtg, modbgn, -1
                     IF( cgmgrp(jkey).EQ. - 1 )THEN
                        ielt = rgcmpe(idset, jkey)
                        IF( cgismu(jkey) )THEN
 
                           DO ie = 1, rgnenr(idset)
                              iy = ie + ies
                              y(iy + (ncomp-1)*npts)
     &                           = y(iy + (ncomp-1)*npts)
     &                           *esavar(ie + ielt)
                           ENDDO
 
                        ELSEIF( cgisco(jkey) )THEN
 
                           DO impar = cgbmpr(jkey), cgempr(jkey)
                              param(impar - cgbmpr(jkey) + 1)
     &                           = sngl(mgpval(impar,'v'))
                           ENDDO
                           CALL udmget( npts, 6, iye, istat )
			   IF( istat .NE. 0 ) CALL xwrite(
     &        ' WARNING: Memory allocation failed in PLTMOD.', 2)

                           CALL conmod(jkey, ear(rgbear(idset)), 
     &                          rgnenr(idset), param, idset, 
     &                          y(ies+1+(ncomp-1)*npts),memr(iye))

                           CALL udmfre( iye, 6, istat )   ! for now the error is not plotted
 
                        ENDIF
                     ENDIF
                  ENDDO
               ENDIF
 
               DO ie = 2, rgnenr(idset)
                  iy = ie + ies
                  IF( y(iy).EQ.0. )y(iy + (ncomp-1)*npts)
     &                = y(iy - 1 + (ncomp-1)*npts)
               ENDDO
 
            ENDIF
 
         ENDDO
 
c  If there was only one additive component then it is the same as
c  the total model so don't make a separate plot vector
 
         IF( ncomp.EQ.ncstrt )ncomp = ncomp - 1
 
         ies = ies + rgnenr(idset)
 
c  Write a row of NOs to indicate the end of a dataset
 
         IF( idset.NE.dgndst() )THEN
            ies = ies + 1
            DO icomp = 1, ncomp
               y(ies + (icomp-1)*npts) = NO
            ENDDO
         ENDIF
 
c  End loop over datasets
 
      ENDDO
 
c  If the plot required is in keV or keV^2 (or A or A^2) rather than photons 
c  then modify all y-axis components. Using the geometric rather than 
c  arithmetic means is a slightly better approximation when the bin sizes 
c  are large.
 
      IF ( nkeV .GT. 0 )THEN
         DO icomp = 3, ncomp
            ioff = (icomp-1)*npts
            DO ie = 1, npts
               IF ( y(ie) .NE. NO ) THEN
                  y(ie+ioff) = y(ie+ioff)
     &               *((y(ie)-y(ie+npts))*(y(ie)+y(ie+npts)))**(nkeV/2.)
               ENDIF
            ENDDO
         ENDDO
      ENDIF
 
c Convert from flux to flux/keV or flux/A as appropriate
 
      DO icomp = 3, ncomp
         DO ie = 1, npts
            IF( y(ie+npts).EQ.0. )THEN
               y(ie + (icomp-1)*npts) = NO
            ELSE
               IF( y(ie+(icomp-1)*npts).NE.NO )y(ie + (icomp-1)*npts)
     &             = y(ie + (icomp-1)*npts)/(2*y(ie+npts))
            ENDIF
         ENDDO
      ENDDO
 
c Set up the PLT commands
 
      nvec = ncomp - 1
 
      cmd(ncmd + 1) = 'LA T Current model'
      IF( ixopt.EQ.1 .OR. ixopt.EQ.2 )THEN
         cmd(ncmd + 2) = 'LA X Energy'
         IF ( nkeV .EQ. 2 )THEN
            cmd(ncmd + 3) = 'LA Y keV (keV/cm' // char(92) // 'u2' // 
     &                      char(92) // 'd s keV)'
         ELSEIF ( nkeV .EQ. 1 )THEN
            cmd(ncmd + 3) = 'LA Y keV/cm' // char(92) // 'u2' // 
     &                      char(92) // 'd s keV'
         ELSE
            cmd(ncmd + 3) = 'LA Y Photons/cm' // char(92) // 'u2' // 
     &                      char(92) // 'd s keV'
         ENDIF
      ELSEIF( ixopt.EQ.3 )THEN
         cmd(ncmd + 2) = 'LA X Angstroms'
         IF ( nkeV .EQ. 2 )THEN
            cmd(ncmd + 3) = 'LA Y A (A/cm' // char(92) // 'u2' 
     &                      // char(92) // 'd s A)'
         ELSEIF ( nkeV .EQ. 1 )THEN
            cmd(ncmd + 3) = 'LA Y A/cm' // char(92) // 'u2' // char(92)
     &                      // 'd s A'
         ELSE
            cmd(ncmd + 3) = 'LA Y Photons/cm' // char(92) // 'u2' // 
     &                      char(92) // 'd s A'
         ENDIF
      ENDIF
      cmd(ncmd + 4) = 'LOG X'
      cmd(ncmd + 5) = 'LOG Y'
      cmd(ncmd + 6) = 'ERR OFF 1'
      ncmd = ncmd + 6
      DO idset = 1, dgndst()
         WRITE(cmd(ncmd+1), '(''LINE STEP '', i2)')1 + (idset - 1)
     &         *(nvec - 1)
         ncmd = ncmd + 1
         WRITE(cmd(ncmd+1), '(''COLOR '', i4, '' on '', i4)')idset, 
     &         1 + (idset - 1)*(nvec - 1)
         ncmd = ncmd + 1
         DO icomp = 3, nvec
            WRITE(cmd(ncmd+1), '(''LST '',i2,'' on '',i2)')icomp - 1, 
     &            icomp + (idset - 1)*(nvec - 1) - 1
            WRITE(cmd(ncmd+2), '(''LINE STEP '',i2)')icomp + (idset - 1)
     &            *(nvec - 1) - 1
            WRITE(cmd(ncmd+3), '(''COLOR '', i4, '' on '', i4)')idset, 
     &            icomp + (idset - 1)*(nvec - 1) - 1
            ncmd = ncmd + 3
         ENDDO
      ENDDO
      iery(1) = 1
      DO icomp = 2, nvec
         iery(icomp) = 0
      ENDDO
      ymax = -1.E32
      ymin = 1.E32
      DO ie = 1, npts
         DO icomp = 3, nvec + 1
            IF( y(ie+(icomp-1)*npts).NE.NO )THEN
               ymax = max(ymax, y(ie+(icomp-1)*npts))
               ymin = min(ymin, y(ie+(icomp-1)*npts))
            ENDIF
         ENDDO
      ENDDO
      ymin = max(ymin, ymax/1.E3)
      ncmd = ncmd + 1
      cmd(ncmd) = 'R Y '
      WRITE(cmd(ncmd)(5:), *) ymin, ymax
      ncmd = ncmd + 1
      cmd(ncmd) = 'R X '
      WRITE(cmd(ncmd)(5:), *) xmin, xmax

c If required write the line ID labels.

      IF ( xpgqid() ) THEN
         labno = 0
         CALL mklnid(xmin, xmax, ymin, ymax, xpglkt(), xpglem(), 
     &               xpglz(), 1000, cmd, ncmd, labno, istat)
      ENDIF

      RETURN
      END
 
 
 
 
 
